Shallow foundations are typically considered as the simplest and most economical foundation for supporting small to medium size structures. They transfer the structural loads to the near surface soil that is mostly unsaturated and fluctuates with climatic condition. Recent studies show that the strength and deformation parameters of soil are influenced by the degree of saturation of the soil. Since the degree of saturation of the near surface soil varies with climatic and geotechnical parameters such as rainfall, water table depth, evapotranspiration, hydraulic conductivity, there can be variation in the strength and deformation parameters of the near surface soil. The two of the design considerations of shallow foundations are the safety against the overall shear failure in the soil and the settlement. This paper focuses on systematically incorporating the site specific climatic and geotechnical parameters in computing ultimate bearing capacity of soil.
The ultimate bearing capacity of a continuous shallow foundation is typically calculated using Terzaghi’s  ultimate bearing capacity equation (Equation (1)) assuming that the soil below the footing fails in general shear failure mode.
where is the cohesion of soil, is the unit weight of the soil, is the effective overburden pressure given by , is the depth from soil surface to bottom of footing, B is the width of footing, and are the dimensionless bearing capacity factors and are functions of the soil friction angle .
The application of Terzaghi’s equation to compute ultimate bearing capacity is limited because it is applicable to shallow footings (Df ≤ B) and subjected to concentric vertical loads. In 1963, Meyerhof  suggested a general bearing capacity equation to overcome the shortcomings of Terzaghi’s equation by introducing shape factors, depth factors, and load inclination factors. The general bearing capacity equation with the adjustment factors for shape, depth, and load inclination is shown in Equation (2).
where are shape, depth, and load inclination factors, respectively. Both Terzaghi’s and Meyerhof’s equations were derived for the failure mechanism and resistance along the failure surface based on saturated soil mechanics principles. However, recent studies show that the shear strength and volume change characteristics of soils vary with its degree of saturation  -  . It is also found that the shallow foundations designed based on the saturated soil mechanics principles are often conservative. Recently, the need for incorporating unsaturated soil mechanics principles in geotechnical engineering practice was summarized by Fredlund  with practical examples. In another related study, Mohamed and Vanapalli  showed that the ultimate bearing capacity of a square model footing on a coarse-grained soil under unsaturated conditions is approximately five to seven times higher than the bearing capacity under fully saturated condition. Therefore, taking into account the influence of the degree of saturation of the soil for computing the bearing capacity of shallow foundation in unsaturated soils can be economical at certain geotechnical, hydrological, and climatic conditions. Therefore, a comprehensive methodology must be developed based on unsaturated soil mechanics principles to incorporate the influence of the degree of saturation for computing ultimate bearing capacity. This paper presents such a methodology that considers site specific rainfall and water table depth in a probabilistic manner.
Unsaturated soil is a three phase medium. It consists of three bulk phases (solid, water, and air) and three interfaces (solid-water, water-air, and air-solid). Among the three interfaces, the air-water interface plays a critical role in the mechanical behavior of unsaturated soils    . The dynamic equilibrium among these bulk phases and interfaces can be expressed in terms of , , and  , where is the total stress, is the pore air pressure, and is the pore water pressure. Fredlund and Morgentern’s study  concluded that although any two of the three possible combinations can be used to describe the state of unsaturated soils, the net stress and the matric suction are the most satisfactory combinations for use in engineering practice. The combination is advantageous because the effects of a change in net stress can be separated from the effects of a change in pore water pressure. Also the combination is advantageous because pore air pressure, one of the two fluid pressures, is atmospheric in most engineering problems under normal loading conditions and therefore no need to calculate them. The matric suction is related to the degree of saturation of the soil through a constitutive relationship called Soil Water Characteristic Curve (SWCC). Therefore, knowing the value of degree of saturation, the matric suction can be computed and used in empirical or semi-empirical relationships for computing ultimate bearing capacity.
The shear strength parameters for a soil with matric suction are defined by Fredlund  as the effective angle of internal friction ( ), the effective cohesion ( ), and the angle of shear strength change with respect to matric suction ( ). The modified bearing capacity equation that considers the effect of suction is shown in Equation (3).
where is the matric suction.
Because of the difficulties in determining for use in Equation (3), researchers have proposed empirical and/or semi-empirical equations for computing ultimate bearing capacity of shallow foundations   . Among the many equations, the semi-empirical equation proposed by Vanapalli and Mohamed  is considered in this study. It should be noted that the framework shown in this study is not attached to a particular bearing capacity equation. Any other bearing capacity equation based on latest information can be used in the framework. The general form of the equation is shown in Equation (4). It should be noted that the Equation (4) consists of both degree of saturation and matric suction as variables in it in addition to the other variables found in the original bearing capacity equation.
where is the air entry value from SWCC, is the average air-entry value, is the effective friction angle, S is the degree of saturation, and ψ is the bearing capacity fitting parameter given by Equation (5).
where Ip is the plasticity index. The average suction in the above bearing capacity equation is given by Equation (6).
where is the matric suction at the bottom of the footing and is the matric suction at a depth equal to 1.5 times the width of the footing (1.5B).
The bearing capacity equation (Equation (4)) for unsaturated soil requires computing the matric suction profile within the influence zone of a shallow foundation. Since the matric suction is related to the degree of saturation through SWCC, one may easily compute matric suction from degree of saturation using SWCC. The degree of saturation in the soil can also be measured by soil sampling and dielectric sensor. Low resolution satellites such as NOAA-AVHRR and TERRA-MODIS can provide daily evapotranspiration fluxes in a clear sky at the 1 km scale  allowing for the determination of suction changes. Numerical models have also been developed for predicting the variation of matric suction and/or degree of saturation with depth for a give hydrological and geotechnical conditions   . Although these deterministic methods can be used to predict the spatial and temporal variation of suction, the application is limited because the rate of infiltration and water table depth varies randomly. Therefore, the problem must be approached in a probabilistic manner considering the extreme events and suitable return periods for the rainfall and water table depth.
This paper describes a probabilistic approach considering the randomness in the rainfall and water table depth using Monte Carlo method for computing design matric suction and degree of saturation within the influence zone of shallow foundation. Then, the design matric suction and degree of saturation are then used for computing the ultimate bearing capacity using Equation (4). In this study, the soil is assumed to be homogeneous with vertical 1-D downward flow of water. Also, this study evaluates the advantage of coupled hydrological-geotechnical approach by comparing the ultimate bearing capacities computed with conventional and proposed methods. A parametric study is also performed to investigate the sensitivity of the key parameters that influence the bearing capacity of the soil.
2. Analysis Procedure
The first step in performing a Monte Carlo simulation is defining the problem in terms of random variables and representing them with suitable probabilistic distributions. The variables in the bearing capacity equation arise from Equation (4). In this study, the base shear strength parameters of the soil are considered as constants. It should be noted that a comprehensive study may consider the shear strength parameters also as random variables. The inherent randomness of these variables can be incorporated into the Monte Carlo simulation technique but to allow for comparisons between the saturated and unsaturated soil bearing capacities the base shear strength parameters were kept as constants. The remaining random variable, i.e. the matric suction and degree of saturation must be defined in terms of its probability density function (PDF) to allow for the bearing capacity to be solved for through Monte Carlo simulations.
The spatial and temporal variations of matric suction are not something that has been recorded over long periods of time as a fundamental climatic data. It must be computed using the site specific recorded rainfall, evapotranspiration, surface runoff, water table location, and other climatic and geotechnical parameters as inputs for a mathematical model that represents the infiltration of water through initially unsaturated soil. Rainfall has been recorded in detail around the United States and many other countries. In the United States, the National Climatic Data Center (NCDC) records daily rainfall values and many sites have data for more than 60 years  . The daily rainfall data obtained from the NOAA (National Oceanic and Atmospheric Administration) and NCDC weather stations are utilized to derive precipitation frequency estimates in this paper. Precipitation frequency estimates are typically obtained by analyzing annual maximum series or partial duration series  . Annual maximum series were used in this study and are constructed by extracting the highest precipitation amount for a particular duration in each successive year of record. A water year starting on October 1 of the previous calendar year and ending on September 30 would be another typical option for selecting the maximum rainfall during a period of time. After the appropriate distribution is selected for the rainfall data, the range within the distribution must also be selected. If the entire distribution is used in the Monte Carlo simulation, the values randomly selected could have unrealistic return periods and may not represent realistic field conditions. In this paper, a range of return periods were tested, but in practice an engineer would have to select the return period appropriate for the project. The randomly selected rainfall events measured in inches are data inputs for the model. In order to determine the matric suction, the rainfall event needs to be quantified in terms of infiltration in units such as cm/hr. The worst case scenario was modeled with no runoff or pooling of water in this paper. However, runoff can be taken into account in the proposed procedure by quantifying the value and subtracting it from the total rainfall event. It was also assumed that the rainfall event would have an infiltration rate equal to the saturated hydraulic conductivity (Ks). The duration of the rainfall event was calculated by dividing the randomly selected rainfall event in inches by the saturated hydraulic conductivity of the soil. Dividing the rainfall event by Ks is a simplified approach compared to determining an infiltration rate with the unsaturated hydraulic conductivity K(h), which varies with time and depth. Besides rainfall, evapotranspiration and surface runoff affect the value of the influx at the ground surface which is the top boundary condition for solving the mathematical model that represents the flow of water through unsaturated soil (Richards equation). This will ultimately affect the temporal and spatial variation of matric suction in soil. In this paper, the surface runoff and evapotranspiration are ignored.
The water table depth is the bottom boundary condition which affects the matric suction and/or degree of saturation of soil in the influence zone of shallow foundation. The U.S. Geological Survey (USGS) and National Water Information System, provides data about the water table depths over a long period of time at various locations  . It is important to note that the water table data can only be determined in unconsolidated sand and gravel aquifers. There are over 26, 135 sites in these types of aquifers with 50 data points or more. The locations of these wells in the USA are shown in Figure 1. Fortunately, numerous wells are located in arid or semi-arid climates within the continental United States that are suitable for applying the proposed method for computing ultimate bearing capacity.
With rainfall and water table depth distributions, the spatial and temporal variations of matric suction and/or degree of saturation can be computed. The method used to calculate the matric suction is described in the next section. Besides the matric suction, the unit weight of the soil is another parameter that varies with changing degree of saturation of the soil. By modeling the flow of water into the soil considering infiltration and the water table depth as the top and bottom boundary conditions, the variation in unit weight of the soil can be weight values within the influence zone of the foundation, the Equation (4) can be solved for the number of Monte Carlo Simulations selected for the study.
Figure 1. Location of USGS wells in the continental USA.
Through numerous studies with different numbers of simulations, the location and shape parameters are checked for convergence. The converged data from the Monte Carlo simulation accurately represents the bearing capacity for the foundation. From the data collected in the Monte Carlo simulation, a cumulative distribution function (CDF) of the bearing capacity can be used to make risk assessments.
Design risk can be quantified by the product of probability of failure and the consequence of failure. But calculating the consequences of failure is difficult for practicing engineers, thus current practice in reliability based foundation design accounts for the consequence of failure indirectly by prescribing different target probabilities of failure  . The reliability index ( ), which is the inverse standard normal cumulative function of the probability of failure is the prescribed value usually published. The Canadian Building Code uses a target = 3.5 for the superstructure and the foundations. AASHTO uses a target = 3.5 for the superstructure and target values from 2.0 to 3.5 for the foundations. For this paper a probability of failure of 10−4 was selected to enable the bearing capacity to be read off the CDF. This is equivalent to a = 3.7, which exceed the Canadian Building Code and AASHTO target values. Figure 2 is a flow chart that outlines the method used to calculate the bearing capacity of unsaturated soil.
3. Calculation of SPATIAL Variation of Matric Suction and Unit Weight
3.1 .Flow of Water through Unsaturated Soil Zone and Spatial and Temporal Variations of Matric Suction and Degree of Saturation
Matric suction is directly related to the hydraulic head (hw) of the soil:
where ua is the atmospheric pressure, y is the gravitational head, and ρw is the density of water. The flow behavior of water in unsaturated soil is complex compared to the saturated soil because of the variation in hydraulic head with time and depth. The variation in hydraulic head with time and depth due to an infiltration event with the ground water table set at the datum can be solved using Richard’s equation in unsaturated soils.
where is the water capacity function.
The parameters in Richards equation can be solved for using the equations developed by van Genuchten  .
Figure 2. Flow chart of the analysis procedure.
where is the residual water content, is the saturated water content, is the approximation of the inverse of the pressure head at which the retention curve becomes the steepest, is the dimensionless water content, and n and m are model constants (typically ). All of these parameters are based on the soil type and are fitting parameters for an empirically determined soil water retention curve (SWRC).
After calculating these parameters, Richards equation can be solved numerically by the finite difference method. The Crank-Nicolson scheme implemented in Bolster and Raffensperger’s (1996) Matlab  program to solve Richards equation was implemented in the Monte Carlo simulation algorithm in this study. The result is the variation in hydraulic head, which as explained above can be used to solve for the term in the bearing capacity equation for the unsaturated soil.
The average air entry value, , is the other type of suction that needs to be solved. It is inversely proportional to the van Genuchten soil parameter α as given by the following equation.
where γw is the unit weight of water.
3.2. Spatial and Temporal Variations of Unit Weight
Since the unit weight is affected by the degree of saturation (S), the variation of unit weight is determined based on the average S in the influence zone of the foundation. In this study, the van Genuchten  soil water characteristic curve (Equation (12)) was used to solve for the moisture content of the soil with depth.
Using the same method used for calculating the , the average moisture in the soil within the influence depth of the footing (1.5B) can be calculated by:
Then, the average degree of saturation (S) in the influence area can be calculated by:
Finally, the variation in the unsaturated unit weight of the soil can be calculated using the average degree of saturation in the influence zone of the shallow foundation using the basic weight-volume relationship given by:
where e is the void ratio, Gs is the specific gravity of the soil solids, and w is the unit weight of water.
4. Sample Application
4.1. Problem Definition
A square 1.07 m × 1.07 m square footing located 0.61 m below the ground surface in Victorville, California was considered to demonstrate the proposed framework for computing ultimate bearing capacity considering site specific rainfall, water table and geotechnical data. This investigation was extended to determine if the depth and size of the footing on the computed ultimate bearing capacity. In addition, the influence of matric suction and degree of saturation related term and the depth factor shown in Equation (16) in the ultimate bearing capacity.
To determine if the size and depth have significant influence on the bearing capacity, three footings were considered and the computed results were compared. For the first example, the footing width (B) was increased from 1.07 to 1.52 m. Another example kept the ratio equal to the initial footing size and depth, thus B = 1.52 m and Df = 0.87 m. The last example allowed B to remain equal to 1.07 m and increase the Df to 0.87 m.
To determine how the return period of the hydrological parameters affects the bearing capacity of the footing, the square 1.07 m located 0.61 m below the ground surface was analyzed by considering return periods between: 1 yr and 2 yrs, 1 yr and 5 yrs, 1 yr and 10 yrs, 1 yr and 50 yrs, 1 yr and 100 yrs, 1 yr and 200 yrs and 1 yr and an infinite number of years for the rainfall and water table depth.
The Victorville, CA location was selected in this study due to its arid climate and the availability of van Genuchten soil water retention curve parameters for the Adelanto Loam located in this region. The van Genuchten parameters for the soil water retention curve of Adelanto Loam were taken from Zhang  are: θs = 0.423, θr = 0.158, α = 0.00321 cm−1, n = 1.26 and Ks = 0.003492 cm/min.
The soil strength parameters were taken from a geotechnical engineering report by Kleinfelder  . The report was from a site about 15 miles from Victorville, CA. The site was a similar distance from the river that passes Victorville, thus it was assumed that the water table would be reasonably similar. The dry unit weight for the soil at a depth of 1.52 m is 16.19 kN/m3. The angle of internal friction for the soil at a depth of 1.52 m is 33 degrees. The cohesion at a depth of 1.52 m is 0. The USCS soil type for the soil at 1.52 m is SM.
4.2. Site Specific Rainfall Data
The rainfall data was taken from the Victorville Pump Station, Victorville, CA, within the climate division CA-07. The station was in service from November 1, 1938 to the present. The elevation of the station is 871 m (2858 ft) above mean sea level. The latitude and longitude of the station are 34˚32'00''N and 117˚17'34''W, respectively. The data for the pump station was processed from an ASCII file that was downloaded from the National Climatic Data Center  . The maximum rainfall in inches during a year was tabulated for each year 1938-2009. Years where not all 365 days were recorded were removed from the data set. This prevents non-rainy season maximum yearly values from affecting the overall distribution. Out of 72 years, a total of 6 years was excluded from the data set. To determine the best fitting distribution for the rainfall data, the probability paper plotting technique was used. The Type II Extreme Largest (Frechet distribution), Type I Extreme Largest (Gumbel distribution), and the Type III Extreme Largest (Weibull distribution) were checked for the best fit. The Gumbel distribution was deemed the best fit based on R2 values. The probability plot of rainfall data for the Gumbel distribution is shown in Figure 3(a).
Using the Gumbel distribution the CDF was transformed into a linear equation shown below, it can be determined that the location parameter, and the shape parameter .
where xi is the annual maximum rainfall data and n is the number of data points.
4.3. Site Specific Water Table Data
The water table data at Victorville, CA was taken from the U.S. Geological Survey National Water Information System: Web Interface  . The water table depth was recorded between the years of 1930 and 1958. To determine the best-fit distribution for the water table data the probability paper plotting technique was used. For this case, the Frechet Distribution (Type II Extreme Largest), Gumbel Distribution (Type I Extreme Largest), Weibull Distribution (Type III Extreme Largest), Normal distribution, and Lognormal distribution were checked for the best fit. The Frechet distribution had the best fit, but for simplicity the Gumbel distribution, the second best fit was used. The probability plot of the water table data for the Gumbel distribution is shown in Figure 3(b). Using Equation (17) the Gumbel distribution parameters were determined as and .
4.4. Results and Discussions
The convergence of the mean and the coefficient of variation for the bearing capacity distributions with the number of simulations are plotted in Figure 4. At 10,000 simulations, there is evidence of a convergence for each of the different example footings. The mean is the location parameter and the coefficient of variation takes into account the shape. With both of these measurements of the distribution, it can be understood that the empirical CDF created from the Monte Carlo simulations accurately represents the bearing capacity for the footings considered in Victorville, CA.
The CDFs created for each of the example footings using 10,000 simulations are plotted in Figure 5. The bearing capacity for the footing at a probability of failure of 10−4 is tabulated in Table 1. To determine if considering the unsaturated soil mechanic principles is advantages for computing the ultimate bearing capacity over the conventional method based on fully saturated soil mechanics principles, the Meyerhof’s equation was used to calculate the bearing capacity of the footing assuming the soil to be fully saturated (Table 1) and compared with
Figure 3. Type I Extreme Largest (Gumbel distribution) for rainfall data and water table data for Victorville, CA location.
Figure 4. Number of simulations versus the mean and coefficient of variation.
the bearing capacity calculated with unsaturated soil mechanics principles. The percent increases in bearing capacity from conventional method to the proposed method are also listed in Table 1 for various footing sizes and depths.
From the results in Table 1, it is evident that the depth factor in the cohesion term has a significant influence in the ultimate bearing capacity of shallow
Figure 5. Bearing capacity versus probability of failure for all systems in Victorville (a).
Table 1. Computed ultimate bearing capacities.
foundation. The footing with the 1.07 m width and the depth of 0.87 m (case 2) showed the highest bearing capacity increase. This is a larger bearing capacity than the larger footing at the same depth. This shows that a smaller depth factor has more influence than a larger footing. It is clearly evident that the bearing capacity of the soil is significantly affected by the matric suction and the variation of unit weight. All the bearing capacities for the different example footings have increased by over 233% when using the proposed method based on unsaturated soil mechanics principles.
The effect of the return period (RP) on the bearing capacity for the case 1 footing is plotted in Figure 6(a). As expected, the decrease in the return period resulted in increases in the ultimate bearing capacity of the footing with a probability of failure of 10−4. Figure 6(b) provides a plot that an engineer could use to select a probability of failure for a foundation, the return period for the rain and water table event, and the bearing capacity for the site.
5. Sensitivity Analyses
5.1. Sensitive Parameters and Analysis Procedures
Sensitivity analyses were conducted to determine the effect of variations in input parameters used in the calculations and also served as additional examples to
Figure 6. Bearing capacity versus probability of failure for 1.06 m square footing at a depth of 0.60 m for different return periods in Victorville (a).
explain the procedure. The first sensitivity analysis was to determine if another site with increased rainfall events would significantly decrease the increased bearing capacity determined from the sample application. The Levelland, Texas was selected as the second site in this study. The required data for the Levelland site was collected following the procedure for the Victorville site. The soil strength parameters were obtained from a geotechnical report made by Amarillo Testing and Engineering, Inc  . The soil parameters are summarized in Table 2. The Gumbel distribution was once again the best fit for the rainfall data in Levelland (Figure 7(a)) and the normal distribution was the best fit for the water table data in Levelland (Figure 7(b)). Compared to the sample application which used van Genuchten parameters published  for a specific soil, a new approach for defining the van Genuchten parameters is used in this sensitivity analysis.
The van Genuchten values determined through the lowest of the hierarchical sequences in the program Rosetta  are printed in a table in the program’s help index. The values in that table were generated by computing the average value for each textural class. Each textural class includes one standard deviation of uncertainty for the van Genuchten parameters. The soil classification in the geotechnical report was used to determine the best class of Levelland Texas soil. The Levelland was considered to be in the sandy clay textural class. The van Genuchten parameters for this soil are in Table 2.
Through this new method of collecting van Genuchten parameters, a second sensitivity analysis was performed. From the sample application it is impossible to determine if the ultimate bearing capacity of footing is controlled by the rainfall and water table distributions or the van Genuchten. In order to compare the effect of van Genuchten parameters, the sample application results which use van Genuchten parameters specifically for the Adelanto Loam are compared to the van Genuchten parameters recorded in the textural class determined from the soil classification in the geotechnical report. The soil in the Victorville site was considered as sandy loam. The van Genuchten parameters for Victorville are also listed in Table 2. The results from the sample application at the Victorville
Figure 7. Type I Extreme Largest (Gumbel distribution) for rainfall data and normal distribution for water table data.
Table 2. USDA Textural Class Average Values of Hydraulic Parameters and Soil Parameters for Victorville and Levelland Sites at Depth of 1.524 m
site are referred to as Victorville (Adelanto Loam) while the results from the van Genuchten parameters from the Rosetta  table are Victorville (Sandy Loam).
A third sensitivity analysis was also performed by testing each van Genuchten parameter determined by the textural class individually. Each parameter was increased by one standard deviation to determine which parameter plays the largest role in affecting bearing capacity of shallow footings. Each parameter’s one standard deviation increase for both the Victorville (Sandy Loam) and the Levelland is summarized in Table 2. This analysis also allows for a discussion of the importance of accurate SWRC in the ultimate bearing capacity calculation using the proposed method.
The final sensitivity analysis observed the influence of the depth factor in the cohesion term on the bearing capacity. Once again the footing was tested with the four different variations of footing size and depth described in the sample application for both the mean values of Levelland and Victorville (Sandy Loam).
5.2. Results and Discussions
The results from the first sensitivity analysis, a change in location are recorded in Table 3. The CDF for the simulation is plotted in Figure 8. It is evident that the bearing capacity of the footing in Levelland, Texas is lower than the bearing capacity in Victorville, California. To determine if the decrease in bearing capacity is primarily due to a change in the soil strength parameters comparisons were made between the saturated and unsaturated soil bearing capacities. Even with increased rainfall, the soil in Levelland had an increase in bearing capacity over 76% (Table 3) through the Monte Carlo simulation method. Thus this method is valuable to use in regions that are not one of the driest in the United States.
The sensitivity analysis of the van Genuchten parameters was tested by comparing the bearing capacities of the Victorville (Adelanto Loam) and the Victorville (Sandy Loam) sites. The bearing capacity for the case 1 footing for the Victorville
Table 3. Computed bearing capacities for Victorville and Levelland with sensitivity analyses.
*Victorville (Sandy Loam).
Figure 8. Bearing capacity versus probability of failure for all systems in Levelland.
(Adelanto Loam) is 1548 KN/m2. The bearing capacity for the same footing for the Victorville (Sandy Loam) is 872 KN/m2. All the footing sizes and depth have similar changes (refer to Table 1 and Table 3). The changes in parameters also affected the distribution of bearing capacities computed form the Monte Carlo simulation. Comparing Figure 4 to Figure 9 it is noticeable that the CDFs are much steeper. Though the rainfall and water table distributions are very important, this sensitivity analysis revealed that the SWRC parameters are equally important in accurately predicting the increase in bearing capacity in this probabilistic method.
Based on the individual sensitivity analysis of each parameter, the parameter alpha in Equation (9) resulted in the highest change in bearing capacity while increasing it by one standard deviation (refer to Table 3). It is reasonable that this parameter has a significant control since it describes the matric suction a soil has when it is almost completely saturated. The smaller the alpha value is; the greater the matric suction is at higher degrees of saturation. The rest of the van Genuchten parameters have similar sensitivities to the increase in bearing capacity when increased by one standard deviation (refer to Table 3). With this probabilistic method it is most important to accurately predict alpha, while fitting a SWRC to the volumetric moisture content and hydraulic head data.
The influence of the depth factor was studied by comparing the percent increase in the ultimate bearing capacity for two footings of different sizes at a specified depth. The difference between the percent increase for the deterministic method with fully saturated soil and the Monte Carlo simulation with unsaturated soil was calculated. Since the depth is the same, the effects of matric suction are constant. Considering this, the percent increase in the ultimate bearing capacity for both cases should only be due to the change in size of the footing. Thus the percent difference between the two methods for increased bearing capacity should be the same; this is not the case in Table 4. The important trend in Table 4 is that as the influence of matric suction increases in the bearing
Figure 9. Bearing capacity versus probability of failure for all systems in Victorville (Sandy Loam).
Table 4. Effects of depth factor on the cohesion term in the bearing capacity equation.
capacity, the influence of the depth factor increases. The difference in percent increase in bearing capacity calculated using the two methods while increasing the footing size reduces by 1.70% for the depth of 0.87 m in Victorville (Sandy Loam). For sites where suction has more influence such as Victorville (Adelanto Loam), the difference in the percent increase in bearing capacity due to an increase in footing size reduces by 5.18% for a footing at a depth of 0.87 m when comparing the two methods. The negative percent increase is due to the depth factor having more control than the increased footing size in Victorville (Adelanto Loam), which is explained in the results section of the sample application section. From these results it is evident that the depth factor has influence on the bearing capacity calculated from the bearing capacity equation.
The method for determining the bearing capacity of a footing in unsaturated soil using Monte Carlo simulation is a useful tool for further understanding how unsaturated soil mechanics can be applied in problems faced by practicing engineers. The sample study gave evidence that considering unsaturated soils in design can increase the bearing capacity of a footing by at least 2.3 times the bearing capacity calculated using Meyerhof’s equation. This study considered homogeneous soil with 1-D flow. However, the proposed method can be easily extended for 2-D and 3-D cases.
The sensitivity analysis reinforced the importance of having accurate SWRC or SWCC parameters when using methods relying on the SWRC or SWCC. The sensitivity analysis also confirmed that sites with additional rainfall can still benefit from an increase in bearing capacity by considering unsaturated soils. In the case of Levelland, Texas more than 76% increase in bearing capacity was observed.
The effect of the depth factor is an important finding from the sensitivity analysis. As the suction increases the value of the cohesion term, the depth factor has a greater influence on the ultimate bearing capacity. This results in smaller factors of increase in bearing capacity when there is an increase in footing size which creates a conservative estimate of the bearing capacity of footings in high matric suction.
 Steensen-Bach, J.O., Foged, N. and Steenfelt. J.S. (1987) Capillary Induced Stresses— Fact or Fiction? Proceedings of the 9th European Conference on Soil Mechanics and Foundation Engineering, Budapest, 83-89.
 Costa, Y.D., Cintra, J.C. and Zornberg, J.G. (2003) Influence of Matric Suction on the Results of Plate Load Tests Performed on a Lateritic Soil Deposit. Geotechnical Testing Journal, 26, 219-226.
 Rojas, J.C., Salinas, L.M. and Seja, C. (2007) Plate-Load Tests on an Unsaturated Lean Clay. Experimental Unsaturated Soil Mechanics. Proceedings of the 2nd International Conference on Unsaturated Soils, Weimar, 445-452. https://doi.org/10.1007/3-540-69873-6_44
 Vanapalli, S.K. and Mohamed, F.M.O. (2007) Bearing Capacity of Model Footings in Unsaturated Soil. In: Schanz, T., Ed., Experimental Unsaturated Soil Mechanics, Springer-Verlag, Berlin, Heidelberg, 483-493. https://doi.org/10.1007/3-540-69873-6_48
 Fredlund, D.G. (2006) Unsaturated Soil Mechanics in Engineering Practice. Journal of Geotechnical and Geoenvironmental Engineering, 132, 286-321.
 Mohamed, F.O. and Vanapalli, S.K. (2006) Laboratory Investigations for the Measurement of the Bearing Capacity of an Unsaturated Coarse-Grained Soil. Proceedings of the 59th Canadian Geotechnical Conference, Vancouver, 219-226.
 Vanapalli, S.K., Fredlund, D.G., Pufahal, D.E. and Clifton, A.W. (1996) Model for the Prediction of Shear Strength with Respect to Soil Suction. Canadian Geotechnical Journal, 33, 379-392. https://doi.org/10.1139/t96-060
 National Climatic Data Center. U.S. Department of Commerce. NOAA Satellite and Information Service. National Environmental Satellite, Data, and Information Service.
 Perica, S., Martin, D., Lin, B., Parzybok, T., Riley, D. and Yekta, M. (2011) Precipitation-Frequency of the United States. NOAA Atlas, Vol. 4.
 U.S. Geological Survey. National Water Information System: Web Interface.
 Kulhawy, F.H., Phoon, K.K. and Wang, Y. (2012) Reliability-Based Design of Foundations—A Modern View. Geotechnical Engineering State of the Art and Practice. Geo Congress Oakland California, 102-121.
 Van Genuchten, M.T. (1980) A Closed Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Science Society of America Journal, 44, 892-898.