Soils and Rocks Analytical and numerical methods for prediction of the bearing capacity of shallow foundations in unsaturated soils

for prediction or estimation of the bearing capacity of unsaturated soils. These studies clearly show that the bearing capacity of unsaturated soils Abstract Bearing capacity of saturated soils can be estimated using effective or total stress approaches extending the concepts proposed by Terzaghi (1943) and Skempton (1948), respectively. Recent studies have shown that similar approaches (i.e., Modified Effective Stress Approach, MESA and Modified Total Stress Approach, MTSA) can be used for interpretation and prediction of the bearing capacity of unsaturated soils by considering the influence of matric suction. However, comprehensive discussion for the application of the MESA and the MTSA in geotechnical engineering practice applications is lacking in the literature. For this reason, in this state-of-the-art paper, the background associated with the MESA and MTSA is first introduced. The analytical and numerical methods available for the prediction of the bearing capacity of unsaturated soils from the literature are revisited. The various available methods are explained by categorizing them into two groups: MESA and MTSA along with their applications using examples. The focus of this state-of-the-art paper is directed towards not only for providing tools for rational understanding but also for better prediction of the bearing capacity of unsaturated soils for extending them in geotechnical engineering practice applications.


Introduction
The soil bearing capacity is key information required in the design of shallow and deep foundations. In conventional engineering practice, bearing capacity of soils in many scenarios is estimated using either analytical or numerical methods. Both these methods require saturated shear strength parameters of soil taking account of their rate of loading and drainage conditions. For example, the bearing capacity equation originally proposed by Terzaghi (1943) requires the effective shear strength parameters, c' and ϕ' (i.e., Effective Stress Approach, ESA) assuming drained condition in terms of pore-water. On the other hand, Skempton (1948) suggested that ϕ u = 0 be used in calculating the bearing capacity of saturated soils under undrained condition (i.e., Total Stress Approach, TSA). In other words, the saturated shear strength parameters must be reliably determined taking account of rate of loading and drainage conditions to minimize the uncertainties associated with the calculation or prediction of the bearing capacity.
Research studies related to the bearing capacity of unsaturated soils highlight that it is significantly different from saturated soils (Steensen-Bach et al., 1987;Oloo et al., 1997;Costa et al., 2003;Mohamed & Vanapalli, 2006;Rojas et al., 2007;Oh & Vanapalli, 2013;Tan et al., 2021). These studies also suggest that it is not reliable to use conventional approaches for estimating the bearing capacity of unsaturated soils. A third of the earth's surface constitutes of arid or semi-arid regions (Fredlund & Rahardjo, 1993). The natural ground water table is relatively at a greater depth from the natural ground surface in these regions. Due to this reason, the base of shallow foundation or major portion of a deep foundation is typically located in the vadose zone where soils are in a state of unsaturated condition. The influence of matric suction (or negative pore-water pressure) should be taken into account for the determination of the bearing capacity of unsaturated soils. Several researchers have developed analytical (e.g., Oloo et al., 1997;Vanapalli & Mohamed, 2007;Oh & Vanapalli, 2013;Vahedifard & Robison, 2016;Vo & Russel, can be estimated by extending the conventional ESA and TSA considering the influence of matric suction, which are referred to as the Modified Effective Stress Approach (i.e., MESA) and the Modified Total Stress Approach (i.e., MTSA).
In this state-of-the-art paper, the background of MESA and MTSA is first discussed in relation to the shear strength of unsaturated soils. The authors then revisit the analytical and numerical methods that are available in the literature to predict the bearing capacity of unsaturated soils. These methods are succinctly explained by categorizing them into two groups; namely, MESA and MTSA highlighting them with examples. In addition, the MESA proposed by Vanapalli & Mohamed (2007) is improved to provide better prediction of the bearing capacity of unsaturated coarse-grained soils for all the three zones; namely boundary effect, transition and residual zones that can be derived from the Soil-Water Characteristic Curve (SWCC). The summarized information in this state-of-the-art paper is useful for the practicing engineers for estimating the bearing capacity of unsaturated soils using limited information, which include the saturated shear strength properties and the SWCC.

