A New p–y Curve for Laterally Loaded Large-Diameter Monopiles in Soft Clays

: In harsh offshore environmental conditions, the monopile foundations supporting offshore wind turbines must be designed for lateral loads such as winds, waves, and currents. The Beam on Nonlinear Winkler Foundation (BNWF) method has been widely used because of its clear concept and lower calculation cost. The selection of a reasonable p–y curve is critical to the calculation accuracy of this method. This paper clariﬁed the defects of widely used API p–y curves for soft clays and then proposed a new p–y curve with better versatility and applicability. The suitability of the proposed p–y curve was validated by comparing it with the calculation results from the three-dimensional ﬁnite element method (3D FEM). Compared with the API p–y curve, the proposed p–y curve can better predict the lateral behavior of large-diameter piles in soft clays, such as the load–displacement curve of the pile head, lateral deﬂection proﬁle, and bending moment proﬁle. The research ﬁndings can provide guidance for the design of monopile foundations supporting offshore wind turbines in soft clays.


Introduction
Natural resources such as coal, petroleum, and natural gas are being consumed at an alarming rate, hence many countries are actively promoting green and sustainable new energy technologies.Wind energy is one of the most promising new energies due to its advantages of cleanness and renewability.In addition, compared with land-based wind energy, offshore wind energy is attractive for its various advantages including high speed, low wind shear, low turbulence, and high output.The urgent demand for wind power is promoting the development of construction technology for offshore wind turbines (OWTs) [1].The design of the foundations supporting OWTs is extremely important because it may account for as much as 30% of the total cost of a typical offshore wind project [2,3].Nonetheless, the large diameter monopile foundation remains by far the most popular foundation type being deployed in more than 75% of the currently installed OWTs [4,5].In harsh offshore environmental conditions, pile foundations must be designed for lateral loads because OWTs are always subjected to such loads caused by winds, waves, and currents.
The Beam on Nonlinear Winkler Foundation (BNWF) method based on the p-y curve is widely used to analyze the behavior of laterally loaded piles.McClelland and Focht [6] first proposed the concept of a p-y curve by comparing the measured curve of lateral soil resistance versus the pile lateral deflection with the soil stress-strain curve in the consolidated undrained tests.Matlock [7] proposed a p-y curve model for piles installed in clay based on limited pile load tests at the Sabine River site, which was adopted by the API RP 2GEO (American Petroleum Institute (API)) [8] and is widely used in industry.In the BNWF approach, the pile is simulated as an elastic beam while the soil is idealized as a series of discrete nonlinear springs along the depth of the pile.The nonlinear p-y springs describe the relationship between the local lateral soil resistance p and the pile lateral relative displacement y.The BNWF method based on the p-y curve is actually a simplified method ignoring the soil continuity.
In recent decades, some new p-y curves for soft clay have been proposed to analyze the lateral behavior of piles.The hyperbolic equation, initially proposed by Kondner [9] and later identified by Georgiadis et al. [10] as a suitable formulation, is used to define the shape of the p-y curves.Dewaikar and Patil et al. [11] developed an improved method for the construction of hyperbolic p-y curves, in which the initial stiffness of the p-y curve is taken to be increasing with depth.Jeanjean [12] proposed a hyperbolic tangent p-y curve and verified its rationality using the centrifugal test and finite element simulation results of pile in clays.Klinkvort and Hededal [13] performed a series of centrifuge tests to investigate the behavior of a rigid pile loaded with a high eccentricity and presented a reformulation of the p-y curve.Rathod et al. [14] developed new p-y curves for piles located on crests of soft clay with different sloping ground surfaces under static lateral loading by carrying out a series of laboratory model tests.Zhu et al. [15] modified the hyperbolic tangent p-y curve by introducing an additional constant and verified its rationality using the field tests of offshore driven piles.Xu et al. [16] found that the hyperbolic method yields larger prediction errors than the API method and proposed a modified p-y approach based on field measurements.The p-y curve method has been developed to better simulate lateral responses based on model tests and numerical simulation results [17][18][19][20][21][22][23][24].The above research mainly focus on the p-y curve under lateral monotonic loads.
The pile cyclic lateral response was also analyzed by developing cyclic p-y hysteresis curves based on the Beam on Nonlinear Winkler Foundation (BNWF) method.In these models, some pile-soil interaction characteristics are simulated using specific components such as spring, dashpot, drag cell, and gap cell to capture the dynamic hysteretic behavior of the soil-pile interface [25][26][27][28][29].However, the introduction of too many components usually makes the implementation of the model more complex.In addition to the p-y method, the three-dimensional finite element method (3D FEM) was also used to analyze the lateral response of piles even though it requires larger costs of calculation [30][31][32][33][34][35][36][37][38].The 3D FEM can better simulate the pile-soil interaction because both the pile and soil are simulated by continuum elements.Compared with the method, although discrete the p-y springs do not rigorously capture the soil continuum behavior, the BNWF method based on the p-y curve has still been widely utilized because of its clear concept and less calculation cost.
At present, the API p-y curve for soft clay [8], which is most widely used in industry, is proposed based on the specific field test of small diameter piles.Its applicability to larger diameter pile foundations, such as large-diameter steel pipe piles for offshore wind turbines, has been questioned.P-y curves with better applicability urgently needs to be developed.This research has two main objectives: (a) clarify the defects of widely used API p-y curve for soft clays and (b) develop a new p-y curve of soft clay with better applicability for large diameter monopile foundations.
In order to achieve the above research purposes, the BNWF method based on the p-y curve and three-dimensional finite element method for analyzing the lateral responses of monopile are developed based on ABAQUS software; then, the accuracy of the proposed two methods are verified using the field test results of pile in soft clays [7].Subsequently, the lateral responses of monopiles with different diameters are simulated employing the two methods.The defects of the API p-y curve are clarified by comparing the calculated results obtained by the two methods.Based on the above findings, a new p-y model for soft clay is proposed and the good applicability of the model to large diameter monopiles is validated by the above three-dimensional numerical simulation results.The novelty of this study is that a new p-y model with variable parameters was proposed, which is more versatile and offers a wider applicability to represent different soft clay and pile diameters.
Compared with the API p-y curve, the proposed p-y curve can better predict the lateral behavior of large diameter monopiles in soft clays.

