Estimation of Distribution Factor for Peak Penetration Resistance Prediction of Spudcan Foundations in Loose to Medium-Dense Sand Overlying Clay

: A series of 3D ﬁnite element (FE) analyses were performed to estimate the peak penetration resistance of spudcan foundations in sand over clay soil proﬁles. Elasto-perfectly plastic models following Mohr–Coulomb and Tresca failure criteria were used for sand and clay layers, respectively. The coupled Eulerian–Lagrangian (CEL) approach was used to simulate the large deformation in soil that occurs during the spudcan penetration. The performance of the numerical model was validated against centrifuge test results. A parametric study with a broad range of strength parameters for sand and clay was performed. The numerical results were used to assess the inﬂuence of sand thickness ( H s ), the diameter of spudcan ( D ), friction angle of sand, and undrained shear strength of clay ( s u ). A wide range of s u was utilized to predict the resistance both of the soft and sti ﬀ clays. The calculated peak resistances are compared with a published analytical model. It is demonstrated that the model highly overestimates the peak resistance for sti ﬀ clays, most likely because it was developed speciﬁcally for soft clays and, therefore, does not account for the inﬂuence of s u . One of the parameters of the model is revised to account both for s u . Comparisons highlight that the modiﬁed model is able to capture the simulated peak penetration resistance for both soft and sti ﬀ clays. A comprehensive suite of 3D FE analyses incorporating the CEL were performed. An elastic-perfectly plastic model with non-associated ﬂow rule following Mohr–Coulomb failure was used for the sand layer. The same model following


Introduction
Spudcan foundations are widely used to support jack-up barges, which are widely used for offshore construction. For the design and operation of jack-up barges, the estimation of the peak bearing capacity (q peak ) during spudcan installation is essential. A major concern is the installation into sand overlying clayey soil, which may produce large differential settlements and potentially lead to a punch through failure.
Because of the evident similarities between a spudcan and a shallow foundation, bearing capacity equations developed for shallow foundations have been used by SNAME [1]. One of the earliest models to estimate the bearing capacity of shallow foundation in sand over clay profile was proposed by Meyerhof [2] and Hanna and Meyerhof [3]. The procedure, termed as the punching shear method, assumes that a sand block with vertical failure surfaces is pushed into the clay layer. Modifications were proposed to the punching shear method to improve the prediction accuracy. In the projection method, the failure surface is assumed to be inclined at an angle of β [4][5][6]. A range of values have been used for β. Okamura et al. [7,8] performed the centrifuge experiments to measure the bearing capacity of shallow foundations on dense sand overlying clay. Based on the model test measurements, a limit equilibrium-based model that combines the concepts of the punching shear and projection models was presented. Park and Park [9] proposed a new model for bucket foundations on sand over clays profiles based on the outputs of non-linear finite element analyses.
Although there are similarities between two types of foundation, the clear difference is that whereas the shallow foundation can be considered to be wished-into place, the spudcan is pushed into the ground from the surface. Centrifuge model tests have been extensively performed to investigate the peak penetration resistance (q peak ) of spudcan foundations. Teh et al. [10] performed centrifuge model experiments to examine the failure mechanism of spudcan foundations at q peak in sand overlying normally consolidated (NC) clay using the particle image velocimetry technique. It was revealed that at q peak , a sand wedge forms between the base of the spudcan and sand-clay interface. The results of the centrifuge model tests were utilized to develop predictive models for q peak . Lee et al. [11] presented a stress dependent model based on the observations of Teh et al. [10]. The model was calibrated against experimental data for limited dense sand overlying clay profiles. Hu et al. [12] modified the model of Lee et al. [11] to account for the embedment depth. The model was validated against centrifuge model tests.
Alternatively, advanced numerical modelling techniques capable of simulating large-deformation penetration have been used to calculate the penetration resistance of spudcan foundations. Tho et al. [13] utilized the coupled Eulerian-Lagrangian (CEL) approach to simulate the penetration of spudcan for both homogenous and layered soil profiles. The parameters which influence the results of penetration resistance of spudcan were identified. Qiu and Henke [14] performed finite element (FE) analyses using the CEL approach on loose sand overlying soft clay. The influences of sand thickness and bottom clay strength on the bearing capacity of spudcan foundations were evaluated. Qiu and Grabe [15] conducted the CEL analyses using the hypo-plastic model for sand and a visco-plastic model for clay. An equation for q peak was not presented. Hu et al. [16] performed FE analyses using the CEL approach to simulate penetration of a spudcan into medium dense sand over NC clay. q peak results from the numerical simulations were compared with the simplified model proposed by Hu et al. [12]. Hu et al. [17] further extended the CEL approach accounting for the strain-softening behavior of dense sand. The simplified model of Hu et al. [12] was again demonstrated to provide reasonable estimate of q peak .
In addition to deterministic analyses, probabilistic assessments were performed to evaluate the capacity of shallow and spudcan foundations. Chwała and Puła [18] conducted analyses for probabilistic prediction of bearing capacity of shallow foundations on homogenous sand overlying spatially variable clay soil. Shu et al. [19] performed random field FE analyses for probabilistic evaluation of bearing capacity of spudcan foundations embedded at various depths. The undrained strength of the soil was spatially varied with depth. Based on numerical outcomes, it was concluded that the non-stationary random field may generate lower values of mean probabilistic bearing capacity. Li et al. [20] adopted the probabilistic investigation approach to evaluate failure envelops of spudcan foundation under combined vertical, horizontal and moment loading space (VHM). It was concluded that presence of spatially variable soils reduces both the bearing capacity and failure envelopes.
The literature review reveals that the previous studies focus on the penetration characteristics of sand overlying soft clay, because of the potential for a punch-through failure. However, there is a need to estimate the penetration resistance also in sand over stiff clay for operation of jack-up barges supported by spudcan foundations. Hu and Cassidy [21] performed CEL analyses to determine the penetration resistance in sand overlying stiff clay. However, the study focused on the resistance in the clay layer and an equation for q peak was not presented. Based on a limited comparison with available centrifuge test results for s u = 32.6 and 43.9 kPa, it was reported that the Hu et al. [12] model provides favorable predictions of q peak . It should be noted that s u values even for this validation do not cover the range encountered in the field.
In this study, q peak of spudcan foundations in sand overlying a soft to stiff clay layer was calculated using three-dimensional (3D) finite element analyses incorporating the CEL approach. A wide range of strength parameters were used, which include the effective friction angle of the sand layer (φ ) and undrained shear strength of the clay layer (s u ). The calculated bearing capacity results are compared with the analytical model. A new empirical correlation for one of the parameters of the conceptual model is proposed to improve the fit for not only soft clays, but also stiff clays. The enhancement in the prediction accuracy is thoroughly investigated through comparisons with the baseline simulations.

