Undrained capacity of circular shallow foundations on two-layer clays under combined VHMT loading

Wind turbines are typically designed based on fatigue and serviceability limit states, but still require an accurate assessment of bearing capacity. Overconsolidated clay deposits in Canada often have a thin layer of crust with a relatively high undrained shear strength. However, existing bearing capacity design methods do not consider surficial crusts. This paper studies the undrained VHMT (vertical, horizontal, moment, and torsional) failure envelope of circular foundations founded on a surficial crust underlain by a uniform soil using finite element analysis. Crust correction factors have been introduced to account for the effects of the stiff layer on the vertical and moment capacities. The same forms of equation that are used for uniform soils, but with different parameters provide satisfactory fits for the failure envelopes for a soil with a crustal layer. An analytical expression for the 4-D VHMT failure envelope is derived, and an application of this method is provided.


Introduction
Estimation of the bearing capacity of shallow foundations under combined loadings can be of great significance. This is particularly important for large onshore structures, such as wind turbines, transmission towers and masts, due to their complex environmental loadings. The majority of offshore foundations are circular or close to circular in form. For wind turbines, apart from the vertical load due to its self-weight, the horizontal load induced by wind is also substantial and a large structural height further leads to a significant moment load on the foundation.
Traditionally, this type of design is based on classical solutions for the uniaxial vertical bearing capacity of shallow foundations using the superposition principle (Terzaghi, 1951). The effects of load inclination and eccentricity are incorporated by introducing a load inclination factor and an effective foundation area (e.g. DNV, GL, 2016). Since load inclination and eccentricity effects are separately considered, this semi-empirical modification of the conventional theory may sometimes be insufficiently accurate for practical design (Gourvenec, 2007). This approach is still common for onshore shallow wind turbine foundation design.
The failure envelope method has been recommended as an alternative to conventional theories in some geotechnical design guidelines (particularly those focused on offshore geotechnics), such as API (2011) and International Standards Organisation (ISO) (2016), due to the load interaction effect between various load components (i.e. combined vertical, horizontal, moment, and torsional loads) being explicitly incorporated (Shen et al., 2017). Failure envelopes for different types of foundations (e.g. strip (Bransby and Randolph, 1998;Fathipour et al., 2022), rectangular (Gourvenec and Randolph, 2003) and circular (Suryasentana et al., 2021) foundations), homogeneous (Taiebat and Carter, 2010) or non-homogeneous (He and Newson, 2022b;Wu et al., 2022) soils, and zero-tension (Shen et al., 2016) or unlimited-tension (Gourvenec and Randolph, 2003) interface conditions have been previously investigated under undrained soil conditions. These studies have been primarily confined to a single layer soil with a uniform or linearly increasing soil undrained shear strength profile. However, an onshore clay deposit often has a thin layer of stiff crust with a relatively high undrained shear strength developed from weathering, desiccation, and chemical process (Lutenegger, 1995). The shear strength of the upper crustal layer can be more than 10 times that of the underlying clay (Lee and Park, 1999). Nakase et al. (1978) and Sagaseta and Arroyo (1982) have demonstrated that the undrained shear strength profile of the crust strongly affects the stability analysis of shallow foundations and embankments because a substantial portion of the failure surface under the structure can be located within the crust. Understanding the effect of a surficial crust on the bearing capacity of shallow foundations is therefore important for optimal foundation design. However, existing design standards (e.g. DNV, GL, 2016) predict the bearing capacity of shallow foundations under idealized soil conditions without considering the existence of the surficial crust.
Studies available in the literature that account for the effect of a high strength surficial layer on foundation bearing capacity are still sparse. Merifield et al. (1999) evaluated the undrained bearing capacity of a centrally, vertically loaded surface strip foundation on a two-layer clay deposit using numerical upper and lower bound analysis. Recently, the same cases for square and circular foundations were further investigated by Merifield and Nguyen (2006) using finite element (FE) analysis. However, this considers only vertical bearing capacity in the absence of horizontal, moment and torsional loads. Park et al. (2010) determined the bearing capacity factor N c for strip and circular foundations resting on a non-homogeneous crust overlying a uniform soil, but the failure envelopes under combined loadings were not addressed. Feng et al. (2015) studied the failure envelope of a rectangular foundation founded on a crustal layer overlying a normally consolidated clay using an unlimited-tension soil-foundation interface. However, unlike skirted shallow foundations of offshore structures, onshore shallow foundations can uplift and separate from the soil under a large overturning moment because the soil-foundation interface cannot provide tensile resistance (i.e. zero-tension interface). Xia et al. (2022) explored the VHM envelope of bucket foundations in a two-layer clay, but the soil is soft-over-stiff rather than stiff-over-soft.
The object of this paper is to investigate the VHMT failure envelope of circular foundations founded on a surficial crust underlain by a uniform soil under a zero-tension interface condition using finite element analysis. Undrained soil conditions have been considered for both the crustal layer and the underlying soil. The effects of the surficial crust on the failure envelope have been studied. An analytical expression for the 4-D VHMT failure envelope is also derived in this paper. The derived 4-D failure envelope has been briefly summarized in a technical note by the authors (i.e. He and Newson, 2022a), but the detailed procedures to derive the envelope were not provided. As a complement, this paper presents the analysis of deriving the final 4-D envelope in more detail and, an example of the use of this method for design is also provided.