Beam on Nonlinear Winkler Foundation (BNWF) Method
The BNWF model is developed based on ABAQUS software [39] to simulate the lateral monotonic load tests on teel pipe piles installed in soft clay at a site near the Sabine River [7].The test pile had an outer diameter of 0.324 m with a wall thickness of 12.7 mm.The total length of the pile was 13.1 m with an embedded depth of 12.8 m.The lateral load was applied at 0. pile diameters.Compared with the API p-y curve, the proposed p-y curve can better predict the lateral behavior of large diameter monopiles in soft clays.

Beam on Nonlinear Winkler Foundation (BNWF) Method
The BNWF model is developed based on ABAQUS software [39] to simulate the lateral monotonic load tests on teel pipe piles installed in soft clay at a site near the Sabine River [7].The test pile had an outer diameter of 0.324 m with a wall thickness of 12.7 mm.The total length of the pile was 13.1 m with an embedded depth of 12.8 m.The lateral load was applied at 0.  where p is the mobilized lateral resistance (kPa); y is the local pile lateral displacement (m); y 50 is the pile lateral displacement at half of ultimate soil resistance (m); and y 50 = 2.5ε 50 D, ε 50 is the strain at one-half the maximum deviator stress in laboratory undrained compression tests of undisturbed soil samples.ε 50 = 0.02 in this study.D is the pile diameter and D = 0.324 m. γ is effective unit weight of soil (kN/m 3 ) and γ = 7.85 kN/m 3 .p u is the ultimate resistance (kPa); p u = 3s u + γz + Js u z/D for z < z r and p u = 9s u for z ≥ z r .s u is the undrained shear strength (kPa) and s u = 14.4 kPa; z is the depth below the original seafloor (m); J is the dimensionless empirical constant with values ranging from 0.25 to 0.5 having been determined using field testing, and J = 0.5 in this study.Figure 2a,b compare the calculated lateral deflection and bending moment profiles along the pile and those evaluated from the field measurements for four different pile head load levels.It can be seen that the calculated responses agree well with the test results for all load levels despite the slight deviations.One reason for the slight deviation is that the soil strength is idealized as uniform strength, but it is not uniform in the field test.Another reason is the limitation of the BNWF method itself, i.e., it ignores the continuity of soil, which will cause the error of calculation.This good agreement demonstrates that the developed BNWF model based on ABAQUS software can well predict the increase in lateral deflection and bending moment profiles with the increasing load level and the location of the maximum bending moment for different load levels.It also demonstrates the ability of the developed model to predict the lateral monotonic behavior of the pile in soft clays.
the pile diameter and D = 0.324 m.  is effective unit weight of soil (kN/m 3 ) and  = z is the depth below the original seafloor (m); J is the dimensionless empirical constant with values ranging from 0.25 to 0.5 having been determined using field testing, and J = 0.5 in this study.
Figure 2a,b compare the calculated lateral deflection and bending moment profiles along the pile and those evaluated from the field measurements for four different pile head load levels.It can be seen that the calculated responses agree well with the test results for all load levels despite the slight deviations.One reason for the slight deviation is that the soil strength is idealized as uniform strength, but it is not uniform in the field test.Another reason is the limitation of the BNWF method itself, i.e., it ignores the continuity of soil, which will cause the error of calculation.This good agreement demonstrates that the developed BNWF model based on ABAQUS software can well predict the increase in lateral deflection and bending moment profiles with the increasing load level and the location of the maximum bending moment for different load levels.It also demonstrates the ability of the developed model to predict the lateral monotonic behavior of the pile in soft clays.