Peak Penetration Resistance Prediction Models in Sand Overlying Clay
Among various penetration resistance models developed specifically for spudcans, two q peak models widely used in practice are introduced in this section. Lee et al. [11] proposed an equation to predict q peak of spudcan foundation in dense sand over NC clay using a simplified conceptual model, illustrated in Figure 1a. An inverted truncated cone of sand is assumed to be pushed into the underlying clay. q peak is mobilized by the shearing resistance along the sliding surface and the bearing capacity of the clay layer. It was assumed q peak is mobilized at the surface, based on the centrifuge tests of Lee et al. [22]. The predictive equation is as follows: where N c0 is the bearing capacity factor for clay at the base of a circular foundation, s u is the undrained strength of clay at the sand-clay interface, q 0 is the overburden pressure, H s is the height of the sand layer, D is the diameter of spudcan, ψ is the dilation angle of sand, and E* is an empirical parameter, and γ s is the unit weight of sand. The following equation proposed by Houlsby and Martin [23] for isotropic clays was used for N c0 : where k is the strength gradient of NC clay. E* is defined as follows: where tanφ* is calculated as: D F is defined as: where σ n is the normal stress along the failure surface and σ z is the average of the vertical stress within the failure surface (frustrum). Because of difficulties in determining the failure surface and normal stress acting on it, the following empirical correlation for D F was presented from a linear regression to fit the centrifuge measurements as follows: The operative friction and dilation angle were recommended to be selected through an iterative process, because they are dependent on the stress level Bolton [24]. The following equation was presented to determine q peak dependent operative friction and dilation angles: where I D = relative density, Q is the particle crushing strength of sand, and φ cv is the critical state friction angle of sand. It was reported that Q = 10 for silica sand. Hu et al. [12] performed spudcan penetration tests on medium-loose sand overlying clay and compared qpeak. It was reported that qpeak is not mobilized at the surface, but at a depth of dpeak. To account for dpeak, the model of Lee et al. [11] was modified as illustrated in Figure 1b. The following modified equation was proposed: Additionally, the relationship of DF was recalibrated and optimized for a wide range of spudcan geometries and sand stiffness. It was reported that a non-linear correlation was found between DF and Hs/D, resulting in the following empirical equation: Both models were developed from centrifuge experiments with a narrow range of su ranging from 11 to 19.1 kPa (Table 1), because the focus of both studies was punch-through failure.

