Analysis of the Bearing Capacity of Shallow Foundation in Unsaturated Soil Using Monte Carlo Simulation

The ultimate bearing capacity of shallow foundation supported by unsaturated soil depends on the degree of saturation of the soil within the influence zone because the strength and deformation parameters of soil are affected by the degree of saturation. As the degree of saturation varies with rainfall, surface runoff, evapotranspiration and other climatic and geotechnical parameters, these parameters must be systematically incorporated for accurately computing the ultimate bearing capacity. In this study, a framework is proposed to compute the ultimate bearing capacity of a shallow footing in unsaturated soil considering site specific rainfall and water table depth distributions. The randomness in rainfall and water table depth is systematically considered using Monte Carlo method. The infiltration of water through the unsaturated zone is modelled using Richards equation considering infiltration and water table location as the top and bottom boundary conditions, respectively. The results show that the bearing capacity calculated using the proposed method is approximately 2.7 times higher than that calculated using the deterministic approach with fully saturated soil parameters.


Introduction
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 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 [1] ultimate bearing capacity equation (Equation ( 1)) assuming that the soil below the footing fails in general shear failure mode.
where c′ is the cohesion of soil, γ is the unit weight of the soil, q is the ef- fective overburden pressure given by ( ) ( ) ( ) where , , and F F F are shape, depth, and load inclination factors, respective- ly.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 [3]- [8].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 [9] with practical examples.In another related study, Mohamed and Vanapalli [10] 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 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 [8] [11] [12].The dynamic equilibrium among these bulk phases and interfaces can be expressed in terms of ( ) , and ( ) , where σ is the total stress, a u is the pore air pressure, and w u is the pore water pressure.Fredlund and Morgentern's study [11] concluded that although any two of the three possible combinations can be used to describe the state of unsaturated soils, the net stress ( ) The shear strength parameters for a soil with matric suction are defined by Fredlund [12] as the effective angle of internal friction ( φ′ ), the effective cohe- sion ( c′ ), and the angle of shear strength change with respect to matric suction ( b φ ).The modified bearing capacity equation that considers the effect of suction is shown in Equation (3).
( ) where ( ) a w u u − is the matric suction.Because of the difficulties in determining b φ for use in Equation (3), re- searchers have proposed empirical and/or semi-empirical equations for computing ultimate bearing capacity of shallow foundations [7] [13].Among the many equations, the semi-empirical equation proposed by Vanapalli and Mohamed [7] 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, ( ) a w ave u u − 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).
( ) ( ) 1.0 0.34 0.0031 where I p is the plasticity index.The average suction in the above bearing capacity equation is given by Equation ( 6).
where ( ) 1 a w u u − is the matric suction at the bottom of the footing and a w u u − 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 [14] 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 [15] [16].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.

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 [17].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 [18].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   probabilities of failure [20].The reliability index ( β ), which is the inverse stan- dard 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 Cana- dian 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.

.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 (h w ) of the soil: where u a 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 dh is the water capacity function.
The parameters in Richards equation can be solved for using the equations developed by van Genuchten [21].
Figure 2. Flow chart of the analysis procedure.
( ) ( ) where r θ is the residual water content, s θ 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 1 1 m n = − ).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 ( ) where γ w is the unit weight of water.

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 [21] 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 ( ) ( ) ( ) 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, G s is the specific gravity of the soil solids, and w is the unit weight of water.

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.
1 0.4 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 The soil strength parameters were taken from a geotechnical engineering report by Kleinfelder [24].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/m 3 .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.

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 [17].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 R 2 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, 0.8472 and the shape parameter 0.5011 where x i is the annual maximum rainfall data and n is the number of data points.

Site Specific Water Table Data
The

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 International Journal of Geosciences  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   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.

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  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 [25].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 [23] 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 [26]  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 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.

N.
Ravichandran et al.DOI: 10.4236/ijg.2017.8100711232 International Journal of Geosciences 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.

D
is the depth from soil surface to bottom of footing, B is the width of footing, and , , and c q N N N γ 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 (D f ≤ B) and subjected to concentric vertical loads.In 1963, Meyerhof [2] 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).
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.
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 .
Figure 1.Location of USGS wells in the continental USA.
soil within the influence depth of the footing (1.5B) can be calculated by:

1 .
52 m.Another example kept the f D B ratio equal to the initial footing size and depth, thus B = 1.52 m and D f = 0.87 m.The last example allowed B to remain equal to 1.07 m and increase the D f 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[23] are: θ s = 0.423, θ r = 0.158, α = 0.00321 cm −1 , n = 1.26 and K s = 0.003492 cm/min.

Figure 3 .
Figure 3. Type I Extreme Largest (Gumbel distribution) for rainfall data and water table data for Victorville, CA location.

Figure 4 .
Figure 4. Number of simulations versus the mean and coefficient of variation.
m D f =0.60 m International Journal of Geosciences

Figure 5 .
Figure 5. Bearing capacity versus probability of failure for all systems in Victorville (a).
m Df=0.60 m International Journal of Geosciences

Figure 6 .
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).
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.
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.
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 Infor- N.Ravichandranet al. DOI: 10.4236/ijg.2017.8100711236 International Journal of Geosciences (Ks).mation System, provides data about the water table depths over a long period of time at various locations [19].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 Solve for the Bearing Capacity of a Footing in Unsaturated SoilsDetermine Location, Size, and Depth of the bottom of the Footing water table data at Victorville, CA was taken from the U.S. Geological Sur- [19]National Water Information System: Web Interface[19].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 Figure3(b).Using n β = .

Table 2 .
[26]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 Table2.USDA Textural Class Average Values of Hydraulic Parameters and Soil Parameters for Victorville and Levelland Sites at Depth of 1.524 m are referred to as Victorville (Adelanto Loam) while the results from the van Genuchten parameters from the Rosetta[26]table are Victorville (Sandy Loam).
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 N. Ravichandran et al.DOI: 10.4236/ijg.2017.8100711245 International Journal of Geosciences site