Three-Dimensional Finite Element Method (3D FEM)
A three-dimensional finite element model is established by employing the finite element software ABAQUS to simulate the interaction between soft clay and the pile foundation considering the specific conditions for Matlock's filed tests.Exploiting the symmetries of the geometry and loading conditions, the model simulates only one half of the system, as shown in Figure 3, to improve the computational efficiency.A sensitivity study indicates that a model of size 20 D × 50 D (the cross section is a semicircle and D is the pile diameter) is sufficient to avoid boundary effects on the simulation results.The pile and

Three-Dimensional Finite Element Method (3D FEM)
A three-dimensional finite element model is established by employing the finite element software ABAQUS to simulate the interaction between soft clay and the pile foundation considering the specific conditions for Matlock's filed tests.Exploiting the symmetries of the geometry and loading conditions, the model simulates only one half of the system, as shown in Figure 3, to improve the computational efficiency.A sensitivity study indicates that a model of size 20 D × 50 D (the cross section is a semicircle and D is the pile diameter) is sufficient to avoid boundary effects on the simulation results.The pile and the soils are all simulated using 5466 8-node linear brick elements with reduced integrations (C3D8R).Normal horizontal constraints are applied to the vertical boundaries and fixed constraints are applied to the bottom boundary.The top boundary is fully free.Lateral loads are applied to the pile according to the loading condition in the above field tests.
the soils are all simulated using 5466 8-node linear brick elements with reduced integrations (C3D8R).Normal horizontal constraints are applied to the vertical boundaries and fixed constraints are applied to the bottom boundary.The top boundary is fully free.Lateral loads are applied to the pile according to the loading condition in the above field tests.The interaction between the monopile and surrounding soil is simulated utilizing the surface-to-surface contact method.The surface of the pile is defined as the master surface (relatively stiff) and the soil surface in contact with the pile is defined as the slave surface (relatively soft).Hard contact is set in the normal direction of the contact surface, which allows separation between the interface elements when they are subjected to tension.The penalty contact method is used in the tangential direction.When the two surfaces are in contact, the interface behavior is governed by Coulomb's friction theory.The critical friction shear stress crit τ at the contact surface can be expressed by crit c τ = μ p × in terms of the frictional coefficient μ and contact pressure c p , with μ = 0.3 in this study.When the shear stress at the contact surface exceeds crit τ , the tangential slip occurs.The behavior of clay is simulated by employing the elastic-perfectly plastic constitutive model with Mises yield criterion.The model parameters include the soil elastic modulus, Poisson's ratio, and shear strength, which can be determined using the laboratory tests or field tests.In this study, the model parameters are set according to the research of Matlock [7] and Hong et al. [37] as shown in Table 1.The material behavior of the hollow steel pipe pile is simulated using the linear elastic constitutive model.The model parameters are set according to the above steel properties.During the finite element simulation, geostatic stress balance analysis shall be carried out first, i.e., gravity is applied to the soil domain to establish the initial stress field.Meanwhile, the initial displacement field is set to 0. Then, lateral monotonic behavior of the pile is analyzed.The interaction between the monopile and surrounding soil is simulated utilizing the surface-to-surface contact method.The surface of the pile is defined as the master surface (relatively stiff) and the soil surface in contact with the pile is defined as the slave surface (relatively soft).Hard contact is set in the normal direction of the contact surface, which allows separation between the interface elements when they are subjected to tension.The penalty contact method is used in the tangential direction.When the two surfaces are in contact, the interface behavior is governed by Coulomb's friction theory.The critical friction shear stress τ crit at the contact surface can be expressed by τ crit = µ • p c in terms of the frictional coefficient µ and contact pressure p c , with µ = 0.3 in this study.When the shear stress at the contact surface exceeds τ crit , the tangential slip occurs.
The behavior of clay is simulated by employing the elastic-perfectly plastic constitutive model with Mises yield criterion.The model parameters include the soil elastic modulus, Poisson's ratio, and shear strength, which can be determined using the laboratory tests or field tests.In this study, the model parameters are set according to the research of Matlock [7] and Hong et al. [37] as shown in Table 1.The material behavior of the hollow steel pipe pile is simulated using the linear elastic constitutive model.The model parameters are set according to the above steel properties.During the finite element simulation, geostatic stress balance analysis shall be carried out first, i.e., gravity is applied to the soil domain to establish the initial stress field.Meanwhile, the initial displacement field is set to 0. Then, lateral monotonic behavior of the pile is analyzed.Figure 4 shows the comparison between the calculated results and test results of the lateral deflection and bending moment profiles along the pile for four different pile head load levels.The calculated results agree with the test results in the overall trend although there is some deviation.One reason for the deviation is that the soil strength is idealized as uniform strength, but it is not uniform in the field test.Another reason is that it is a simplification to use the elastic-perfectly plastic constitutive model to simulate the stress-strain responses of the soil; when more advanced constitutive models are adopted, better calculation results can be expected.This comparison demonstrates that the developed three-dimensional finite element model based on ABAQUS software can well predict the lateral monotonic behavior of the pile in soft clays.Figure 4 shows the comparison between the calculated results and test results of the lateral deflection and bending moment profiles along the pile for four different pile head load levels.The calculated results agree with the test results in the overall trend although there is some deviation.One reason for the deviation is that the soil strength is idealized as uniform strength, but it is not uniform in the field test.Another reason is that it is a simplification to use the elastic-perfectly plastic constitutive model to simulate the stressstrain responses of the soil; when more advanced constitutive models are adopted, better calculation results can be expected.This comparison demonstrates that the developed three-dimensional finite element model based on ABAQUS software can well predict the lateral monotonic behavior of the pile in soft clays.