Numerical Model
Finite element analyses using the CEL approach were performed using Abaqus/Explicit [25] to simulate penetrations of spudcan foundations with conical shapes with extruded spigots, the dimensions of which are reported in Hu et al. [12]. A 3D finite element model with an axisymmetric quarter domain was used for both the spudcan and soil, as shown in Figure 2. The spudcan was discretized with Lagrangian mesh using hexahedron elements (C3D8R). A rigid body constraint was assigned to prohibit relative deformation of the spudcan. The soil was modelled using the Eulerian domain hexahedron elements with reduced integration and hourglass control (EC3D8R). In the analyses, the soil mesh remained intact, whereas the soil materials were allowed to flow into or out of the Eulerian element. A layer of void elements above the top surface of sand layer was used to allow heaving of soil. The bottom of the soil domain was fixed in all translational degrees of freedom. Hu et al. [12] performed spudcan penetration tests on medium-loose sand overlying clay and compared q peak . It was reported that q peak is not mobilized at the surface, but at a depth of d peak .
To account for d peak , the model of Lee et al. [11] was modified as illustrated in Figure 1b. The following modified equation was proposed: Additionally, the relationship of D F was recalibrated and optimized for a wide range of spudcan geometries and sand stiffness. It was reported that a non-linear correlation was found between D F and H s /D, resulting in the following empirical equation: Both models were developed from centrifuge experiments with a narrow range of s u ranging from 11 to 19.1 kPa (Table 1), because the focus of both studies was punch-through failure.

Numerical Model
Finite element analyses using the CEL approach were performed using Abaqus/Explicit [25] to simulate penetrations of spudcan foundations with conical shapes with extruded spigots, the dimensions of which are reported in Hu et al. [12]. A 3D finite element model with an axisymmetric quarter domain was used for both the spudcan and soil, as shown in Figure 2. The spudcan was discretized with Lagrangian mesh using hexahedron elements (C3D8R). A rigid body constraint was assigned to prohibit relative deformation of the spudcan. The soil was modelled using the Eulerian domain hexahedron elements with reduced integration and hourglass control (EC3D8R). In the analyses, the soil mesh remained intact, whereas the soil materials were allowed to flow into or out of the Eulerian element. A layer of void elements above the top surface of sand layer was used to allow heaving of soil. The bottom of the soil domain was fixed in all translational degrees of freedom. The horizontal displacement constraints were applied at the lateral boundaries. The size of the computational domain is consistent with the Hu et al. [17] model. A uniform fine mesh for 1D around spudcan periphery was used, whereas a gradually increasing coarser mesh was extended to edges of the domain. The size of the smallest element used was 0.025D, which was selected from a sensitivity study, as displayed in Figure 3a.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 18 The horizontal displacement constraints were applied at the lateral boundaries. The size of the computational domain is consistent with the Hu et al. [17] model. A uniform fine mesh for 1D around spudcan periphery was used, whereas a gradually increasing coarser mesh was extended to edges of the domain. The size of the smallest element used was 0.025D, which was selected from a sensitivity study, as displayed in Figure 3a. The interface between the spudcan foundation and soil was set to follow the Coulomb friction law. The friction coefficient was defined as αtanϕ′, where α is the roughness factor of the footing. The reported value of α ranges from 0.3 to 0.6 [1, 14,15,26]. α was set to 0.5, following the recommendation of Hu et al. [16]. A sensitivity study was performed to quantify the influence of the penetration rate. The response was shown to converge at rates lower than 0.3 m/s, as shown in Figure 3b. The penetration rate used was 0.2 m/s, which is consistent with the value applied in previous studies [16,17]. Linear elastic-perfectly plastic models following the Mohr-Coulomb and Tresca failure criteria were used for sands and clays, respectively. A drained behavior was assumed for the sand layer, whereas clay was assumed to deform in an undrained condition. The computational process is composed of two steps. In the first step, the gravity load is applied to the model to induce geostatic stresses. In the second step, a constant penetration velocity is applied at the load reference point (RP) of the spudcan. The penetration resistance is calculated by dividing the vertical reaction force at RP with the cross-section area of the spudcan. The interface between the spudcan foundation and soil was set to follow the Coulomb friction law. The friction coefficient was defined as αtanφ , where α is the roughness factor of the footing. The reported value of α ranges from 0.3 to 0.6 [1, 14,15,26]. α was set to 0.5, following the recommendation of Hu et al. [16]. A sensitivity study was performed to quantify the influence of the penetration rate. The response was shown to converge at rates lower than 0.3 m/s, as shown in Figure 3b. The penetration rate used was 0.2 m/s, which is consistent with the value applied in previous studies [16,17]. Linear elastic-perfectly plastic models following the Mohr-Coulomb and Tresca failure criteria were used for sands and clays, respectively. A drained behavior was assumed for the sand layer, whereas clay was assumed to deform in an undrained condition. The computational process is composed of two steps. In the first step, the gravity load is applied to the model to induce geostatic stresses. In the second step, a constant penetration velocity is applied at the load reference point (RP) of the spudcan. The penetration resistance is calculated by dividing the vertical reaction force at RP with the cross-section area of the spudcan.