Material models and interface conditions
A linear elastic perfectly plastic constitutive relationship with a Mohr-Coulomb (M-C) failure criterion was used to model the soil behavior. For undrained soil conditions, the M-C criterion degenerates to the Tresca criterion, which is defined by three soil parameters: the undrained Young's modulus, E u , Poisson's ratio, m, and undrained shear strength, s u .
As shown in Figure 1a, two fundamental parameters of the surficial crust layer may affect the bearing capacity of a surface foundation on a crusted soil: the average crust undrained shear strength, s ut , and the crust thickness, t c . As suggested by Lee and Park (1999), in some cases the average shear strength of the upper crust layer can be more than 10 times that of the underlying clay. In this paper, the ratio, s u0 /s ut (s u0 refers to the undrained shear strength of the underlying soil), has been varied from 0.2 to 1.0 (i.e. 0.2, 0.4, 0.6, 0.8, and 1.0), to represent a typical range of soil profiles. This approach addresses cases of strong-over-soft clays adopted by Merifield et al. (1999). The special case of s u0 /s ut = 1 relates a homogeneous soil. As shown by Bjerrum (1973), the thickness of a crust layer usually ranges from 1 to 8 m depending upon the hydrogeology (i.e. well-drained versus poorly drained). In eastern Canada, the crust thickness is generally 1-5 m and is often of the order of 3 m (Lefebvre et al., 1987). Typically, the diameter (D) of an onshore wind turbine foundation is very large (.15 m), and these foundations are becoming larger with the increase of power output and tower height. The foundation diameter used in this paper is 19 m, representing the typical dimension for current onshore wind turbines in Canada. Therefore, this study has considered models with t c /D ranging from 0.1 to 0.3 (i.e. 0.1, 0.2, and 0.3) to span most cases of practical interest. This range also covers that used by Feng et al. (2015).
As an initial study, the undrained shear strength of the underlying soil (i.e. s u0 ) is considered to be constant with depth for simplicity. In the analysis, s u0 was held constant at 100 kPa which is a typical value for the underlying soil in Ontario Canada. The undrained soil Poisson's was taken to be 0.495. Since the soil Young's modulus, E u , affects only the evolution of the load-displacement relationship but not the calculated collapse loads (Yun and Bransby, 2007), a sufficiently large E u /s u0 ratio equal to 10,000 was selected to minimize mesh distortion (Abyaneh et al., 2015). The foundation was assumed to be a rigid body, and a load reference point (LRP) was used to apply prescribed displacements or loads, located at the bottom center of the foundation.
Similar to Shen et al. (2016), the FE model considered a zero-tension, rough base that allows separation of the foundation from the soil. The zero-tension rough base can be modeled using a Coulomb model with a friction coefficient of 20 (Shen et al., 2016).