Numerical Model
Two methods for analyzing horizontally loaded piles, namely the BNWF method and 3D finite element method, have been established and validated, respectively, in Section 2 above.It should be noted that the BNWF model is actually a two-dimensional simplified calculation model of the pile-soil interaction; the accuracy of the calculation results employing this model mainly depends on the rationality of the p-y curve, hence it has certain limitations when using the currently questionable API p-y curve.In comparison, the three-dimensional finite element method has a better applicability for different pile diameters and site soil characteristics, although the calculation cost is relative larger.In this chapter, the above two methods will be used to simulate the lateral behavior of monopiles with different diameters.Using 3D finite element simulation results as a relatively correct reference standard, the limitation of the current API p-y curve is clarified by comparing the calculation results obtained from two methods.

Numerical Model
Two methods for analyzing horizontally loaded piles, namely the BNWF method and 3D finite element method, have been established and validated, respectively, in Section 2 above.It should be noted that the BNWF model is actually a two-dimensional simplified calculation model of the pile-soil interaction; the accuracy of the calculation results employing this model mainly depends on the rationality of the p-y curve, hence it has certain limitations when using the currently questionable API p-y curve.In comparison, the three-dimensional finite element method has a better applicability for different pile diameters and site soil characteristics, although the calculation cost is relative larger.In this chapter, the above two methods will be used to simulate the lateral behavior of monopiles with different diameters.Using 3D finite element simulation results as a relatively correct reference standard, the limitation of the current API p-y curve is clarified by comparing the calculation results obtained from two methods.
Four diameter piles (1 m, 2 m, 3 m, 5 m) were selected to build both two-dimensional BNWF models and 3D finite element models.The total length of all piles is 60 m with an embedded depth of 50 m and the wall thickness is 30 mm.A lateral load of 2 MN is applied at the pile head (10 m above the mud surface).The three-dimensional finite element models of the four diameter piles are shown in Figure 5.In the 3D numerical model, the constitutive models and material parameters of pile and soil, pile-soil contact condition, boundary conditions, etc. are consistent with those in Section 2 above, which will not be repeated here due to space limitations.The setting of the two-dimensional model is also consistent with those in Section 2 except for the different pile sizes.
embedded depth of 50 m and the wall thickness is 30 mm.A lateral load of 2 MN is appl at the pile head (10 m above the mud surface).The three-dimensional finite element mo els of the four diameter piles are shown in Figure 5.In the 3D numerical model, the co stitutive models and material parameters of pile and soil, pile-soil contact conditi boundary conditions, etc. are consistent with those in Section 2 above, which will not repeated here due to space limitations.The setting of the two-dimensional model is a consistent with those in Section 2 except for the different pile sizes.

Comparison of Predicted Lateral Behavior for Two Methods
Figure 6 compares the lateral response of piles with various diameters obtained usi the BNWF method based on the API p-y curve and 3D FEM when the pile head lo equals 250 kN.It can be seen from Figure 6a that the nonlinearity of pile head load-d placement curve becomes more significant as the pile diameter decreases.The displa ment of the pile head is greater for the smaller pile diameter under the same lateral lo The prediction result of displacement based on the BNWF method is larger than that the 3D FEM, which indicates that the BNWF method underestimates the bearing capac of piles.As shown in Figure 6b, the lateral deflection of the pile along the depth (ran from −10 to 10) increases significantly as the pile diameter decreases, which is because pile becomes more flexible due to the decrease in the stiffness.The lateral deflection p file obtained by the BNWF method is larger than 3D FEM and the difference between

Comparison of Predicted Lateral Behavior for Two Methods
Figure 6 compares the lateral response of piles with various diameters obtained using the BNWF method based on the API p-y curve and 3D FEM when the pile head load equals 250 kN.It can be seen from Figure 6a that the nonlinearity of pile head load-displacement curve becomes more significant as the pile diameter decreases.The displacement of the pile head is greater for the smaller pile diameter under the same lateral load.The prediction result of displacement based on the BNWF method is larger than that of the 3D FEM, which indicates that the BNWF method underestimates the bearing capacity of piles.As shown in Figure 6b, the lateral deflection of the pile along the depth (range from −10 to 10) increases significantly as the pile diameter decreases, which is because the pile becomes more flexible due to the decrease in the stiffness.The lateral deflection profile obtained by the BNWF method is larger than 3D FEM and the difference between the two calculation results is more significant for the smaller diameter piles.It can be observed from Figure 6c that the bending moment of the pile along the depth increases significantly as the pile diameter increases and the location of the maximum bending moment moves downward gradually.The bending moment profile obtained by the BNWF method is larger than 3D FEM and the difference is more significant for the larger diameter piles.From the above comparison, it can be concluded that the BNWF method based on the API p-y curve will overestimate the deflection and bending moment and underestimate the ultimate bearing capacity of the large diameter pile.
, x FOR PEER REVIEW 8 of 17 pile becomes more flexible due to the decrease in the stiffness.The lateral deflection profile obtained by the BNWF method is larger than 3D FEM and the difference between the two calculation results is more significant for the smaller diameter piles.It can be observed from Figure 6c that the bending moment of the pile along the depth increases significantly as the pile diameter increases and the location of the maximum bending moment moves downward gradually.The bending moment profile obtained by the BNWF method is larger than 3D FEM and the difference is more significant for the larger diameter piles.
From the above comparison, it can be concluded that the BNWF method based on the API p-y curve will overestimate the deflection and bending moment and underestimate the ultimate bearing capacity of the large diameter pile.