Validation of Model
Hu et al. [12] performed centrifuge experiments on medium-loose sand overlying normally consolidated clay with sand thickness H s of 6, 5, and 3.2 m and H s /D ratio between 0.16 and 1.0. Hu et al. [16] performed numerical simulations to model the complete penetration resistance profile of a spudcan in sand overlying clay. Two experimental penetration resistance-depth outputs (L1SP1 and L3SP1) measured from centrifuge experiments performed on medium-loose sand overlying normally consolidated clay were used to validate the numerical model. H s /D for L1SP1 and L3SP1 were 1.0 and 0.53, respectively. In the numerical simulations, the properties of the soil were assigned following those provided in Hu et al. [16] were used. The elastic modulus (E) and Poisson's ratio (ν) of sand were set to 25 MPa and 0.3, respectively. φ of sand was set to 31 • , as reported in [16]. s u was set to 12.96 kPa, and the strength gradient k was set to 1.54 kPa/m. The centrifuge test outputs and numerical simulations results are compared in Figure 4. Overall, good agreement between the experiment recordings and CEL outputs are shown. As for the L3SP1 test, an overestimation of the penetration resistance in the clay layer is revealed. However, the peak resistances in the sand layer predicted in both tests show favorable correspondences.

Validation of Model
Hu et al. [12] performed centrifuge experiments on medium-loose sand overlying normally consolidated clay with sand thickness Hs of 6, 5, and 3.2 m and Hs/D ratio between 0.16 and 1.0. Hu et al. [16] performed numerical simulations to model the complete penetration resistance profile of a spudcan in sand overlying clay. Two experimental penetration resistance-depth outputs (L1SP1 and L3SP1) measured from centrifuge experiments performed on medium-loose sand overlying normally consolidated clay were used to validate the numerical model. Hs/D for L1SP1 and L3SP1 were 1.0 and 0.53, respectively. In the numerical simulations, the properties of the soil were assigned following those provided in Hu et al. [16] were used. The elastic modulus (E) and Poisson's ratio (ν) of sand were set to 25 MPa and 0.3, respectively. ϕʹ of sand was set to 31°, as reported in [16]. su was set to 12.96 kPa, and the strength gradient k was set to 1.54 kPa/m. The centrifuge test outputs and numerical simulations results are compared in Figure 4. Overall, good agreement between the experiment recordings and CEL outputs are shown. As for the L3SP1 test, an overestimation of the penetration resistance in the clay layer is revealed. However, the peak resistances in the sand layer predicted in both tests show favorable correspondences.  Hu and Cassidy [21] performed centrifuge experiments on medium dense sand overlying stiff clays with a sand thickness Hs of 6 and 4 m, and Hs/D ranging from 0.33 to 1.0. Two experimental penetration resistance profiles (SP3 and SP6) measured from centrifuge experiments performed on medium dense sand overlying stiff clay were used to further validate the numerical model. In the numerical simulations, the value of E and ν of sand were set to 50 MPa and 0.3, respectively [16,26]. ϕʹ of sand was set to 37°, as reported in [21]. The values of su were set to 43.9 and 32.6 kPa, for SP3 and SP6, respectively. The strength gradient k was set to 3.4 kPa/m. The centrifuge test outputs and numerical simulations results are compared in Figure 5. Overall, good agreement between the experiment recordings and CEL outputs are shown. The peak resistances in the sand layer predicted in both tests show favorable correspondences. Hu and Cassidy [21] performed centrifuge experiments on medium dense sand overlying stiff clays with a sand thickness H s of 6 and 4 m, and H s /D ranging from 0.33 to 1.0. Two experimental penetration resistance profiles (SP3 and SP6) measured from centrifuge experiments performed on medium dense sand overlying stiff clay were used to further validate the numerical model. In the numerical simulations, the value of E and ν of sand were set to 50 MPa and 0.3, respectively [16,26]. φ of sand was set to 37 • , as reported in [21]. The values of s u were set to 43.9 and 32.6 kPa, for SP3 and SP6, respectively. The strength gradient k was set to 3.4 kPa/m. The centrifuge test outputs and numerical simulations results are compared in Figure 5. Overall, good agreement between the experiment recordings and CEL outputs are shown. The peak resistances in the sand layer predicted in both tests show favorable correspondences.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 18 Figure 5. Comparison of penetration resistance profiles measured from centrifuge tests (sand over stiff clay) and calculated from FE analyses.