Geometry and mesh
The analysis in this paper was conducted using the finite element software ABAQUS (Dassault Syste`mes, 2016). To avoid the effects of model boundaries on the development of the failure mechanisms, the mesh length, L, and mesh height, H, were taken as 120 and 50 m, respectively, following the recommendation of Deshpande (2016).
A mesh convergence study was carried out for a number of cases. A typical outcome is shown in Figure 2. The model solution with Mesh 3 takes about six times longer than that using Mesh 2, but the difference between the ultimate vertical capacities using Meshes 2 and 3 is only around 2%. Thus, this study adopted Mesh 2 in the following analyses. Figure 3 shows the half-model of the 3-D model using Mesh 2 (He and Newson, 2022a). This mesh was composed of approximately 36,000 8-noded brick elements (i.e. first-order, ABAQUS C3D8R with reduced integration and hourglass control). Following Gourvenec and Randolph (2003), the soil regions in the vicinity of the foundation edge and the horizontal thin soil layer close to the interface were carefully refined to capture the intense stress concentration close to the foundation edge and the large plastic shear strains at the interface. The cylindrical circumference of the soil mass was constrained to prevent out-of-plane translations, and the bottom of the soil mass was fixed in the three orthogonal directions.

Sign conventions and loading paths
The sign conventions for the loads are shown in Figure 1b. The horizontal and moment loads were considered to be in the same plane.
Probe tests and swipe tests were used to track the failure envelopes under different loading conditions. One probe test provides only a single point on a failure envelope. A swipe test is able to track the whole or a portion of a failure envelope. Refer to He and Newson (2022b) for more details of probe tests and swipe tests. Two typical failure envelopes obtained using the swipe and probe tests are shown in Figure 4.

Finite element results
As shown in Figure 5, the overall strategy for assessing the undrained bearing capacity of a shallow foundation under combined VHMT loads is decomposed into three steps: (i) determining uniaxial ultimate capacities, for example V ult , H ult , M ult , and T ult ; (ii) normalizing H-V, M-V, and T-V failure envelopes with the corresponding uniaxial ultimate capacities; and (iii) normalizing M-H, H-T, and M-T failure envelopes with the corresponding maximum values. Using the results from the above steps, the failure envelope under combined VHMT loadings can be derived.

Pure uniaxial capacities
The ultimate loads for vertical, horizontal and torsional modes are referred to as the corresponding uniaxial loadcarrying capacities in the absence of the other loading modes. Note that the torsional and moment loads on an onshore wind turbine foundation result from the effects of horizontal loading (i.e. horizontal load non-coplanarity and tower height). As a foundation with a zero-tension interface cannot resist any moment loading without vertical loads, the ultimate moment capacity is represented by the maximum moment load under vertical loading only.
Analyses have shown that the uniaxial horizontal and torsional capacities of a surface foundation are purely related to the undrained shear strength of the surface soil: H ult ¼ As ut and T ult ¼ ADs ut =3, which has also been confirmed by previous studies (e.g. Finnie and Morgan, 2004;He and Newson, 2020). Therefore, only the vertical and moment capacities have been investigated in this section. For validation purposes, the current results for a single-layer clay are compared to those from previous studies (Shen et al., 2016;Taiebat and Carter, 2010) in Table 1. It can be seen that the comparison is satisfactory.
The crust correction factor, s cr , defined as the ratio of the dimensionless capacity for a crusted soil to that of a uniform soil, is introduced to characterize the effect of a surficial crust: where N v s u0 =s ut ¼ V ult = As ut ð Þ for the vertical capacity and N m s u0 =s ut ¼ M ult = ADs ut ð Þfor the moment capacity. Note that the crust correction factors for horizontal and torsional loading modes remain at unity. The variations of s crV and s crM with regard to s u0 /s ut and t c /D are shown in Figure 6 (He and Newson, 2022a). The vertical and moment factors significantly increase with s u0 /s ut and gradually converge to unity as s u0 /s ut approaches unity. A quadratic polynomial equation with respect to s u0 /s ut is proposed to estimate the relationships: The coefficient, f, is a function of t c /D, defined as f t c =D ð Þ¼À0:97t c =D À 0:27 for the vertical capacity and f t c =D ð Þ¼ À1:18t c =D t c =D + 0:18 for the moment capacity. As seen in Figure 6, the curve fitting is in close agreement with the FE results.
To better understand the effects of s u0 /s ut and t c /D on the collapse mechanism at V ult and M ult , Figure 7 shows the contours of the maximum plastic shear strain increment for a uniform soil (Figure 7a and d), the case of s u0 /s ut = 0.2 and t c /D = 0.2 (Figure 7b and e), and the case of s u0 /s ut = 0.6 and t c /D = 0.3 ( Figure  7c and f). The geometry of the collapse mechanism can be effectively visualized using the maximum plastic shear strain increment (Loukidis et al., 2008): where De p xx , De p yy , De p zz , Dg p xy , Dg p xz and Dg p yz are the plastic normal and shear strain increments in Cartesian coordinates. Figure 7a-c shows the failure mechanisms at V ult . As shown in Figure 7a, the failure mechanism for a uniform soil compares well to that of general shear failure (white dashed lines) found from the method of characteristics (Martin, 2003). The comparison between Figure 7a-c shows that the depth of the active triangular zone remains almost the same, although the active triangular zone for the case of s u0 /s ut = 0.2 and t c /D = 0.2 is slightly curved. However, the shear fan zone for the case of s u0 /s ut = 0.2 and t c /D = 0.2 extends to a depth of 1D, but ends at the  base of the crust. In addition, the case of s u0 /s ut = 0.2 and t c /D = 0.2 has no passive zone close to the ground surface, indicating a local shear failure mode. This is because the relatively strong top crust acts as rigid column that restricts both upward and lateral deformations within the crustal layer, while this restriction in turn increases the depth of the failure zone within the bottom layer. A partial shear failure mechanism also appears to be initiating in the lower layer for the two crusted cases. This phenomenon was also observed for both square and circular foundations by Merifield and Nguyen (2006). Figure 7c shows that the failure mechanism at V ult for the case of s u0 /s ut = 0.6 and t c /D = 0.3 lies in between Figure 7a and b, since this case has a moderately strong surficial crust. The collapse mechanisms at M ult (at V/V ult = 0.50) are shown in Figure 7d-f. For the uniform soil, a combined scoop-wedge mechanism is observed (similar to that found by Bransby and Randolph, 1998). Compared to the uniform soil case, the failure zone for the case of s u0 /s ut = 0.2 and t c /D = 0.2 is primarily confined to the underlying layer, and the crust behaves as a rigid column due to the relatively high strength. A wedge mechanism still exists on the right hand side of the foundation, but the scoop seems to have been suppressed. Moreover, the failure pattern for the case of s u0 /s ut = 0.6 and t c /D = 0.3 lies in between the other two cases. Therefore, the main effect of a crustal layer on the collapse zone is the suppression of the surface failure and the scoop portion of the mechanisms, and an increase of the depth of the failure zone within the underlying layer.   Figure 8a exhibit an expansion of the curves with increasing s u0 /s ut . Note also that a stiffer H-V failure envelope is observed for a smaller value of s u0 /s ut , that is a lower s u0 /s ut ratio corresponds to a larger normalized H-V envelope, as shown in Figure 8b. For rectangular foundations on a soil with a crust using an unlimited-tension interface, Feng et al. (2015) proposed an analytical equation (see equation (4)) as a function of the loading angle, u. The curves with u = 0 and 90°are presented in Figure 8b. In addition, a curve fit to Green's exact solution (equation (5)) is also widely used to describe H-V envelopes. The comparison shows that equation (4) with u = 0 and equation (5) provide conservative predictions of the FE results. For simplicity, equation (5) is considered to model the current H-V failure envelope by ignoring its slight dependence on s u0 /s ut . The cases of s u0 /s ut = 0.6 and t c /D = 0.1-0.3 are presented in Figure 8c and d to show the effect of the normalized crust thickness t c /D on the H-V failure envelopes. As shown in Figure 8c, the envelope size increases with t c /D, but the rate of increase appears to gradually decrease. It can also be seen that equation (5) also provides a relatively conservative evaluation of the normalized curves, although dispersion of the curves caused by t c /D can be observed. Figure 9 shows the effects of s u0 /s ut and t c /D on the M-V failure envelopes. A significant expansion of the failure loci with increasing s u0 /s ut and t c /D is observed in Figure 9a and c, although the rate of expansion gradually decreases with s u0 /s ut . The failure envelopes in terms of loads normalized by their ultimate values in Figure 9b and d fall in a very tight band, with the shape following the parabolic function given by equation (6).

Moment-vertical loading
The equation proposed by Feng et al. (2015) for rectangular foundations on a soil with a crust under an unlimitedtension interface is also incorporated in Figure 9, which differs considerably from the current FE results due to the unlimited-tension interface.  Figure 10 shows the effects of s u0 /s ut and t c /D on the T-V failure envelopes. Similar to the horizontal capacity, the torsional capacity is also determined only by the surface soil strength, therefore, the T-V failure envelopes share similar features with increasing s u0 /s ut and t c /D to the H-V envelopes. The analytical equations proposed by Feng et al. (2014) (see equation (7)), Abyaneh et al. (2015) (see equation (8)) and Feng et al. (2015) (see equation (9)) are compared to the FE results in Figure 10b and 10d. The comparison shows that the formula of Abyaneh et al. (2015) gives reasonable and conservative predictions for the normalized failure envelope.

Moment-horizontal loading
The effects of s u0 /s ut and t c /D on the M-H failure envelopes are illustrated in Figure 11. As the failure envelopes for V/V ult = 0.25, 0.50, and 0.75 are similar in shape (only sizes are different), only the results for V/V ult = 0.50 are presented for brevity. The equation proposed by Feng et al. (2015) for rectangular foundations on a soil with a crust under an unlimited-tension interface condition is also compared. Figure 11a and b show that with the increase of s u0 /s ut and t c /D, an expansion of the failure envelopes is observed, but a similar shape of the curves is expressed. This feature may assist in eliminating the dependence on s u0 /s ut and t c /D levels by normalizing the failure envelopes by their corresponding maximum values, as shown in Figure 11c and d. Moreover, the current M-H curves are almost symmetrical about H = 0, while the curve proposed by Feng et al. (2015) is oblique due to the effect of the unlimited-tension interface. In this study, a unique equation, expressed as: is used to simulate the normalized envelopes. The comparison shows that this simple expression provides reasonable fits for various s u0 /s ut and t c /D levels.

Horizontal-torsional loading
Similar to the M-H envelopes shown in Figures 11 and 12 shows the effects of s u0 /s ut and t c /D on the H-T failure envelopes at V/V ult = 0.50. The dimensionless envelopes shown in Figure 12a and b exhibit a high similarity regardless of the values of s u0 /s ut and t c /D. As shown in Figure 12c and d, the equation proposed by Feng et al. (2015) lies slightly inside the normalized envelopes. In addition, equation (11) proposed by Finnie and Morgan (2004) is considered to fit the H-T failure envelopes. The dimensionless powers, l = 1.5 and n = 1.95, in equation (11) yield satisfactory fits, as can be seen from Figure 12c and d.

Moment-torsional loading
The effects of s u0 /s ut and t c /D on the ultimate load-carrying capacity under moment and torsional loading at V/V ult = 0.50 are shown in Figure 13. As shown in Figure 13a and b, the M-T envelope expands with increasing s u0 /s ut due to the increase of the maximum moment. Normalization by the corresponding maximum moment and torsional loads also leads to a relatively narrow band of the curves. The comparison shows that the expression given by Feng et al. (2015) seems to be slightly unconservative. The analytical relationship of equation (12) with the two dimensionless parameters equal to 2 is compared to the FE results in Figure 13c and d, where favorable predictions are observed. to derive the final VHMT envelope expression. Due to the very similar shape of the M-H failure envelope (only the sizes are different), it is reasonable to assume that under the condition of T 6 ¼ 0, equation (14) is still applicable for the M-H failure envelope normalized by the corresponding maximum values, H 0 max and M 0 max (reduce to H max and M max in equation (14) if T = 0). Therefore, equation (14) can be replaced by a more generalized form: Mathematical manipulations of equations (13), (15), and (16) allow the formulation of an analytical expression for the failure envelope in VHMT loading space in terms of V =V ult , H=H ult , M=M ult , and T =T ult : In practical design, the design loads (factored loads and materials), VHMT, can be directly substituted into the left-hand side of equation (17); values less than 1 represent a sufficient ultimate limit design and vice versa. To visualize the shape of the 4-D normalized failure surface, three 3-D failure surfaces in terms of V =V ult , H=H ult , M=M ult , and T =T ult are presented in Figure 14 (He and Newson, 2022a). The specific curves obtained from the FE results are also incorporated for comparison. For the VHT and VMT failure surfaces, the portion of T\0 is presented considering the symmetry about the plane of T = 0.

Application for design
As an example of the application of the method, the ultimate limit design of a shallow foundation of a typical onshore wind turbine is addressed using the failure envelope method. The diameter of the surface circular foundation is 19 m, and the undrained shear strength estimated with a cone penetration test (CPT) is shown in Figure 15. Due to the existence of fissuring and stress-strain compatibility, the high undrained shear strength of the upper crust is unlikely to be fully mobilized at failure (Lefebvre et al., 1987). Rochelle et al. (1974) recommended to use the full (i.e. maximum strength value), mid-depth (i.e. strength value at mid-depth of the crust), and minimum (i.e. strength value just below the crust) strength of the upper crust for limit state analysis. For the s u profile shown in Figure 15, the mid-depth strength is close to the maximum value. Lefebvre et al. (1987) adopted the mean value of the maximum strength of the upper crust and the minimum strength of the underlying soil. In this analysis, the undrained shear strength of the underlying soil is held constant at s u0 = 80 kPa, and the minimum (s ut = 80 kPa), maximum (s ut = 368 kPa), and mean (s ut = 224 kPa) undrained shear strength of the crust are adopted for comparison. The minimum strength case reduces the soil to a single-layer material. The crust thickness is assumed to be t c = 6.8 m. The assumed undrained shear strength profiles are shown in Figure 15. In addition, the factored ultimate limit state loads from International Electrotechnical Commission (2005)  The corresponding uniaxial bearing capacities are computed based on the factored undrained shear strength (as recommended by DNV, GL, 2016, partial safety factor for materials, g m = 1.25 for total stress analysis), as summarized in Table 2. The safety check for the ultimate limit design using the failure envelope method shows that the three cases yield acceptable ultimate capacity designs. However, the soil with a surficial crust gains extra margins of safety compared to the single-layer one. Figure 16 shows the comparison between the reduced M-H failure envelopes and the factored design loads. These M-H envelopes are obtained by substituting the factored design vertical and torsional loads into equation (17) in order to incorporate the effects of vertical and torsional loads. For comparison, the traditional approach using the load inclination factor and the effective foundation area (DNV, GL, 2016) is also incorporated for the case of a single-layer soil. It can be seen that the traditional approach is more conservative than the failure envelope method, which has also been confirmed by Shen et al. (2016). In addition, the considerable difference in the absolute sizes of the three failure envelopes indicates that ignoring the contribution of the upper crust may significantly underestimate the bearing capacity and lead to a great overdesign of the foundation.