Limitations of API p-y Curve for Soft Clay
The pile-soil interaction is simulated through the surface-surface contact mode for the 3D finite element model.The soil resistance value P at a certain depth of the pile shaft is calculated by extracting the contact force on the contact surface.Figure 7a shows the typical distribution of the calculated contact normal force (CNF) and the contact shear

Limitations of API p-y Curve for Soft Clay
The pile-soil interaction is simulated through the surface-surface contact mode for the 3D finite element model.The soil resistance value P at a certain depth of the pile shaft is calculated by extracting the contact force on the contact surface.Figure 7a shows the typical distribution of the calculated contact normal force (CNF) and the contact shear force (CSF) around the pile.Figure 7b shows the diagrammatic sketch for the x-direction component of the CNF (CNF1) and the x-direction component of the CSF (CSF1) at a certain position on the pile cross-section.Figure 7b,c shows the typical distribution of CNF1 and CSF1 around the pile.It can be observed that CNF1 is the largest in the x-axis direction and the smallest at the edges on both sides of the section (almost 0); on the contrary, CSF1 is the smallest (almost 0) in the x-axis direction and the largest at the edges on both sides of the section.
Sustainability 2022, 14, x FOR PEER REVIEW certain position on the pile cross-section.Figure 7b,c shows the typical distribu CNF1 and CSF1 around the pile.It can be observed that CNF1 is the largest in th direction and the smallest at the edges on both sides of the section (almost 0); on t trary, CSF1 is the smallest (almost 0) in the x-axis direction and the largest at the e both sides of the section.The soil resistance P at a certain depth of the pile shaft can be determined by t of CNF1 and CSF1 around the pile, as represented by: ( ) where n is the number of nodes outside the pile cross-section at a certain depth of shaft.Then, the p-y curve can be determined based on the numerical simulation r Figure 8 compares the API p-y curve with the p-y curve obtained using the 3 element method.These p-y curves correspond to different depths (depth = 1 m, 3 and 10 m) for various pile diameters (diameter = 1 m, 2 m, 3 m, and 5 m).For thes diameter piles, there are significant difference between the API p-y curve and curve obtained using the 3D finite element method.The comparison demonstrates API p-y curve significantly underestimates the ultimate soil reaction of soft clays degree of the underestimation increases as the soil depth increases; (2) the API psignificantly underestimates the initial stiffness of soft clays.This conclusion can plain why the BNWF method based on the API p-y curve overestimates the de The soil resistance P at a certain depth of the pile shaft can be determined by the sum of CNF1 and CSF1 around the pile, as represented by: where n is the number of nodes outside the pile cross-section at a certain depth of the pile shaft.Then, the p-y curve can be determined based on the numerical simulation results.Figure 8 compares the API p-y curve with the p-y curve obtained using the 3D finite element method.These p-y curves correspond to different depths (depth = 1 m, 3 m, 5 m, and 10 m) for various pile diameters (diameter = 1 m, 2 m, 3 m, and 5 m).For these larger diameter piles, there are significant difference between the API p-y curve and the p-y curve obtained using the 3D finite element method.The comparison demonstrates: (1) the API p-y curve significantly underestimates the ultimate soil reaction of soft clays and the degree of the underestimation increases as the soil depth increases; (2) the API p-y curve significantly underestimates the initial stiffness of soft clays.This conclusion can well explain why the BNWF method based on the API p-y curve overestimates the deflection and bending moment while it underestimates the ultimate bearing capacity of the large diameter pile.

New p-y Formula and Parameter Influence
In order to improve the shortcomings of the API p-y curve, a new p-y curve is proposed in this section, as shown in the formula: where max G is the maximum shear modulus and max G is usually considered to be a multiple of u s for soft clay according to the engineering experience (e.g., typical range of max G = from 300 u s to 1800 u s ); parameters A and B control the shape of the p-y curve.max p is the ultimate lateral soil resistance, which is calculated based on the expression proposed by Murff and Hamilton [40] for shear strength profiles approximately linearly increasing with depth, i.e.,:

New p-y Formula and Parameter Influence
In order to improve the shortcomings of the API p-y curve, a new p-y curve is proposed in this section, as shown in the formula: where G max is the maximum shear modulus and G max is usually considered to be a multiple of s u for soft clay according to the engineering experience (e.g., typical range of G max = from 300 s u to 1800 s u ); parameters A and B control the shape of the p-y curve.p max is the ultimate lateral soil resistance, which is calculated based on the expression proposed by Murff and Hamilton [40] for shear strength profiles approximately linearly increasing with depth, i.e.,: where s u is the undrained shear strength of soil at the point in question; N p is a nondimensional lateral bearing factor; γ is the soil submerged unit weight; z is the soil depth below the original seafloor; D is the pile outside the diameter; s u0 is the shear strength intercept at the seafloor; and s u1 is the rate of increase in shear strength with depth.
In order to investigate the influence of the relevant parameters on the p-y curve, the soil properties are assumed as follows: s u = 14.4 kPa (s u0 = 14.4 kPa, s u1 = 0 in Equation ( 4)); G max = 500 s u = 7.2 MPa; pile diameter D = 3 m; p max = 170 kPa according to Equation (4) when soil depth z = 15 m.
A series of A values are set to investigate its influence on the shape of the p-y curve (B = 1 for all cases).Figure 9a,b show the p-y curves and the normalized p-y curves (p is normalized by p max and y is normalized by D) for various A values.It can be observed that the decrease in the soil stiffness (tangent modulus) becomes slower as the displacement increases and the variation of the p-y from linear to nonlinear becomes more gradual and starts earlier as A increases.It is found that the p-y (or normalized p-y) curve has insignificant change when A is less than 0.01 and hence A usually ranges from 0.1 to 10 to fit the actual p-y curves.It is worth noting that the current model can be simplified into Kondner's hyperbolic p-y model [9] when A = 1, as represented by: where u s is the undrained shear strength of soil at the point in question; p N is a non- dimensional lateral bearing factor; '  is the soil submerged unit weight; z is the soil depth below the original seafloor; D is the pile outside the diameter; 0 u s is the shear strength intercept at the seafloor; and 1 u s is the rate of increase in shear strength with depth.
In order to investigate the influence of the relevant parameters on the p-y curve, the soil properties are assumed as follows: A series of A values are set to investigate its influence on the shape of the p-y curve (B = 1 for all cases).Figure 9a p = y For this case, the initial maximum tangent modulus of the p-y curve is B • G max /D.A series of B values are set to investigate its influence on the shape of the p-y curve (A = 1 for all cases).Figure 10a,b show the p-y curves and the normalized p-y curves for various B values.Contrary to the effect of parameter A on the p-y curve, the decrease in the soil stiffness (tangent modulus) becomes faster as the displacement increases, and the variation of p-y from linear to nonlinear becomes more abrupt as B increases.When parameter B exceeds 10, it has little effect on the p-y (or normalized p-y) curve and hence B usually ranges from 1 to 10 to fit the actual p-y curves.Obviously, the proposed p-y model with variable A and B is more versatile and offers wider applicability to represent different soft clays and pile diameters.
various B values.Contrary to the effect of parameter A on the p-y curve, the decrease in the soil stiffness (tangent modulus) becomes faster as the displacement increases, and the variation of p-y from linear to nonlinear becomes more abrupt as B increases.When parameter B exceeds 10, it has little effect on the p-y (or normalized p-y) curve and hence B usually ranges from 1 to 10 to fit the actual p-y curves.Obviously, the proposed p-y model with variable A and B is more versatile and offers wider applicability to represent different soft clays and pile diameters.

Effect of Pile Diameter and Soil Depth on p-y Curve
The soil layer parameters are set as in Section 4.1, then different D values are set to investigate the effect of pile diameter on the shape of the p-y curve (A = 1, and B = 3).Figure 11a

Effect of Pile Diameter and Soil Depth on p-y Curve
The soil layer parameters are set as in Section 4.1, then different D values are set to investigate the effect of pile diameter on the shape of the p-y curve (A = 1, and B = 3).Figure 11a,b show the p-y curves and the normalized p-y curves for various D values.It can be observed from Figure 11a that the decrease in the soil stiffness (tangent modulus) becomes slower as the displacement increases and the ultimate lateral soil resistance p max becomes smaller as D increases.However, the normalized p-y curves for different pile diameters are completely consistent, as shown in Figure 11b.
variation of p-y from linear to nonlinear becomes more abrupt as B increases.When parameter B exceeds 10, it has little effect on the p-y (or normalized p-y) curve and hence B usually ranges from 1 to 10 to fit the actual p-y curves.Obviously, the proposed p-y model with variable A and B is more versatile and offers wider applicability to represent different soft clays and pile diameters.