Shear strength of unsaturated soils
The engineering behavior of saturated soils has been successfully interpreted and implemented in practice using the effective stress, σ' as a tool. Effective stress is an independent stress state variable that is defined as the difference between total stress, σ and pore-water pressure, u w (i.e., σ' = σ -u w ; Terzaghi, 1936). The shear strength parameters in terms of effective stress can be determined by combining the effective stress with the Mohr-Coulomb failure criteria. The effective shear strength parameters have been widely used in conventional geotechnical engineering practice such as the design of foundations/retaining walls and stability analysis of slope/excavation by ignoring the influence of matric suction. Such an approach has been used in practice based on the assumption that ignoring the influence of matric suction leads to a conservative design. However, various case studies showed that stability analysis of natural soil slopes (Hong Kong Government, 1972) or unsupported cuts (Richard et al., 2021) conducted by ignoring the failure mechanism attributed to the influence of matric suction can result in life losses. Bishop (1959) extended effective stress equation for saturated soils proposed by Terzaghi (1943) to unsaturated soils by introducing a parameter, χ along with two stress state variables. The shear strength of unsaturated soils can be estimated extending Bishop (1959) where τ us is shear strength of unsaturated soil, c' is effective cohesion, ϕ' is effective internal friction angle, σ' us is effective stress of unsaturated soil, (σ -u a ) is net normal stress, (u a -u w ) is matric suction, u a is pore-air pressure, u w is pore-water pressure, and χ is parameter that is a function of degree of saturation Research studies that followed after the pioneering work of Bishop (1959) suggested that the net normal stress, (σ -u a ) and matric suction, (u a -u w ) should be considered as independent stress state variables for rational interpretation of the engineering behavior of unsaturated soils. In saturated soils, there are identical changes in water content and void ratio associated with changes in (σ -u w ) because they are uniquely related; however, the mechanical behavior of unsaturated soils associated with similar changes in the (σ -u a ) and (u a -u w ) are not identical and hence their influence cannot be rationally interpreted using a single state effective stress equation for unsaturated soils (Bishop & Blight, 1963). Fredlund et al. (1978) analyzed the measured shear strength of unsaturated soils available in the literature and validated the 'two independent stress state variables' approach using Equation 2.
where ϕ B is angle indicating the rate of increase in shear strength relative to the matric suction For the suction values less than the air-entry value, shear strength contribution due to matric suction, ϕ B = ϕ'; therefore, shear strength increases linearly with increasing matric suction. Beyond the air-entry value, ϕ B contribution gradually decreases with an increase in matric suction; due to this reason, ϕ B is less than ϕ'; therefore, shear strength of unsaturated soils varies nonlinearly with respect to suction (Escario & Sáez, 1986, Gan et al., 1988. The typical shear strength behavior of unsaturated soils for conventional soils such as the sand, silt, glacial till and low plastic clays is illustrated in Figure 1. As can be seen, the shear strength can be either overestimated or underestimated when a constant ϕ B is used depending on the range of suction and soil type. Vanapalli (2009) provides a comprehensive summary of various models that are available in the literature to estimate or predict the nonlinear variation of shear strength of unsaturated soils. Many of the shear strength models that were proposed in the literature use the Soil-Water Characteristic Curve (SWCC), which is the relationship between degree of saturation (alternatively volumetric water content or water content) and soil suction. The matric suction in a soil specimen at failure, however, can be different from the initial matric suction due to the influence of various parameters such as volume change, stress history, shearing rate and drainage condition of poreair and pore-water. This discrepancy of matric suction can contribute to some differences between the measured and estimated or predicted shear strength of unsaturated soils. Lu et al. (2010) suggested a different approach which is consistent with conventional soil mechanics principles. Their approach is based on single stress state variable (i.e., effective stress) supporting the use of Equation 3 that takes the same form as Terzaghi's effective stress by using the suction stress, σ S .
where S is degree of saturation, S e [= (S -S r )/(1 -S r )] is effective degree of saturation, S r is residual degree of saturation A closed-form expression for suction stress, σ S is shown in Equation 4 by eliminating the degree of saturation using van Genuchten (1980) SWCC model.
where n and α are van Genuchten (1980) empirical fitting parameters of unsaturated soil properties.

Modified effective stress approach
Various models have been proposed to estimate the bearing capacity of unsaturated soils (see Table 1). Oloo et al. (1997) proposed a model that can be used to estimate the bearing capacity of a shallow foundation in unsaturated soils considering the influence of matric suction (Equation 5). The model was developed extending the Terzaghi (1943)'s effective stress approach assuming both pore-air and porewater are in drained condition during the loading stage. This approach was referred to as the Modified Effective Stress Approach (MESA) by Oh & Vanapalli (2013). Oloo et al. (1997)'s model uses a constant ϕ B for the suction values greater than the air-entry value. However, this approach can either over-or underestimate bearing capacity depending on the range of suction. Such a behavior is consistent with the reason explained using Figure 1 for the shear strength of unsaturated soils. Vanapalli & Mohamed (2007) improved the Oloo et al. (1997)'s model by adopting the nonlinear shear strength models proposed by Vanapalli et al. (1996) and Fredlund et al. (1996) (Equation 6) (hereafter referred to as VM model). Terzaghi (1943)'s original bearing capacity equation was developed based on a plane strain condition for continuous footings. Several research studies suggest that ϕ' obtained for plane strain condition is 1.1 times greater than that obtained from a conventional triaxial test under axisymmetric condition. The same approach was also used by Danish code of practice DS 415 (DSCE 1984) and Steensen-Bach et al. (1987). Vanapalli & Mohamed (2007) also suggested that 1.1ϕ' can provide better prediction of bearing capacity for both saturated and unsaturated conditions. Table 2 summarizes the bearing capacity and shape factors used by Vanapalli & Mohamed (2007) and other researchers (Terzaghi, 1943;Meyerhof, 1963;Vesić, 1973). The comparison between the measured and predicted bearing capacity values using the model footing test results for saturated condition (Steensen-Bach et al., 1987;Mohamed & Vanapalli, 2006) showed better agreement when predicted using the bearing capacity and shape factors proposed by Vanapalli & Mohamed (2007), as shown in Figure 2.
Vahedifard & Robinson (2016) proposed a unified method for estimating the bearing capacity of unsaturated soils considering the variation of degree of saturation associated with steady state flow conditions (Equation 7). The contribution of matric suction on the bearing capacity was taken into account by introducing the average suction stress, which is the matric suction at the centroid of matric suction profile in a certain depth (i.e., from the base of shallow foundation to the depth of 1.5B or 2B, Oh & Vanapalli, 2013). The effective degree of saturation under constant infiltration and evaporation rate, q is calculated using Equation 8 (Yeh, 1989;Lu & Likos, 2004). Vo & Russel (2016) proposed charts that can be used to estimate the bearing capacity of smooth and rough footings in unsaturated soils (Equation 9) extending the research by Martin (2004). Soil suction at the depth, z under constant surface infiltration, q is estimated using Equation 10. Ghasemzadeh & Akbari (2019) proposed a method to predict the bearing capacity of footings placed on unsaturated soil extending the limit equilibrium method. The bearing capacity equation consists of two terms; bearing capacity for saturated condition and contribution of matric suction towards the bearing capacity due to uniform and nonlinear matric suction distribution. Garakani et al. (2020) proposed an analytical solution for the bearing capacity of shallow foundation in unsaturated soil (Equation 14). The influence of suction was taken into account by adopting a correction factor C f that depends on the soil Ghasemzadeh & Akbari (2019) Constant suction distribution (u a -u w ) AVE : average matric suction Garakani et al. (2020) ( ) N ct , N qt , N γt : bearing capacity factors based on the unified strength theory (N s ) cst , (N s ) l : bearing capacity factor in terms of suction for constant and linear suction distribution, respectively q: constant infiltration at the ground surface q s : surface surcharge s: soil suction F, V: dimensionless ratios z': height above a horizontal ground water table α and n: van Genuchten (1980) fitting parameters λ m : rate of decrease of matric suction with depth γ: unit weight of soil γ m : modified average unit weight ψ: fitting parameter for bearing capacity ξ c , ξ q , ξ γ : shape factors ξ s : shape factor for suction (=ξ c ) ϕ': effective internal friction angle ϕ t ': unified internal friction angle ϕ t B : unified internal friction angle due to the contribution of suction ϕ B : internal friction angle due to the contribution of suction ϕ R : internal friction angle due to the contribution of suction in residual zone Table 2. Summary of the bearing capacity and shape factors used by Vanapalli & Mohamed (2007) and other researchers. Vesić (1973) Vesić (1973) K pγ = passive earth pressure used by Terzaghi (1943) to calculate N γ (differs from the passive horizontal stress coefficient defined by Rankine's limit state); P γmin = minimum passive pressure used in fining the critical failure surface by Kumbhokjar (1993). The bearing capacity models listed in Table 1 were developed based on the failure surface under a continuous footing following the approach by Prandtl (1921) and Terzaghi (1943). The contribution of matric suction was considered as a form of total cohesion or bearing capacity factor. The change in the mechanical properties of unsaturated soils can be predominantly attributed to the environmental factors such as the rainfall. The weather-imposed conditions primarily affect the matric suction component rather than the osmotic suction. Leong & Abuel-Naga (2018) conducted unconfined compression tests on compacted soil specimens. The specimens were prepared using distilled water and sodium chloride solutions to investigate the influence of matric suction and osmotic suction on the shear strength, respectively. They concluded that the influence of osmotic suction on the shear strength of compacted specimens is negligible. According to Fredlund & Rahardjo (1993), the influence of osmotic suction on the soil behavior may be significant only when the salt content is altered by chemical contamination or chemical change. For such a scenario, it is necessary to consider the osmotic suction as part of stress state or as an independent stress state variable in estimating the bearing capacity of unsaturated soils.
In this state-of-the-art-paper, the bearing capacity model proposed by Vanapalli & Mohamed (2007, VM model) is revisited. First, the relationship between the shear strength behavior and the SWCC of sands are highlighted using Donald (1957) results, who performed direct shear tests for four different sands at different matric suction values. The variation in the bearing capacity of sands with respect to matric suction would be similar to that of shear strength of the respective unsaturated sands. The results in Figure 3 highlight that the shear strength increases with an increase in the matric suction up to a certain value and then start decreasing with a further increase in matric suction. This is attributed to the reason that the contribution of matric suction towards shear strength decreases as matric suction approaches the residual suction value. This characteristic behaviour was also explained extending mathematical framework by Vanapalli et al. (1996). More comprehensive explanations related to variations in the shear strength of unsaturated soils are available in Vanapalli (2009). Figure 4 shows the SWCC and the comparison between the measured bearing capacity values (Mohamed & Vanapalli, 2006) and those predicted using the VM model for Unimin 7030 industrial sand. The maximum average matric suction value used in the model footing tests was 6 kPa due to limitations of the tank size (i.e., height) used in the testing program; therefore, no comparisons are available beyond this suction value. However, the predicted bearing capacity clearly shows the bearing capacity starts decreasing as matric suction value approaches the residual zone of desaturation.
The VM model was developed by extending the approach shown in Figure 5. For a matric suction value equals to or less then air-entry value, the bearing capacity increases linearly with increasing matric suction. Hence, the bearing capacity for air-entry value can be calculated using Equation 17 assuming the soil is in a state of saturated condition.
where (u a -u w ) b is air-entry value, B is width of footing, γ is unit weight of soil, N c , N γ are bearing capacity factors, and ξ c , ξ γ are shape factors The degree of saturation reduces significantly and reaches a value close to zero at high suction values in the residual zone. In this zone, the contribution of matric suction towards the bearing capacity significantly reduces, which is similar to the shear strength behavior of sands. The contribution from suction which is represented by second term in Equation 18 becomes negligible. In other words, for a suction value greater than air-entry value, the predicted bearing capacity obtained VM model is equal to or greater than the one calculated with Equation 17 regardless of matric suction value. To overcome this limitation of the VM model, in this paper, the VM model is improved by adding  For a matric suction greater than air-entry value, the bearing capacity is represented as the sum of BC sat , BC AVE and BC unsat as shown in Equation 18.
The model footing test results presented by Safarzadeh & Aminfar (2020) were revisited to check the validity of the revised VM model. The tests were carried out with a square footing (B × L = 150 mm × 150 mm) in a poorly graded Goomtapeh sand. The effective cohesion, effective internal friction angle and dilatancy angle are 0 kPa, 33.6°, and 6°, respectively. Their experimental results were chosen in the study since one test was carried out at a relatively high suction value that is close to the dry condition (i.e., 30 kPa). Figure 6 shows the SWCC and the variation of bearing capacity with respect to matric suction. As shown in Figure 6b, the measured bearing capacity at 30 kPa of matric suction value (i.e., 73 kPa) is slightly higher than that of saturated condition (i.e., 48 kPa). This indicates that the bearing capacity decreases significantly as the degree of saturation approaches zero. This behavior was not captured in the original VM model due to the reason explained earlier. The revised VM model provides good agreement between the measured and predicted bearing capacity values over entire range of SWCC (i.e., with three different zones; namely, boundary, transition and residual zone). The discrepancy in the predicted bearing capacity values between the original and revised model can be noticed at the matric suction of approximately 10 kPa. This is the point where the bearing capacity starts decreasing rapidly in the residual zone of desaturation. Terzaghi (1943) bearing capacity equation was originally developed assuming general shear failure for drained conditions with respect to pore-water. There are also scenarios where shallow foundations are loaded at a relatively fast rate on/in a low permeable soils (i.e., undrained condition). The bearing capacity should be calculated extending the ϕ u = 0 approach (i.e., total stress approach; Skempton, 1948) for such soils. Application of two different approaches depending on the soil type and drainage condition in saturated soils also leads to the argument to the use of MESA for calculating the bearing capacity of unsaturated fine-grained soils. Failure mode from the in-situ plate load tests on unsaturated fine-grained soils represents punching failure, rather than the general failure (Oloo 1994;Schnaid et al., 1995;Consoli et al., 1988;Costa et al., 2003;Rojas et al., 2007). For punching failure condition, i) well defined failure is not observed from the stress versus settlement behavior of model footing tests, ii) no heave is observed on the soils outside the loaded areas, and iii) the vertical displacement of a footing is caused mainly by the compression of the soil directly below the footing as well as the vertical shearing of the soil around the footing perimeter (Vesić, 1963(Vesić, , 1973Oh & Vanapalli, 2013). Oh & Vanapalli (2013) also pointed out that the drainage conditions of porewater and pore-air cannot be clearly defined during in-situ plate load tests on unsaturated fine-grained soils. There are no guidelines in the literature for the duration of plate loading for achieving equilibrium condition. For this reason, Oh & Vanapalli (2013) suggested that it is reasonable to assume that the pore-air is in drained condition, while pore-water is in undrained condition for plate load tests performed in unsaturated fine-grained soils. In other words, the bearing capacity of unsaturated fine-grained soils is more dependent on the compressibility of the soil below a footing (i.e., punching failure) and the constant water content (CW) test can be regarded as the most reasonable test method to simulate the loading and the drainage condition for the unsaturated fine-grained soils. However, since the CW tests are time-consuming and require elaborate equipment (Rahardjo et al., 2004), the unconfined compression test results for unsaturated soils can be used in the MTSA to accommodate both punching failure mode and drainage conditions of pore-air and pore-water in unsaturated fine-grained soils. Oh & Vanapalli (2013) extended this background for their experimental investigations studies and interpreted the bearing capacity of unsaturated fine-grained soils using Equation 20 extending the Skempton (1948) total stress approach. This approach has been referred to as the Modified Total Stress Approach (MTSA). Oh & Vanapalli (2013) suggested using N c(unsat) = 5.14 for achieving good comparisons between the measured and predicted bearing capacity values in unsaturated fine-grained soils. It is interesting to note that the factor, 5.14 is the same as N c for saturated condition with a zero internal friction angle (i.e., undrained loading).

Modified total stress approach
where, s is shear strength based on unconfined compressive strength, ξ [= 1+0.2(B/L)] is shape factor (B = Breadth, L = length) proposed by Meyerhof (1963) and Vesić (1973) for undrained condition, N c is bearing capacity factor, q u is unconfined compressive strength and subscript unsat is unsaturated condition Meyerhof (1974) studied the ultimate bearing capacity of a clay deposit overlain by a homogeneous thick bed of a sand (Figure 7a). For the scenario where the H/B ratio is relatively small, the sand deposit is compressed in a shape of truncated pyramid (i.e., punching failure), which typically follows after a general failure in the clay deposit. If the bearing capacity of sand bed is much greater than that of the underlying clay deposit, the net ultimate bearing capacity of a shallow foundation can be estimated using Equation 21.
where, q b is bearing capacity of the bottom soil layer, q t is bearing capacity of the top soil layer, C a is adhesive force (= c' a H), c' a is adhesive cohesion, H is distance from the base of shallow foundation to the bottom of lower layer, B is width of the shallow foundation, L is length of the shallow foundation, (1 + B/L) is shape factor, D f is embedded depth, K s is punching shear coefficient, ϕ' 1 is effective internal friction angle of the top layer, γ 1 is unit weight of the top layer, and Q u is ultimate load. Oh & Vanapalli (2013) suggested that Equation 21 can be further simplified as Equation 22, assuming that i) bearing pressure at the top of lower layer is q b + γ 1 D f and ii) bearing capacity of the upper layer is governed only by the compressibility of the soil block A-A'-B-B' without lateral deformation (i.e., K s in Equation 21 is zero) (Figure 7b).
The in-situ plate load tests result in an unsaturated cohesive soil by Larson (1997) showed that majority of settlements are within a depth of 2B below the plate, with most of the settlements (i.e., about 87%) within the depth of B. These results justify the use of 1.5B as a significant depth (Poulos & Davis, 1974). Hence, Equation 22 can be further simplified as Equation 23 with the total cohesion, c a for a square footing on a single unsaturated fine-grained soil assuming the bearing capacity is governed by its compressibility (i.e., q b = γ 1 D f = 0). The total cohesion is equal to half of unconfined compressive strength, q u of an unsaturated finegrained soil as shown in Equation 20. Figure 7. (a) Failure of soil below footing on dense layer above weak layer (modified after Meyerhof (1974) and Das (2015)) and (b) simplified failure mechanism.
It is interesting to note that the constant '6' in Equation 23 is close to the value calculated from Equation 20 (i.e., 5.14 × 1.2 = 6.17) for a square footing. Consoli et al. (1988) conducted a series of in situ plate load tests in an unsaturated cohesive soil (lightly cemented homogeneous sandy-silt red clay, I p = 20%) using three different sizes of steel circular plates (0.3, 0.45, and 0.6 m in dimeter) and concrete square footings (0.4, 0.7 and 1 m). Oh & Vanapalli (2013) revisited the plate load test results and estimated bearing capacity values based on the stress versus settlement behaviours. For each of these tests, the bearing capacity was defined as a stress corresponding to the intersection of the tangents to the initial and final portion of stress versus settlement behavior within the settlement equals to 10% of width (or diameter) of a footing (or plate). The bearing capacity predicted using the MTSA (i.e., Equation 20) was 155 kPa, which was significantly lower compared to those from the plate load tests. The discrepancy gradually decreased with an increase in the plate size and becomes negligible when the plate size is of 1 m (see Figure 8). The discrepancy between the measured and predicted bearing capacity can be attributed various factors such as the disturbance of test specimen used for unconfined compression test and the scale effect (i.e., representative or average matric suction value can be different depending on the footing size). The summarized discussion validates the use of MTSA (i.e., Equation 20) for the prediction of the bearing capacity of fine-grained soils for full-size foundations used in engineering practice.
The MTSA (i.e., Equation 20) indicates that the variation of bearing capacity of unsaturated fine-grained soils with respect to matric suction can be obtained by estimating the undrained shear strength of a soil at a targeted matric suction value. Oh & Vanapalli (2018a) proposed a model to estimate the variation of undrained shear strength of unsaturated fine-grained soils with respect to matric suction using the unconfined compressive shear strength of unsaturated soil.
where c u (unsat) and c u(sat) = undrained shear strength of unsaturated and saturated soil (i.e., half of unconfined compressive strength), respectively, P a is atmospheric pressure, S is degree of saturation, and μ, ν are fitting parameters. Oh & Vanapalli (2018a) analyzed seven sets of unconfined compression test results for fine-grained soils. Reasonable estimates were obtained between the measured and predicted undrained shear strengths values for soils within the range of plasticity indices between 8% and 60% using the fitting parameter, ν = 2 (ν = 1 for unsaturated cohesionless soil). The fitting parameter, μ is a function of plasticity index as shown in Figure 9.
The pore-air and pore-water pressure along the shear zone is different from those at the test specimen's boundaries. However, due to the limitations of laboratory testing techniques, the shear strength and bearing capacity models use the pore-air and pore-water pressure applied at the specimen's boundaries for initial or failure condition. The research studies by Chae et al. (2010) and Zhou et al. (2016) suggest that the unconfined compressive strength of unsaturated soils can be predicted by using the soil suction at failure condition. Oh & Vanapalli (2018a) also performed investigations along similar lines by comparing the measured undrained shear Figure 8. Comparison between the measured bearing capacity values and those predicted using the MTSA (data from Consoli et al., 1988).

Figure 9.
Relationship between plasticity index, I p and the fitting parameter, μ in Equation 24 (Oh & Vanapalli, 2018a;Babu et al., 2005;Chen, 1984;Pineda & Colmenares, 2005;Ridley & Burland, 1993;Vanapalli et al., , 2000. strengths with those predicted using Equation 24 using the experimental results of Cunningham et al. (2003). Figure 10 summarizes these results highlighting the degree of saturation and soil suction for initial (subscript i) and failure (subscript f) condition. The undrained shear strength predicted using both [S i , (u a -u w ) i ] and [S f , (u a -u w ) f ] showed good agreement with the measured values when plotted against (u a -u w ) f . This indicates that Equation 24 can be used to estimate the variation of undrained shear strength with respect to matric suction without the information of degree of saturation and soil suction at failure. Such trends in results can be attributed to the reason that the increase in soil suction at failure leads a decrease in the degree of saturation and vice versa, which minimizes the difference in undrained shear strength obtained

Numerical approaches
Various numerical approaches can be used for predicting the bearing capacity by simulating the vertically applied stress versus settlement behaviors. Abed & Vermeer (2006) were some of the earliest investigators who estimated bearing capacity of soils taking account of saturated and unsaturated soil conditions using the Barcelona Basic Model (BBM) and the Modified Cam Clay (MCC) model, respectively. Ghorbani et al. (2016) used an extended MCC to simulate load versus settlement behaviours of statically loaded rigid footings in both saturated and unsaturated soils. Elastic-perfectly plastic Mohr-Coulomb model was used by several researchers to estimate the behaviors of shallow foundations in unsaturated soils (Tang et al., 2016;Baker, 2004;Serrano et al., 2005;Oh & Vanapalli, 2011;Cheng et al., 2021). In this section, details of numerical methods proposed by Oh & Vanapalli (2011, MESA) and Oh & Vanapalli (2013, MTSA) to estimate the bearing capacity of unsaturated sand and fine-grained soils, respectively, are summarized. Mohamed & Vanapalli (2006) performed a series of model footing (B × L = 100 mm × 100 mm) tests on sand. The load was applied on the model footing placed on the surface of the sand housed in a specially designed tank (900 mm × 900 mm × 900 mm). The water table in the tank was adjusted such that the tests were carried out to achieve four different average matric suction values (i.e., 0, 2, 4, and 6 kPa). Figure 11 shows the schematic of testing set up and the matric suction distribution with depth with the water table at the depth of 600 mm from the soil surface. Matric suction values were measured using Tensiometers embedded at four different depths. The difference between the measured and assumed hydrostatic matric suction distribution was not significant. Hence, the four average matric suction values were obtained with assumed hydrostatic matric suction distributions. Oh & Vanapalli (2011) conducted numerical analyses to estimate the bearing capacity extending the MESA using the model footing test results by Mohamed & Vanapalli (2006). Figure 12 shows the soils properties and boundary conditions used in the numerical analysis. The total cohesion was estimated using Equation 25 (Vanapalli et al., 1996). The coefficient of earth pressure at rest, K 0 was calculated using Equation 26 (Jaky, 1944) and the Poisson's ratio, μ s was assumed equal to 0.3. The initial tangent elastic modulus, E i for the saturated condition was obtained using Equation 27, which was used to estimate the variation of E i with respect to matric suction (Equation 28 24) undrained shear strengths using the data from Cunningham et al. (2003). Degree of saturation and soil suction for both initial and failure conditions were used in the prediction of the undrained shear strength (modified after Oh & Vanapalli (2018a)).

Estimating stress versus settlement behaviour extending MESA
where S is degree of saturation and κ is fitting parameter that is a function of plasticity index (Garven & Vanapalli, 2006 where K 0 is coefficient of earth pressure at rest, E i(sat) is initial tangent elastic modulus for saturated condition and (Δδ/Δq) is the slope of the settlement versus footing pressure in elastic range.
where E i(unsat) is initial tangent elastic modulus for unsaturated condition, P a is atmospheric pressure, and α› and β› are fitting parameters. Oh et al. (2009) studies suggest that β› = 1 and 2 can be used for coarse-and fine-grained soils, respectively. Vanapalli & Oh (2010) investigated plate load and model footing test results in unsaturated soils and proposed that the fitting parameter, α› is a function of plasticity index, I p as shown in Figure 13. Although α› = 2.5 provided best comparison between the measured and estimated E i(unsat) for sandy soils, Oh et al. (2009) recommended α› values between 1.5 and 2 for conservative elastic settlement estimation in engineering practice applications. Figure 14 shows the applied vertical stress versus surface settlement behaviors obtained from the model footing test and numerical analysis for an average matric suction value of 4 kPa. Figure 15 shows the comparison between the measured bearing capacity values and those predicted using the VM model (Equation 18) and numerical method. The bearing capacity values predicted using the VM model with ϕ' and numerical method with 1.1ϕ' show better agreement when compared to the measured values. This can be attributed to the different failure mechanism between VM model (i.e., limit equilibrium method) and numeral method (i.e., Mohr-Coulomb criteria) (Oh & Vanapalli, 2011).   (2010)).

Figure 14.
Comparison of the applied stress versus surface settlement behaviors obtained from model footing tests (Mohamed & Vanapalli, 2006) and numerical analysis for the average matric suction values of 4 kPa (modified after Oh & Vanapalli (2011)).

Estimating stress versus settlement behaviour extending MTSA
Oh & Vanapalli (2018b) proposed a modeling technique to simulate the vertical stress versus surface settlement behaviors of shallow foundations in unsaturated fine-grained soils extending the MTSA. As discussed earlier, the bearing capacity of saturated fine-grained soil is calculated extending ϕ u = 0 approach. Oh & Vanapalli (2018b) suggested that the same approach can also be extended to model the vertical stress versus surface settlement behaviors of shallow foundations on unsaturated fine-grained soils. This approach is consistent with the concept used in developing Equations 20 and 23. The bearing capacity of unsaturated fine-grained soils can be predicted with Equation 20 using the unconfined compressive strength of unsaturated fine-grained soil.
The variation of elastic modulus with respect to matric suction was estimated by modifying Equation 28 for undrained condition. In other words, the elastic modulus for undrained condition was first estimated using the slope of the stress versus strain behaviour from unconfined compression test results and then Equation 28 was used to estimate the variation of elastic modulus with respect to suction for undrained condition. Figure 16 shows the mesh and boundary condition used in the numerical analysis to simulate the applied stress versus surface settlement behaviours of model footing tests from Oh & Vanapalli (2013). Figure 17 shows the comparison between the measured applied stress versus surface settlement behaviour and those predicted with numerical method. Good comparisons were obtained when the applied stress versus surface settlement behaviors were estimated with Poisson's ratio, μ s = 0.45 and 1 for saturated and unsaturated conditions, respectively. The bender element test results by Lee & Santamarina (2005) also supports μ s = 0.1 for unsaturated soils. The numerical analysis results also showed that the applied stress versus surface settlement behaviors is not affected by the coefficient of earth pressure at rest, K 0 . Figure 18 shows comparison between the measured bearing capacity values and those estimated from MESA (revised VM model), MTSA (Equation 20), and numerical method extending the MTSA. The bearing capacity values predicted with numerical analysis extending the MTSA provided the best agreement when compared with the measured bearing capacity values.  Figure 17. Comparison between the measured stress versus settlement behaviours (Oh & Vanapalli, 2013) and those predicted with numerical method for various Poisson's ratio (modified after Oh & Vanapalli (2018b)).

Summary and conclusions
Many attempts have been made to improve the Terzaghi (1943) bearing capacity model extending effective stress approach (ESA) since its development. These models are valuable for the practicing engineers to estimate the bearing capacity of saturated soils considering soil type (i.e., coarseand fine-grained soils), drainage condition (i.e., drained and undrained condition), shape and embedment depth of footings. However, there are many scenarios where the influence of soil suction has to be taken into account for better estimation of the bearing capacity of unsaturated soils. More recently, researchers developed analytical models for this reason considering the influence of matric suction. Some of these analytical models have been successfully used in numerical analysis to simulate the applied stress versus settlement behaviours of shallow foundations in unsaturated soils. Unlike saturated soils, bearing capacity of unsaturated soils have been predicted extending the Modified Effective Stress Approach (MESA) regardless of soil type and drainage condition of pore-air and pore-water. Oh & Vanapalli (2013) provided a different approach and validated the use of Modified Total Stress Approach (MTSA) for unsaturated fine-grained soils. In this state-of-the-art paper, authors revisited the bearing capacity models for unsaturated soils and discuss the use of models by categorizing them into two groups; MESA and MTSA.
The model footing test results on unsaturated cohesionless soils show that the bearing capacity at high suction values approaches the bearing capacity values for saturated condition. This is because of the contribution of matric suction towards the bearing capacity is negligible at high suction values (i.e., in the residual zone of desaturation). To accommodate this behavior, the bearing capacity model proposed by Vanapalli & Mohamed (2007, VM model) extending the MESA is revised introducing additional criteria. The revised VM model reasonably captures the decrease in the bearing capacity in the residual zone of SWCC, which was verified using model footing test results for unsaturated sand. Oh & Vanapalli (2013) showed that the bearing capacity model for saturated soils under undrained condition originally proposed by Skempton (1948) can be used for unsaturated soils by simply replacing the undrained shear strength for saturated condition with a value of half the unconfined compressive strength for unsaturated soil. The MTSA is simpler and provides better estimates when compared to those estimated with MESA based bearing capacity models. This may be attributed to the reasons that the MTSA was developed taking account of typical failure mode (i.e., punching failure) and the drainage condition that is more appropriate for unsaturated fine-grained soils beneath the shallow foundations.
Several constitutive models and modeling techniques are available in the literature that can be used to estimate deformation in unsaturated soils considering both drained and undrained conditions. However, these constitutive models require various soil parameters that require extensive laboratory test results. Typically, laboratory tests for unsaturated soils are time consuming and cumbersome since the pore-air and porewater should be controlled separately. On the other hand, the numerical modelling techniques proposed by Oh & Vanapalli (2011) and Oh & Vanapalli (2018b) extending the MESA and MTSA, respectively are simple and requires only the laboratory test results for saturated conditions. These numerical techniques are validated by comparing measured bearing capacity values with those obtained with numerical methods.
This study focuses on the prediction of ultimate bearing capacity (i.e., ultimate limit state) considering the influence of matric suction. The bearing capacity models for unsaturated soils use soil suction as a state variable, which suggests, the estimated bearing capacity values are governed by soil suction. Hence, it is important to use appropriate climate models to predict the variation of bearing capacity of unsaturated soils with respect to soil suction in geotechnical engineering practice. However, in several scenarios, settlement behaviour (i.e., serviceability limit state) is the governing parameter in the design of shallow foundations. Therefore, both bearing capacity and settlement should be taken into account in the rational design of shallow foundations since foundation settlement can exceed the allowable settlement prior to the bearing capacity reaching its ultimate value. According to Tidlund (2021), uncertainty in geotechnical engineering is much more common in comparison to other civil engineering fields. Uncertainty in the estimation of foundation settlement in many scenarios is more likely compared to that of bearing capacity. The estimated foundation settlement using an idealized settlement calculation model without sufficient data can be 'nothing more than an educated guess based on Figure 18. Comparison between the measured bearing capacity values (Oh & Vanapalli, 2013) and those estimated using different approaches (modified after Oh & Vanapalli (2018b)). experience and the judgements of the engineer (Tidlund, 2021). For this reason, Szavits-Nossan (2006) suggested that the observational method is a more reliable tool in the design of shallow foundations, when it is governed by serviceability.
Sufficient site investigation and laboratory tests are required to determine mechanical properties (e.g., effective shear strength parameters, unconfined compressive strength and hydraulic conductivity for saturated condition), active zone based on seasonal base, and the SWCC for reliable use of the MESA or MTSA. More importantly, reliable climate data of the construction site is required for implementing the MESA or MTSA.
The bearing capacity of unsaturated soils can be predicted using either the MESA or MTSA depending on the soil type and drainage conditions of pore-air and pore-water. As per discussions in this paper, the MESA is recommended for unsaturated coarse-grained soils since both pore-air and pore-water are under drained conditions during loading stages. However, in many scenarios of engineering practice for unsaturated fine-grained soils, there are uncertainties with the drainage conditions of pore-air and pore-water during loading stages due to the low coefficient of permeability of the soil. Except for the scenarios, where load is applied at a significantly slow rate, it is reasonable to assume that the drained conditions of both pore-air and pore-water are somewhere between drained and undrained condition. Hence, authors suggest that the bearing capacity of unsaturated finegrained soils should be estimated using both the MESA and MTSA; however, the final decision be made by engineers based on their engineering judgement.