Peak Penetration Resistance in Sand Overlying Clay
A parametric study was performed to investigate the influence of various factors on the bearing capacity of spudcan. A total of 120 FE simulations were performed, the matrix of which is summarized in Table 2. The unit weights of sand and clay were set to 11 kN/m 3 and 7.5 kN/m 3 , respectively. For the sand layer, E and ν were set to 25 MPa and 0.3, respectively. For the clay layer, ν = 0.49 and E was set to 500 su [16]. ϕ′ of sand assigned were = 31° and 35°. The dilation angle was calculated from the following empirical formulation [9,27,28]: The undrained strength of clay at the interface of sand-clay were set to 10, 20, 30, 40, and 60 kPa. The strength gradient k was set to 2 kPa/m. Four sizes were used for the spudcan foundation, where D was varied from 6 to 14 m. The protruded spigot and conical angles of spudcan were fixed to 76° and 13°, respectively. The sand layer height was varied from 4 to 14 m. Hs/D ranges from 1.0 to 0.28. qpeak were defined as the peak penetration resistance calculated within the sand layer regardless of penetration depth, because the results demonstrated that in stiff clays, the peak resistance does not always develop at a depth of 0.12 Hs as reported in Hu et al. [12]. Figure 6 displays the plastic shear strain increment contours calculated at qpeak for Hs/D = 0.7 and various combinations of ϕ′ and su. Comparisons show that an increase in su produces corresponding