Effect of Pile Diameter and Soil Depth on p-y Curve
The soil layer parameters are set as in Section 4.1, then different D values are set to investigate the effect of pile diameter on the shape of the p-y curve (A = 1, and B = 3).Figure 11a Different z values are set to investigate the effect of soil depth on the shape of the p-y curve (A = 1, and B = 3).Figure 12a,b shows the p-y curves and the normalized p-y curves for various z values.As expected, the ultimate lateral soil resistance p max becomes larger and the increasing rate of p max decreases gradually as the soil depth increases, as shown in Figure 12a.However, the normalized p-y curves for different soil depths have an insignificant difference, as shown in Figure 12b.Hence, it can be approximately considered that the normalized p-y curve at different depths is unique.y curve (A = 1, and B = 3).Figure 12a,b shows the p-y curves and the normalized p-y curves for various z values.As expected, the ultimate lateral soil resistance max p becomes larger and the increasing rate of max p decreases gradually as the soil depth increases, as shown in Figure 12a.However, the normalized p-y curves for different soil depths have an insignificant difference, as shown in Figure 12b.Hence, it can be approximately considered that the normalized p-y curve at different depths is unique.

Validation of New p-y Formula
The normalized p-y curves corresponding to different depths for various pile diameters can be obtained according to the above 3D finite element analysis results, as shown in Figure 13.The proposed p-y curves (A = 1 and B = 3) and API p-y curves are plotted also in the figure.It can be observed that the proposed p-y curves can generally fit those obtained from 3D numerical simulation well.However, as a comparison, the API p-y curve significantly underestimates the soil stiffness.It should be noted that the proposed normalized p-y curve does not correspond to a specific soil depth, which is due to the fact that the soil depth has an insignificant effect on the shape of the normalized p-y curve, as described in Section 4.2 above.This comparison reflects the superiority of the proposed p-y curve over the API p-y curve.

Validation of New p-y Formula
The normalized p-y curves corresponding to different depths for various pile diameters can be obtained according to the above 3D finite element analysis results, as shown in Figure 13.The proposed p-y curves (A = 1 and B = 3) and API p-y curves are plotted also in the figure.It can be observed that the proposed p-y curves can generally fit those obtained from 3D numerical simulation well.However, as a comparison, the API p-y curve significantly underestimates the soil stiffness.It should be noted that the proposed normalized p-y curve does not correspond to a specific soil depth, which is due to the fact that the soil depth has an insignificant effect on the shape of the normalized p-y curve, as described in Section 4.2 above.This comparison reflects the superiority of the proposed p-y curve over the API p-y curve.The lateral responses of the piles with various diameters are predicted employing the BNWF method based on the proposed API p-y curve.The prediction results are compared with the 3D finite element results on Section 3.2, as shown in Figure 14.The prediction results including the load-displacement curve of the pile head, lateral deflection profile,  The lateral responses of the piles with various diameters are predicted employing the BNWF method based on the proposed API p-y curve.The prediction results are compared with the 3D finite element results on Section 3.2, as shown in Figure 14.The prediction results including the load-displacement curve of the pile head, lateral deflection profile, and bending moment profile agree well with the 3D finite element results.This comparison demonstrates that, compared with API p-y curve, the proposed p-y curve can better predict the lateral behavior of piles in soft clays and has better applicability to large-diameter piles.

Summary and Conclusions Main Conclusions
(1) The BNWF method based on the p-y curve and 3D FEM for analyzing the lateral responses of monopiles are developed based on ABAQUS software and the accuracy of the two methods is validated using the field test results of the piles in soft clays.(2) The lateral responses of monopiles with different diameters are simulated by employing the two methods.It is found that the BNWF method based on the API p-y curve overestimates the deflection and bending moment while it underestimates the ultimate bearing capacity of the large diameter pile.
(3) The API p-y curve significantly underestimates the ultimate soil reaction of soft clays and the degree of the underestimation increases as the soil depth increases; it also significantly underestimates the initial stiffness of the soft clays.These findings are supported by previous research results [12,15].(4) A new p-y model with variable parameters A and B was proposed that is more versatile and offers wider applicability to represent different soft clay and pile diameters.
The current model can be simplified into the hyperbolic p-y model when A =1.

Significance and Application
As the most popular foundation type for supporting offshore wind turbines, large diameter monopile foundations must be designed for lateral loads such as winds, waves, and currents.The Beam on Nonlinear Winkler Foundation (BNWF) method has been widely used because of its clear concept and lower calculation cost.Compared with the three-dimensional finite element method, the BNWF method is easier to implement numerically.This method can be realized using ABAQUS software or, alternatively, by using other special calculation software for pile foundation such as LPILE or by compiling simple MATLAB programs.Hence, it has better applicability in engineering practice.The selection of a reasonable p-y curve is critical to the calculation accuracy of this method.This paper proposes a new p-y curve with better versatility and applicability.Compared with the API p-y curve, the proposed p-y curve can better predict the lateral behavior of piles in soft clays, such as the load-displacement curve of the pile head, lateral deflection profile, and bending moment profile, and has better applicability to large-diameter piles.However, it should be noted that the proposed p-y curve is mainly applicable to soft clay and not to other types of soil such as stiff clays and sands.When applied to engineering practice, the proposed p-y model can play a good role in the design of monopile foundations supporting offshore wind turbines in soft clays.
3 m above the ground surface.The average undrained shear strength of clays s u = 14.4 kPa, effective unit gravity γ = 7.85 kN/m 3 , elastic modulus of pile material E P = 210 GPa, density ρ = 1820 kg/m 3 , and Poisson's ratio υ = 0.33 were all utilized.The schematic diagram of the BNWF model is shown in Figure 1.The pile is discretized by 128 two-dimensional beam column elements (B21 element) and the length of each element L e = 0.5 m.A nonlinear spring element is set at each element node to simulate the pile-soil interaction.The elastic modulus of the nonlinear spring is set according to the API p-y curve for soft clay.The spring resistance force F = p • L e • D (D is the outer diameter of the pile).The continuous variation of the elastic modulus of spring reflects the nonlinear relationship between the spring resistance force F (or lateral soil resistance) and the spring displacement y (the pile lateral relative displacement).The API p-y formula and relevant parameter values for this study are as follows: Sustainability 2022, 14, x FOR PEER REVIEW 3 of 17

Figure 1 .
Figure 1.Schematic diagram of BNWF model for laterally loaded pile in soft clays.

Figure 1 .
Figure 1.Schematic diagram of BNWF model for laterally loaded pile in soft clays.  

s
is the undrained shear strength (kPa) and u s= 14.4 kPa;

Figure 2 .
Figure 2. Comparison between calculated and measured test results based on the BNWF method.(a) Lateral deflection profile (b) Bending moment profile.

Figure 2 .
Figure 2. Comparison between calculated and measured test results based on the BNWF method.(a) Lateral deflection profile (b) Bending moment profile.

Figure 3 .
Figure 3. Three-dimensional finite element model for laterally loaded pile in soft clays.

Figure 3 .
Figure 3. Three-dimensional finite element model for laterally loaded pile in soft clays.

Figure 4 .
Figure 4. Comparison between calculated and measured test results based on 3D FEM.(a) Lateral deflection profile (b) Bending moment profile.

Figure 4 .
Figure 4. Comparison between calculated and measured test results based on 3D FEM.(a) Lateral deflection profile (b) Bending moment profile.

Figure 5 .
Figure 5. 3D finite element model for laterally loaded pile with different diameter in soft clays.

Figure 5 .
Figure 5. 3D finite element model for laterally loaded pile with different diameter in soft clays.

Figure 6 .
Figure 6.Comparison of lateral response of piles obtained by BNWF method (API p-y curve) and 3D FEM.(a) Load-displacement curve of pile head.(b) Lateral deflection profile.(c) Bending moment profile.

Figure 6 .
Figure 6.Comparison of lateral response of piles obtained by BNWF method (API p-y curve) and 3D FEM.(a) Load-displacement curve of pile head.(b) Lateral deflection profile.(c) Bending moment profile.

Figure 7 .
Figure 7.Typical distribution of the calculated contact normal force (CNF) and contact sh (CSF).

FFigure 7 .
Figure 7.Typical distribution of the calculated contact normal force (CNF) and contact shear force (CSF).

Figure 9 .Figure 9 .
Figure 9.Effect of parameter A on the proposed p-y curve.(a) p-y curve (b) normalized p-y curve.

Figure 10 .
Figure 10.Effect of parameter B on the proposed p-y curve.(a) p-y curve (b) normalized p-y curve.

Figure 11 .Figure 10 .
Figure 11.Effect of pile diameter D on the proposed p-y curve.(a) p-y curve (b) normalized p-y curve.

Figure 10 .
Figure 10.Effect of parameter B on the proposed p-y curve.(a) p-y curve (b) normalized p-y curve.

Figure 11 .Figure 11 .
Figure 11.Effect of pile diameter D on the proposed p-y curve.(a) p-y curve (b) normalized p-y curve.

Figure 12 .
Figure 12.Effect of soil depth z on the proposed p-y curve.(a) p-y curve (b) normalized p-y curve.

Figure 12 .
Figure 12.Effect of soil depth z on the proposed p-y curve.(a) p-y curve (b) normalized p-y curve.

Figure 14 .
Figure 14.Comparison of the lateral responses of piles obtained using the BNWF method (proposed p-y curve) and 3D finite element method.(a) Load-displacement curve of pile head.(b) Lateral deflection profile (c) Bending moment profile.

( 1 )
The BNWF method based on the p-y curve and 3D FEM for analyzing the lateral responses of monopiles are developed based on ABAQUS software and the accuracy of the two methods is validated using the field test results of the piles in soft clays.(2)The lateral responses of monopiles with different diameters are simulated by employing the two methods.It is found that the BNWF method based on the API p-y curve overestimates the deflection and bending moment while it underestimates the ultimate bearing capacity of the large diameter pile.(3) The API p-y curve significantly underestimates the ultimate soil reaction of soft clays and the degree of the underestimation increases as the soil depth increases; it also significantly underestimates the initial stiffness of the soft clays.These findings are supported by previous research results[12,15].

Figure 14 .
Figure 14.Comparison of the lateral responses of piles obtained using the BNWF method (proposed p-y curve) and 3D finite element method.(a) Load-displacement curve of pile head.(b) Lateral deflection profile (c) Bending moment profile.

Table 1 .
Constitutive model parameters and relevant calculation parameters.