Peak Penetration Resistance in Sand Overlying Clay
A parametric study was performed to investigate the influence of various factors on the bearing capacity of spudcan. A total of 120 FE simulations were performed, the matrix of which is summarized in Table 2. The unit weights of sand and clay were set to 11 kN/m 3 and 7.5 kN/m 3 , respectively. For the sand layer, E and ν were set to 25 MPa and 0.3, respectively. For the clay layer, ν = 0.49 and E was set to 500 s u [16]. φ of sand assigned were = 31 • and 35 • . The dilation angle was calculated from the following empirical formulation [9,27,28]: The undrained strength of clay at the interface of sand-clay were set to 10,20,30,40, and 60 kPa. The strength gradient k was set to 2 kPa/m. Four sizes were used for the spudcan foundation, where D was varied from 6 to 14 m. The protruded spigot and conical angles of spudcan were fixed to 76 • and 13 • , respectively. The sand layer height was varied from 4 to 14 m. H s /D ranges from 1.0 to 0.28. q peak were defined as the peak penetration resistance calculated within the sand layer regardless of penetration depth, because the results demonstrated that in stiff clays, the peak resistance does not always develop at a depth of 0.12 H s as reported in Hu et al. [12].  Figure 6 displays the plastic shear strain increment contours calculated at q peak for H s /D = 0.7 and various combinations of φ and s u . Comparisons show that an increase in s u produces corresponding change in the slope of the failure surface for both φ = 31 • and 35 • cases. Whereas the failure surface are tilted outwards with respect to the centerline for s u = 10 kPa, they are inclined inwards for the case of s u = 60 kPa. This shift produces the area of influence for the clay to be smaller for stiff clay relative to the soft clay. For s u = 10 kPa, the lateral spreading of the soft clay can also be observed. It is shown to influence the trace of the plastic strain within the sand layer. Profiles with larger φ produce higher levels of shear strain in the sand layer, displaying that the mobilized resistance in the sand layer increases with an increase its shear strength.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 18    Figures 6 and 7 demonstrate that q peak is influenced by H s /D, φ , and s u . As for φ , the influence is not as significant compared with the parameters H s /D and s u , because only loose to medium dense sand with φ ranging from 31 to 35 • is considered in this study. Figure 8 shows the results of all calculated q peak . q peak is revealed to be positively correlated with H s /D, D, φ , and s u . The calculated q peak are compared with those estimated by theoretical models of Lee et al. [11] and Hu et al. [12] in Figure 9.   It should be noted that both the models of Lee et al. [11] and Hu et al. [12] require an iterative procedure to determine the operative friction and dilation angles, as summarized in the previous section. Additionally, the relative density of the sand layer needs to be defined. In this evaluation, the relative densities for φ = 31 • and 35 • were defined as 44% and 60%, respectively, based on the empirical correlation of Mujtaba et al. [29]. The operative friction angles determined from this iterative procedure vary from 33.5 • to 36.1 • and 32.1 • to 34.1 • for φ = 31 • and 35 • cases, respectively.
The comparisons show that s u has an important influence on the predicted q peak . It is also demonstrated that the theoretical models provide a reasonable estimate up to s u = 20 kPa, which fall within the range of data from which the models were developed. This overestimation is caused by the empirical correlations for D F , which are only dependent on H s /D. q peak increases linearly with s u , thus not capturing the reduction in the mobilized area of the bottom clay layer which produces a lower capacity.  Following the approaches of [11,12], an empirical correlation for D F was developed through a trial-and-error process. D F was varied until q peak predicted by the theoretical model of Hu et al. [12] most favorably fit the calculated output. It should be noted that the iterative procedure to determine the operative friction and dilation angles were not used, because it requires an additional assumption for the relative density, which is most often not provided. Instead, the friction and dilation angles used in the simulations were applied as operative angles, because comparisons demonstrate that this additional step has marginal influence on the calculated response for the cases used in this study. The back-calculated D F values are plotted against H s /D for the range of selected s u in Figure 10. The D F relationship proposed by Hu et al. [12] is also plotted in Figure 10 for comparison purposes. It is demonstrated that D F is strongly dependent on both H s /D and s u . D F is shown to be negatively correlated to both H s /D and s u . For all s u values, the variation in D F is more significant for H s /D < 0.5. It is also illustrated that the relationships are non-linear. A new relationship for D F accounting for the influence of s u in addition to H s /D is as follows: where p a is the atmospheric pressure. The equation is applicable for the range applied in this study, which is 0.28 ≤ H s /D ≤ 1.0.
It should be noted that both the models of Lee et al. [11] and Hu et al. [12] require an iterative procedure to determine the operative friction and dilation angles, as summarized in the previous section. Additionally, the relative density of the sand layer needs to be defined. In this evaluation, the relative densities for ϕ′ = 31° and 35° were defined as 44% and 60%, respectively, based on the empirical correlation of Mujtaba et al. [29]. The operative friction angles determined from this iterative procedure vary from 33.5° to 36.1° and 32.1° to 34.1° for ϕ′ = 31° and 35° cases, respectively.
The comparisons show that su has an important influence on the predicted qpeak. It is also demonstrated that the theoretical models provide a reasonable estimate up to su = 20 kPa, which fall within the range of data from which the models were developed. This overestimation is caused by the empirical correlations for DF, which are only dependent on Hs/D. qpeak increases linearly with su, thus not capturing the reduction in the mobilized area of the bottom clay layer which produces a lower capacity.
Following the approaches of [11,12], an empirical correlation for DF was developed through a trial-and-error process. DF was varied until qpeak predicted by the theoretical model of Hu et al. [12] most favorably fit the calculated output. It should be noted that the iterative procedure to determine Figure 9. Comparison of calculated and estimated peak bearing resistance: (a,b) stress-dependent model Lee et al. [11], (c,d) modified stress model Hu et al. [12].
The proposed equation is compatible with the q peak model of Hu et al. [12]. Figure 11 compares the predicted peak penetration resistances with the numerically calculated outputs. The residuals between the numerically calculated and the predicted peak resistances are within 20%. The comparisons highlight that the proposed equation provides reliable estimates of peak penetration resistance even for stiff clays with s u up to 60 kPa.
The proposed equation is also compared with the Hu et al. [30] D F relationship. Hu et al. [30] proposed an analytical model for calculation of the post-peak penetration resistance in sand over clay soils using measurements from centrifuge tests. The post-peak failure was reported to occur at a depth of d post-peak = 0.3 H s . A database of 66 centrifuge tests was used to calibrate a new D F relationship dependent on H s /D: where d is the penetration depth. The applicable range of d is from 0.3H s to H s . The reported equation is not dependent on s u . Figure 12 compares the calculated values of D F using Equation (12) where pa is the atmospheric pressure. The equation is applicable for the range applied in this study, which is 0.28 ≤ Hs/D ≤ 1.0. The proposed equation is compatible with the qpeak model of Hu et al. [12]. Figure 11 compares the predicted peak penetration resistances with the numerically calculated outputs. The residuals between the numerically calculated and the predicted peak resistances are within 20%. The comparisons highlight that the proposed equation provides reliable estimates of peak penetration resistance even for stiff clays with su up to 60 kPa.  The proposed equation is also compared with the Hu et al. [30] DF relationship. Hu et al. [30] proposed an analytical model for calculation of the post-peak penetration resistance in sand over clay soils using measurements from centrifuge tests. The post-peak failure was reported to occur at a depth of dpost-peak = 0.3 Hs. A database of 66 centrifuge tests was used to calibrate a new DF relationship dependent on Hs/D: where d is the penetration depth. The applicable range of d is from 0.3Hs to Hs. The reported equation is not dependent on su. Figure 12 compares the calculated values of DF using Equation (12) and the empirical equation of Hu et al. [30]. The equation of Hu et al. [30] approximately falls between the calculated values for su = 30 and 40 kPa. Because the equation is not dependent on su, it does not capture the range of factors from soft to stiff clays well.

Conclusions
The objective of this paper was to present a procedure to estimate qpeak of spudcan foundations in sand over both soft and stiff clay profiles. A comprehensive suite of 3D FE analyses incorporating the CEL approach were performed. An elastic-perfectly plastic model with non-associated flow rule following the Mohr-Coulomb failure criterion was used for the sand layer. The same model following the Tresca failure criterion was used for the clay layer. The numerical model was validated against the centrifuge tests performed and reported in a previous study.
A broad range of parameters including Hs, D, ϕ′ and su were employed in the analyses to calculate qpeak. The plastic shear strains that develop in the sand and clay layers show that the penetration response is critically influenced by both Hs/D and su. The area of influence of the clay layer

Conclusions
The objective of this paper was to present a procedure to estimate q peak of spudcan foundations in sand over both soft and stiff clay profiles. A comprehensive suite of 3D FE analyses incorporating the CEL approach were performed. An elastic-perfectly plastic model with non-associated flow rule following the Mohr-Coulomb failure criterion was used for the sand layer. The same model following the Tresca failure criterion was used for the clay layer. The numerical model was validated against the centrifuge tests performed and reported in a previous study.
A broad range of parameters including H s , D, φ and s u were employed in the analyses to calculate q peak . The plastic shear strains that develop in the sand and clay layers show that the penetration response is critically influenced by both H s /D and s u. The area of influence of the clay layer increases with a decrease in H s /D. An increase in s u produces a corresponding change in the slope of the failure surface. Whereas the failure surface is tilted outwards with respect to the centerline for soft clay, it is inclined inwards for the case of stiff clay. This shift causes the area of influence of the clay layer to be smaller for the stiff clay relative to the soft clay.
The calculated results are compared with two widely used analytical models. It is demonstrated that the previously proposed models provide consistent overestimations beyond s u = 20 kPa, which is outside the range used to calibrate these models. The main reason for this misfit is the independence of the key parameter D F of these models on s u . The calculated outputs are used to develop a new predictive equation for D F that is conditioned on both H s /D and s u . Using the modified equation for D F and plugging it in one of the available analytical models, favorable fits are achieved for the entire range of s u covering both soft and stiff clays. Funding: This research was a part of the project titled 'Development and demonstration of a decommission device for piles in shallow sea following international specifications', funded by the Ministry of Oceans and Fisheries, Korea.

Conflicts of Interest:
The authors declare no conflict of interest.