Scour E ﬀ ects on the Lateral Behavior of a Large-Diameter Monopile in Soft Clay: Role of Stress History

: Scouring of soil around large-diameter monopile will alter the stress history, and therefore the sti ﬀ ness and strength of the soil at shallow depth, with important consequence to the lateral behavior of piles. The existing study is mainly focused on small-diameter piles under scouring, where the soil around a pile is analyzed with two simpliﬁed approaches: (I) simply removing the scour layers without changing the strength and sti ﬀ ness of the remaining soils, or (II) solely considering the e ﬀ ects of stress history on the soil strength. This study aims to investigate and quantify the scour e ﬀ ect on the lateral behavior of monopile, based on an advanced hypoplastic model considering the inﬂuence of stress history on both soil sti ﬀ ness and strength. It is revealed that ignorance about the stress history e ﬀ ect (due to scouring) underestimates the extent of the soil failure wedge around the monopile, while overestimates soil sti ﬀ ness and strength. As a result, a large-diameter pile (diameter D = 5 m) in soft clay subjected to a souring depth of 0.5 D has experienced reductions in ultimate soil resistance and initial sti ﬀ ness of the p-y curves by 40% and 26%, and thus an increase of pile head deﬂection by 49%. Due to the inadequacy to consider the stress history e ﬀ ects revealed above, the existing approach (I) has led to non-conservative estimation, while the approach (II) has resulted in an over-conservative prediction.


Introduction
Scour is a process of soil erosion and can often occur around the foundations of offshore structures [1][2][3][4][5]. Currently, monopiles are the most widely employed foundation for offshore wind turbines, and its slenderness ratio of embedded pile length to pile diameter (L/D) are relatively small (typically in the range of [4][5][6][7][8]. Scour reduces the pile embedded length and changes the stress history of the remaining soils, which significantly influences pile responses and the natural frequency of wind turbines [5]. Thus, scour should be well considered during the design of wind turbines. Although extensive research efforts have been paid on sour effect on pile lateral responses [6-9], most of these studies have been largely limited to small-diameter piles, with ignorance of the stress history effect under scour conditions. The response of the laterally loaded pile under scouring is usually analyzed by two simplified approaches (I): simply removing the scour layers without changing the strength and stiffness of the remaining soils [7,8], or (II) solely considering the effects of stress history on the soil strength [10][11][12]. The approach (I) ignores the stress history effect due to scouring, which overestimates the undrained shear strength of the remaining soils (as shown in Figure 1), and Figure 1), and thus leads to a non-conservative estimation of monopile response. The influence of altering the stress history by scouring on the stiffness of soil, which governs the natural frequency and fatigue of monopile supported wind turbine, has not been taken into account in both approaches (I) and (II). At present, the p-y method is widely used for the analysis of laterally loaded piles. In this method, the pile is considered as an elastic beam and the soil is represented by a series of discrete py springs. The most widely used p-y curves were proposed by Matlock [13], which had been adopted in design codes. The Matlock p-y curves are formulated as follows: where p is lateral soil resistance per unit length of a pile; y is lateral pile deflection; yc is the lateral displacement at half the maximum soil stress, which can be determined by: where 50 is the strain at one-half the maximum stress. pu is the ultimate soil resistance per length, which is equal to the smaller value of pu1 and pu2 calculated by: where ' γ is the effective unit weight; su is the soil undrained shear strength; J is a constant value; z is the depth below the post-scour mudline. Lin et al. [10] modified Matlock p-y curves to consider the stress history effect. This was achieved by modifying the ultimate soil resistance, pu. The modification of pu depends on the change of the undrained shear strength and the effective unit weight of the remaining soils after scour, as follows [14]: At present, the p-y method is widely used for the analysis of laterally loaded piles. In this method, the pile is considered as an elastic beam and the soil is represented by a series of discrete p-y springs. The most widely used p-y curves were proposed by Matlock [13], which had been adopted in design codes. The Matlock p-y curves are formulated as follows: y y c ) 1/3 for y ≤ 8y c p u for y > 8y c (1) where p is lateral soil resistance per unit length of a pile; y is lateral pile deflection; y c is the lateral displacement at half the maximum soil stress, which can be determined by: where ε 50 is the strain at one-half the maximum stress. p u is the ultimate soil resistance per length, which is equal to the smaller value of p u1 and p u2 calculated by: where γ is the effective unit weight; s u is the soil undrained shear strength; J is a constant value; z is the depth below the post-scour mudline. Lin et al. [10] modified Matlock p-y curves to consider the stress history effect. This was achieved by modifying the ultimate soil resistance, p u . The modification of p u depends on the change of the undrained shear strength and the effective unit weight of the remaining soils after scour, as follows [14]: γ sc = 1 + e int 1 + e int + C ur log[ where s int u and s sc u are the soil undrained shear strength before and after scour, respectively; σ int and σ sc are the vertical effective stress before and after scour, respectively; OCR is the overconsolidated ratio of soil; Λ is a parameter (approximately 0.8); γ int and γ int are the effective unit weight before and after scour, respectively; e int is soil void ratio before scour; S d is scour depth; C c and C ur denote the compression and swelling indexes obtained from the oedometer tests.
By substituting Equations (5) and (6) into Equations (3) and (4), the equations for p u considering the stress history effect can be rewritten as follows: Lin's method [10] offered novel insights into the stress history effects due to scouring. Subsequently, Zhang et al. [12] and Liang et al. [13] further developed Lin's [10] p-y method to consider not only the stress history effect, but also scour-hole dimension and vertical load effect, respectively. However, the above-mentioned methods for considering the stress history effect are all based on Matlock p-y curves. It is well known that Matlock p-y curves were developed from a full-scale field lateral loading test on long, and flexible pile with a diameter of 0.324 m and the ratio of pile embedded length to diameter of 39.5. Its validity for large diameter (typically 4-6 m) monopiles of offshore wind turbines has been questioned by many researchers [15][16][17][18]. Thus, DNVGL [19] recommends that any proposed design method for large diameter monopiles should be validated by other means, such as by finite element calculations [18]. In addition, due to the stress-dependent behavior of soil, scour-induced stress loss decreases the soil stiffness of the remaining soils, which directly affects the initial stiffness of p-y curves. However, those methods for considering the stress history effect are based on the Matlock p-y curves, which make use of a parabolic curve shape, and thus the initial tangent modulus is infinite. In other words, those methods for considering the stress history effect are unable to reflect the influence on the initial stiffness of p-y curves.
The objective of this paper is to present a three-dimensional finite element method to investigate the stress history effect on the response of a monopile supported offshore wind turbine under scour conditions in soft clay. In this study, a full-scale offshore wind turbine was chosen as the reference structure. An advanced hypoplastic clay model that considers dependency of soil stiffness and strength on stress history was adopted. The three-dimensional model was validated with a published centrifuge test, and then analyses were performed to examine the stress history effect on monopile lateral responses at two different scour depths under the ultimate limit state (ULS) condition. Suitability of the p-y method, proposed by Lin et al. [10] for considering stress history was also assessed.

Three-Dimensional Finite Element Model
A three-dimensional finite element model of a full-scale monopile-supported wind turbine in clay was developed using the software ABAQUS. The properties of the wind turbine and monopile were taken from Ma et al. [20] and Shirzadeh et al. [21]. A schematic diagram of the overall structure is shown in Figure 2a. The monopile has a diameter of 5 m, a wall thickness of 0.07 m and an embedded length of 30 m. The diameter of the tower at the top and bottom are 3.4 and 4.4 m, respectively. Figure 2b shows an isometric view of the finite element mesh and the boundary conditions adopted in this study. The lateral boundary of the finite element mesh was constrained by roller supports, while the bottom boundary was fixed against translation in all directions. The lateral boundary was 10 D (D is the pile diameter) from the center of the pile, which is sufficient to eliminate the boundary effect [22].
The soil and the whole wind turbine structure were modeled using Eight-node brick with pore pressure (C3D8P) elements and Eight-node brick (C3D8) elements, respectively. The monopile foundation was assumed to be linear elastic with typical properties of Young's modulus of E p = 210 GPa and Poisson's ratio of ν p = 0.3. An advanced hypoplastic clay model (to be presented in the following subsection) was used to represent the soil behavior. The interaction between the pile and the soil was simulated based on the Coulomb friction law. The frictional coefficient µ = 0.31 was adopted in this study based on the equation proposed by Randolph and Wroth [23]. The detachment between the pile and the clay was allowed [24].
The tower stiffness, which could have a potential impact on the simulation, was properly considered by adopting the geometry and material of the tower of a 3 MW offshore wind turbine [20,21]. The referred wind turbine was founded on a monopile with a diameter of 5 m, being identical to the pile diameter adopted in this numerical investigation. embedded length of 30 m. The diameter of the tower at the top and bottom are 3.4 and 4.4 m, respectively. Figure 2b shows an isometric view of the finite element mesh and the boundary conditions adopted in this study. The lateral boundary of the finite element mesh was constrained by roller supports, while the bottom boundary was fixed against translation in all directions. The lateral boundary was 10D (D is the pile diameter) from the center of the pile, which is sufficient to eliminate the boundary effect [22].
The soil and the whole wind turbine structure were modeled using Eight-node brick with pore pressure (C3D8P) elements and Eight-node brick (C3D8) elements, respectively. The monopile foundation was assumed to be linear elastic with typical properties of Young's modulus of Ep = 210 GPa and Poisson's ratio of νp = 0.3. An advanced hypoplastic clay model (to be presented in the following subsection) was used to represent the soil behavior. The interaction between the pile and the soil was simulated based on the Coulomb friction law. The frictional coefficient μ = 0.31 was adopted in this study based on the equation proposed by Randolph and Wroth [23]. The detachment between the pile and the clay was allowed [24].
The tower stiffness, which could have a potential impact on the simulation, was properly considered by adopting the geometry and material of the tower of a 3 MW offshore wind turbine [20,21]. The referred wind turbine was founded on a monopile with a diameter of 5 m, being identical to the pile diameter adopted in this numerical investigation.

An Advanced Hypoplastic Model for Clay Considering Stress History Effect
An advanced hypoplastic clay model proposed by Masin [25] with considering stress history effect and small-strain stiffness was selected in this study to represent the soil constitutive model. A general formulation of the hypoplastic model can be written as [26][27][28]: where  T and D represent the objective stress rate and the Euler stretching tensor, respectively.
The hypoelastic tensor L is

An Advanced Hypoplastic Model for Clay Considering Stress History Effect
An advanced hypoplastic clay model proposed by Masin [25] with considering stress history effect and small-strain stiffness was selected in this study to represent the soil constitutive model. A general formulation of the hypoplastic model can be written as [26][27][28]: where • T and D represent the objective stress rate and the Euler stretching tensor, respectively. The hypoelastic tensor L is where I is a fourth-order identity tensor; a, c 1 , and c 2 are the three parameters, and can be calculated by (12) where ϕ' c is a model parameter that denotes critical friction angle; ν is a model parameter that controls the soil shear stiffness at large strain; Parameter α can be calculated by where λ* and κ* are model parameters defining the slope of the isotropic virgin compression and unloading line in the ln(1 + e) versus ln (p ' ) plane, respectively. (e and p ' denote void ratio and mean effective stress, respectively). The N is a second-order constitutive tensor and can be expressed as where Y and m denote the degree of non-linearity and tensorial quantity, respectively, and have the following equations with factor F given by where where p r = 1 kPa is a reference stress; N is a model parameter defining the position of the isotropic virgin compression line in the ln(1 + e) versus ln (p ' ) plane. The model is formulated based on incremental equations, which is distinctively different from the conventional elastoplastic framework, i.e., decomposing strains into elastic and plastic parts. The stress rate • T varies nonlinearly with the strain rate D due to the nonlinear form given by the Euclidian norm D , and thus there is no need to define yield surface when predicting nonlinear behavior.
In summary, the basic hypoplastic clay model requires five parameters, i.e., ϕ' c , N, λ*, κ*, and ν. The parameters are equivalent to those defined in the modified Cam clay model.
The basic hypoplastic clay model can predict the monotonic behavior of clay in the medium to large strain range. In order to consider the small-strain stiffness and stress history effect of clay, the so-called intergranular strain concept [29] is combined in the enhancement of the basic model. The intergranular strain δ is used as a new tensorial state variable and the normalized magnitude of δ is: where R is a model parameter that denotes the size of the elastic range. The direction of the intergranular strain δ is: The general stress-strain relation can be re-written as: where u is a fourth-order tensor that represents stiffness, and can be calculated using the following equation: where χ is a model parameter that controls the rate of stiffness degradation; m rat is a model parameter that can be quantified by the ratio between initial small-strain stiffness upon a 90 • strain path reversal and the initial stiffness upon a 180 • strain reversal; m R represent the initial small-strain stiffness upon a 180 • strain path reversal. m R can be calibrated to fit the initial stiffness G 0 , which is formulated by Wroth and Houlsby [30], as follows: where A g and n g are model parameters that reflect the stress dependency of small-strain stiffness. The evolution equation for the δ is: where β r is a model parameter that controls the rate of stiffness degradation. In summary, the advanced hypoplastic clay model adopted in this study consists of 11 model parameters in total. Five out of the 11 parameters, i.e., ϕ' c , N, λ*, κ* and ν are for the basic model. The six other parameters, i.e., R, m rat , β r , χ, A g , and n g are for the intergranular concept.

Parameter Calibration and Model Validation
The basic model parameters of kaolin clay, i.e., ϕ ' c , N, λ*, and κ* were obtained from Powrie [31] and Al-Tabbaa [32]. The parameters R, m rat , β r , and χ were calibrated against data reported by Benz [33] on small-strain stiffness of kaolin clay, as shown in Figure 3a. In order to calibrate the remaining parameters ν, A g , and n g , an undrained cyclic triaxial test was carried out. The kaolin clay sample was consolidated under an isotropic confining stress of 200 kPa, followed by 100 cycles of undrained cyclic compression. More details about the triaxial test can be found in He [34]. The confining stresses in the afore-mentioned elemental tests (for calibrating model parameters) generally do not exceed 200 kPa, except one case of 300 kPa in Benz [33]'s tests. This range of effective confining stress (i.e., p' ≤ 200 kPa) is relevant to that considered in the numerical investigation reported herein, i.e., p' value of the soil (γ' = 8 kN/m 3 , K 0 = 0.625, where K 0 is the coefficient of lateral earth pressure) along the 30 m deep monopile fall within 180 kPa. It is worth noting that the Kaolin clay can differ a lot depending on the manufacture. The aforementioned databases for the calibration of the model parameters and the trixial test were all based on the same type of Kaolin, i.e., Speswhite Kaolin clay.
Results of the cyclic triaxial test and calibration are shown in Figure 3b for comparison. It can be found that the hypoplastic model can reasonably capture the soil behavior. All of the 11 model parameters used in this study are summarized in Table 1.  Figure 3b for comparison. It can be found that the hypoplastic model can reasonably capture the soil behavior. All of the 11 model parameters used in this study are summarized in Table 1.  Given the scarcity of published experimental results on large diameter rigid piles (e.g., D = 5 m, as simulated in this study), the hypoplastic clay model was verified against centrifuge test results on a semi-rigid pile in soft clay (Hong et al. [35]). The model pile has a diameter of D = 0.8 m in prototype. Its embedded length (L) and load eccentricity (h) are 13.2 and 2 m in prototype, respectively. The slenderness ratio (L/D) of the model pile is therefore 16.5. A three-dimensional finite element model was established to simulate the centrifuge model. Figure 4 shows the computed and measured monotonic lateral load-deflection relationship at the pile head. A good agreement can be found between the computed and measured results. Hence, it can be concluded that the three dimensional  Given the scarcity of published experimental results on large diameter rigid piles (e.g., D = 5 m, as simulated in this study), the hypoplastic clay model was verified against centrifuge test results on a semi-rigid pile in soft clay (Hong et al. [35]). The model pile has a diameter of D = 0.8 m in prototype. Its embedded length (L) and load eccentricity (h) are 13.2 and 2 m in prototype, respectively. The slenderness ratio (L/D) of the model pile is therefore 16.5. A three-dimensional finite element model was established to simulate the centrifuge model. Figure 4 shows the computed and measured monotonic lateral load-deflection relationship at the pile head. A good agreement can be found between the computed and measured results. Hence, it can be concluded that the three dimensional finite element model adopted in this study is capable of modelling the pile-soil interaction. Although the hypoplastic model has reasonably predicted the lateral behavior of the semi-rigid pile, it is desirable to perform tests on large diameter rigid piles in the future for further verifying the predictive capability of the model against such piles. finite element model adopted in this study is capable of modelling the pile-soil interaction. Although the hypoplastic model has reasonably predicted the lateral behavior of the semi-rigid pile, it is desirable to perform tests on large diameter rigid piles in the future for further verifying the predictive capability of the model against such piles.

Load Case
An ultimate design load case, i.e., the Extreme Turbulence Model (ETM) wind load at rated wind speed combined with the 50-year Extreme Wave Height (EWH) [36], was considered in this study. The environmental site conditions adopted in this study are summarized in Table 2 [20]. The wind load acting on the hub Fhub was estimated as [37,38]: where ρa is the density of air; AR is the rotor swept area; CT is the thrust coefficient, and U is the wind speed.
The wind load acting on the tower z tower F of height z was calculated as [36]: where z tower A is the wind pressure area on the tower of height z; Cs is shape coefficient which equals to 0.5 for the tubular tower; Vz denotes the average wind speed as a function of height z. The normal wind speed profile is given by the power law [36]:

Load Case
An ultimate design load case, i.e., the Extreme Turbulence Model (ETM) wind load at rated wind speed combined with the 50-year Extreme Wave Height (EWH) [36], was considered in this study. The environmental site conditions adopted in this study are summarized in Table 2 [20]. The wind load acting on the hub F hub was estimated as [37,38]: where ρ a is the density of air; A R is the rotor swept area; C T is the thrust coefficient, and U is the wind speed. The wind load acting on the tower F z tower of height z was calculated as [36]: where A z tower is the wind pressure area on the tower of height z; C s is shape coefficient which equals to 0.5 for the tubular tower; V z denotes the average wind speed as a function of height z. The normal wind speed profile is given by the power law [36]: where V hub is the wind speed at the height of the hub z hub ; α is the power law exponent, which is assumed to be 0.2. Wave forces on the structure were calculated using the Morison's equation based on the linear Airy wave theory [36]: .
x dz (30) where dF wave is the horizontal wave load on a vertical element dz of the monopile at level z; C M is the inertia coefficient; C D is the drag coefficient; ρ is the mass density of the sea water; D t is the diameter of each section; .
x and ..
x are the wave-induced velocity and acceleration in the horizontal direction. Since current force is relatively small compared to wind and wave force, thus loads due to current are not considered for analysis.

Numerical Modelling Procedure
The validated 3D model was then developed to investigate the stress history effect on the lateral response of the monopile under scour conditions. The clay was assumed to be normally consolidated before scour. For the purpose of the case studies considered in this paper, two scour depths similar to Lin et al. [10], i.e., S d = 0.2 D (1 m) and S d = 0.5 D (2.5 m) were examined. The detailed procedures are as follows: Procedures of modelling scour ignoring stress history: (1) Compared with the no scour model (as shown in Figure 5a), maintaining the pile tip depth constant, a soil condition after scour was developed first, as shown in Figure 5b. (2) Initial K 0 stress of the soil was generated by a spatial calculation method available in ABAQUS in a Geostatic step. In this step, an equivalent pressure that equals the vertical stress of the scour layer was then applied on the soil surface. Due to the equivalent pressure, the soil stress of the remaining soil after scour keep remained unchanged. Therefore, the soil shear strength and other soil properties of the remaining soil after scour were assumed to be the same with those before scour. This operation can model scour ignoring the stress history effect, which is similar to the method often used in practice, i.e., just simply removing the scour layer while keeping the soil properties of the remaining soil unchanged. (3) Wished-in-place pile installation was achieved by changing appropriate elements to a linear elastic material of the pile. Pile installation effect was not considered for a reasonable simplification. (4) The loads described in Section 2.4 were applied on the structure.
Procedures of modelling scour considering the stress history effect: (1) A model without scour was first developed, and then the initial K 0 soil stress was achieved.
(2) Defining a special step for forming scour. In this step, the scour layer was removed by adding keywords in the ABAQUS input file, i.e., Model change, remove, as shown in Figure 5c. This operation models the unloading process when scouring, and thus takes accout of the stress history effect.

Undrained Shear Strength after Scour
The distributions of the overconsolidation ratio (OCR) of the clay at two different scour depths are shown in Figure 6a. After scour, the normally consolidated clay becomes overconsolidated. The OCR increases with increasing scour depth and decreases with soil depth and gradually approaches a normally consolidated condition at a greater depth.
The properties and stresses of the soil elements at different depths after scour were extracted from the model to calculate the undrained shear strength. Figure 6b shows the undrained shear strength (su) of the remaining clay after scour. When the stress history effect is ignored, the soil vertical stress and the OCR of the remaining soil keep unchanged. Therefore, the undrained shear strength of the remaining soil is almost the same as that in the condition of no scour, and could be fitted by a linear line, i.e., su = 1.66 z. However, the undrained shear strength of the remaining soil which considers the stress history effect is found to be decreased when compared with that of ignoring the stress history. In this study, ᴧ = 0.78 in Equation (5) provides the best agreement with the computed results, as shown in Figure 6b.

Undrained Shear Strength after Scour
The distributions of the overconsolidation ratio (OCR) of the clay at two different scour depths are shown in Figure 6a. After scour, the normally consolidated clay becomes overconsolidated. The OCR increases with increasing scour depth and decreases with soil depth and gradually approaches a normally consolidated condition at a greater depth.
The properties and stresses of the soil elements at different depths after scour were extracted from the model to calculate the undrained shear strength. Figure 6b shows the undrained shear strength (s u ) of the remaining clay after scour. When the stress history effect is ignored, the soil vertical stress and the OCR of the remaining soil keep unchanged. Therefore, the undrained shear strength of the remaining soil is almost the same as that in the condition of no scour, and could be fitted by a linear line, i.e., s u = 1.66 z. However, the undrained shear strength of the remaining soil which considers the stress history effect is found to be decreased when compared with that of ignoring the stress history. In this study, Λ = 0.78 in Equation (5) provides the best agreement with the computed results, as shown in Figure 6b.   Table 3. It should be noted that, at any given scour depth, the percentage increases in pile head deflection presented in Table 3 are relative to the values of that with ignoring the stress history effects.

Lateral Load-Deflection Response
(a)   Table 3. It should be noted that, at any given scour depth, the percentage increases in pile head deflection presented in Table 3 are relative to the values of that with ignoring the stress history effects.   Table 3. It should be noted that, at any given scour depth, the percentage increases in pile head deflection presented in Table 3 are relative to the values of that with ignoring the stress history effects.

Lateral Load-Deflection Response
(a)  Under the scour conditions, considering the stress history effect results in 13% (Sd = 0.2D) and 49% (Sd = 0.5D) higher pile-head deflection compared with the case in which the stress history effect is ignored. The percentage increase in pile-head deflection between considering and ignoring the stress history effect is found to increase with increasing scour depth. Therefore, ignoring the stress history of the remaining soil is likely to cause an unconservative analysis of the laterally loaded pile under scour conditions. A similar conclusion is also made by Lin et al. [10] and Zhang et al. [12]. In addition, compared to the result of no scour, considering scour and the resulted stress history effect leads to a 16% (Sd = 0.2D) and 64% (Sd = 0.5D) increase in lateral pile-head deflection. It is recommended that scour and the accompanying stress history effect should be well treated when designing the monopile-supported wind turbines in clay.
For comparison, the calculated lateral load-deflection relationships at pile head by using the modified p-y curves proposed by Lin et al. [10] (see Equations (1), (2) and (5)-(8)) at scour depth of Sd = 0.2D and Sd = 0.5D are also included in the Figure 7a and Figure 7b, respectively. Since scour has an insignificant effect on the change of the effective unit weight of the remaining soil [10,12], thus it was ignored in this study.
As shown in the figures, at any given lateral pile-head displacement, the calculated force at the pile head based on Lin's [10] p-y method is lower than those computed. The difference is likely attributed to the factor that Lin's [10] p-y curves is developed based on Matlock's p-y method [11]. It has been well recognized that Matlock's p-y method [11] is mainly applicable to small-diameter flexible piles, but could significantly underestimate the lateral resistance of a large-diameter rigid pile, due to the ignorance of resistances from base shear force, base moment and skin friction [15][16][17][18]. On  Under the scour conditions, considering the stress history effect results in 13% (S d = 0.2 D) and 49% (S d = 0.5 D) higher pile-head deflection compared with the case in which the stress history effect is ignored. The percentage increase in pile-head deflection between considering and ignoring the stress history effect is found to increase with increasing scour depth. Therefore, ignoring the stress history of the remaining soil is likely to cause an unconservative analysis of the laterally loaded pile under scour conditions. A similar conclusion is also made by Lin et al. [10] and Zhang et al. [12]. In addition, compared to the result of no scour, considering scour and the resulted stress history effect leads to a 16% (S d = 0.2 D) and 64% (S d = 0.5 D) increase in lateral pile-head deflection. It is recommended that scour and the accompanying stress history effect should be well treated when designing the monopile-supported wind turbines in clay.
For comparison, the calculated lateral load-deflection relationships at pile head by using the modified p-y curves proposed by Lin et al. [10] (see Equations (1), (2) and (5)-(8)) at scour depth of S d = 0.2 D and S d = 0.5 D are also included in the Figures 7a and 7b, respectively. Since scour has an insignificant effect on the change of the effective unit weight of the remaining soil [10,12], thus it was ignored in this study.
As shown in the figures, at any given lateral pile-head displacement, the calculated force at the pile head based on Lin's [10] p-y method is lower than those computed. The difference is likely attributed to the factor that Lin's [10] p-y curves is developed based on Matlock's p-y method [11]. It has been well recognized that Matlock's p-y method [11] is mainly applicable to small-diameter flexible piles, but could significantly underestimate the lateral resistance of a large-diameter rigid pile, due to the ignorance of resistances from base shear force, base moment and skin friction [15][16][17][18]. On the other hand, all these factors have been implicitly considered in the 3D finite element analyses reported herein. Nevertheless, the percentage difference between ignoring and considering the stress history effect are comparable. As can be seen in Table 3, by using Lin's [10] p-y method, the calculated pile-head deflection when considering stress history effect increases by 25% (S d = 0.2 D) and 75% (S d = 0.5 D) compared with that of ignoring stress history effect. The differences are higher than those computed by 3D FE analysis, i.e., 13% for S d = 0.2 D and 49% for S d = 0.5 D. The comparison shows that Lin's [10] p-y method overestimates the percentage difference in pile-head deflection between ignoring and considering stress history effect.
It would be not possible to model lateral behavior of piles under unsymmetrical and irregularly shaped scour with p-y curves. Instead, finite element method, as adopted in this study, has an advantage to account for these effects. Figure 8 shows the bending moment profiles at different scour depths. Generally speaking, the computed maximum bending moments in the pile are slightly lower than that calculated by Lin's [10] p-y method. This difference may be due to the inherent difference between the 3D FE analysis and the p-y method. When considering the stress history effect, the location of maximum bending moment shifts toward to a greater depth and results in 1% (S d = 0.2 D) and 2% (S d = 0.5 D) higher maximum bending moment compared with the case in which stress history is neglected. It can also be found that the percentage difference in the maximum bending moment, between considering and ignoring stress history effect increases insignificantly with increasing scour depth. The results indicate that the stress history effect may have a minor influence on the maximum bending moment in the pile. Besides, when the scour depth increases to 0.2 D and 0.5 D, the maximum bending moment increases by approximately 2% and 5%, respectively, compared with that under no scour condition. To some extent, the scour and the stress history effect on the maximum bending moment in the pile can also be ignored. the other hand, all these factors have been implicitly considered in the 3D finite element analyses reported herein. Nevertheless, the percentage difference between ignoring and considering the stress history effect are comparable. As can be seen in Table 3, by using Lin's [10] p-y method, the calculated pile-head deflection when considering stress history effect increases by 25% (Sd = 0.2D) and 75% (Sd = 0.5D) compared with that of ignoring stress history effect. The differences are higher than those computed by 3D FE analysis, i.e., 13% for Sd = 0.2D and 49% for Sd = 0.5D. The comparison shows that Lin's [10] p-y method overestimates the percentage difference in pile-head deflection between ignoring and considering stress history effect. It would be not possible to model lateral behavior of piles under unsymmetrical and irregularly shaped scour with p-y curves. Instead, finite element method, as adopted in this study, has an advantage to account for these effects. Figure 8 shows the bending moment profiles at different scour depths. Generally speaking, the computed maximum bending moments in the pile are slightly lower than that calculated by Lin's [10] p-y method. This difference may be due to the inherent difference between the 3D FE analysis and the p-y method. When considering the stress history effect, the location of maximum bending moment shifts toward to a greater depth and results in 1% (Sd = 0.2D) and 2% (Sd = 0.5D) higher maximum bending moment compared with the case in which stress history is neglected. It can also be found that the percentage difference in the maximum bending moment, between considering and ignoring stress history effect increases insignificantly with increasing scour depth. The results indicate that the stress history effect may have a minor influence on the maximum bending moment in the pile. Besides, when the scour depth increases to 0.2D and 0.5D, the maximum bending moment increases by approximately 2% and 5%, respectively, compared with that under no scour condition. To some extent, the scour and the stress history effect on the maximum bending moment in the pile can also be ignored.

Soil Displacement Field
The computed soil displacement fields, as well as displacement vectors are shown in Figures 9  and 10 at the scour depth of S d = 0.2 D and S d = 0.5 D, respectively. Two distinct soil flow mechanisms can be clearly identified for the large diameter monopile, namely a wedge mechanism near the ground surface and rotational soil flow near the pile toe. Similar failure mechanisms were also observed by Hong et al. [35] and Schroeder et al. [39]. When considering the stress history effect, the width of the wedge failure zone on the ground surface extends from 1.6 D to 1.8 D (S d = 0.2 D) and 1.7 D to 2.5 D (S d = 0.5 D). Meanwhile, the wedge failure zone is observed to extend to a greater depth, i.e., from 0.53 L (L = pile embedded length before scour) to 0.57 L (S d = 0.2 D) and 5.7 L to 6.7 L (S d = 0.5 D). As expected, the differences in the width and depth of the wedge failure zone between considering and ignoring stress history effect increase with increasing scour depth. As for the rotation center of the plane rotation zone, when considering stress history, it moves downward from 0.76 L to 0.78 L (S d = 0.2 D) and 0.78 L to 0.82 L (S d = 0.5 D). As a conclusion, soil failure mechanism of the large diameter monopile consists of two parts, namely wedge failure at shallow and rational soil flow at depth. Ignoring the stress history effect underestimates the width and depth of the wedge failure zone, while overestimates the location of the rotational soil flow zone.

Soil Displacement Field
The computed soil displacement fields, as well as displacement vectors are shown in Figures 9  and 10 at the scour depth of Sd = 0.2D and Sd = 0.5D, respectively. Two distinct soil flow mechanisms can be clearly identified for the large diameter monopile, namely a wedge mechanism near the ground surface and rotational soil flow near the pile toe. Similar failure mechanisms were also observed by Hong et al. [35] and Schroeder et al. [39]. When considering the stress history effect, the width of the wedge failure zone on the ground surface extends from 1.6D to 1.8D (Sd = 0.2D) and 1.7D to 2.5D (Sd = 0.5D). Meanwhile, the wedge failure zone is observed to extend to a greater depth, i.e., from 0.53L (L = pile embedded length before scour) to 0.57L (Sd = 0.2D) and 5.7L to 6.7L (Sd = 0.5D). As expected, the differences in the width and depth of the wedge failure zone between considering and ignoring stress history effect increase with increasing scour depth. As for the rotation center of the plane rotation zone, when considering stress history, it moves downward from 0.76L to 0.78L (Sd = 0.2D) and 0.78L to 0.82L (Sd = 0.5D). As a conclusion, soil failure mechanism of the large diameter monopile consists of two parts, namely wedge failure at shallow and rational soil flow at depth. Ignoring the stress history effect underestimates the width and depth of the wedge failure zone, while overestimates the location of the rotational soil flow zone.

p-y Curves Derived from Finite Element Simulation Results
To investigate the stress history effect on p-y curves and to assess the suitability of the modified p-y curves proposed by Lin et al. [10] for considering the stress history effect. The pile-soil contact force and the corresponding lateral pile displacement were extracted from the numerical results to deduce the p-y curves. Figure 11a,b show the extracted typical p-y curves of the remaining soil after scour at the depth of 0.5 and 1.5 m (measured from the post-scour mudline), respectively. The p-y curves proposed by Lin et al. [10] are also presented in the figure for comparison. It can be seen that both computed and calculated p-y curves reflect a similar trend that is at any given lateral displacement, the p-y curves of considering the stress history effect have much lower lateral soil resistance than that of ignoring the stress history effect. The reduction in lateral soil resistance increases the pile-head deflection and the maximum bending moment in the pile.

p-y Curves Derived from Finite Element Simulation Results
To investigate the stress history effect on p-y curves and to assess the suitability of the modified p-y curves proposed by Lin et al. [10] for considering the stress history effect. The pile-soil contact force and the corresponding lateral pile displacement were extracted from the numerical results to deduce the p-y curves. Figure 11a,b show the extracted typical p-y curves of the remaining soil after scour at the depth of 0.5 and 1.5 m (measured from the post-scour mudline), respectively. The p-y curves proposed by Lin et al. [10] are also presented in the figure for comparison. It can be seen that both computed and calculated p-y curves reflect a similar trend that is at any given lateral displacement, the p-y curves of considering the stress history effect have much lower lateral soil resistance than that of ignoring the stress history effect. The reduction in lateral soil resistance increases the pile-head deflection and the maximum bending moment in the pile. It should be noted that at any given depth, the p-y curves proposed by Lin et al. [10] show much lower soil resistance than the computed p-y curves. This is because the p-y curves proposed by Lin et al. [10] are based on Matlock p-y curves which underestimate the ultimate lateral soil resistance. The underestimation has been identified by many researchers based on pile tests and numerical modelling [17,18,40,41]. It should also be noted that the soil resistance of the p-y curves proposed by Lin et al., [10] presented in the figure is far from reaching the ultimate soil resistance. In other words, the p-y curves proposed by Lin et al. [10] significantly overestimate the required lateral displacement to reach the ultimate soil resistance. Due to the much softer p-y curves, the calculated pile-head deflections are much larger than that computed by 3D FE analysis (see Table 3).
The limited number of three-dimensional analyses reported herein (see Figure 11 for example), which aim to illustrate the typical influences of changing stress history (by sour) on soil-pile interaction, are still not sufficient for rigorously formulating a modification of the p-y method. It will It should be noted that at any given depth, the p-y curves proposed by Lin et al. [10] show much lower soil resistance than the computed p-y curves. This is because the p-y curves proposed by Lin et al. [10] are based on Matlock p-y curves which underestimate the ultimate lateral soil resistance. The underestimation has been identified by many researchers based on pile tests and numerical modelling [17,18,40,41]. It should also be noted that the soil resistance of the p-y curves proposed by Lin et al., [10] presented in the figure is far from reaching the ultimate soil resistance. In other words, the p-y curves proposed by Lin et al. [10] significantly overestimate the required lateral displacement to reach the ultimate soil resistance. Due to the much softer p-y curves, the calculated pile-head deflections are much larger than that computed by 3D FE analysis (see Table 3).
The limited number of three-dimensional analyses reported herein (see Figure 11 for example), which aim to illustrate the typical influences of changing stress history (by sour) on soil-pile interaction, are still not sufficient for rigorously formulating a modification of the p-y method. It will be the authors' future pursuit to propose modified p-y method considering different pile diameters, sour depths and geometries, by performing several series of comprehensive numerical parametric study. Figure 12a,b present the distributions of the ultimate soil resistance along pile length at scour depth of 0.2 D and 0.5 D, respectively. The ultimate soil resistance calculated based on Lin's [10] p-y curves is also included in the figure. It should be pointed out that for 3D FE analysis, the soil resistance at deep depth may not be fully mobilized as the lateral displacement was quite small. Thus, those p-y curves were fitted by a hyperbolic function to obtain the ultimate soil resistance [42,43]. As expected, consideration of the stress history of the remaining soil results in a decrease in the ultimate soil resistance when compared with the results obtained when ignoring the stress history effect. A comparison between Figure 12a,b shows that, the difference in the ultimate soil resistance between considering and ignoring the stress history effect increases with increasing the scour depth. be the authors' future pursuit to propose modified p-y method considering different pile diameters, sour depths and geometries, by performing several series of comprehensive numerical parametric study. Figure 12a,b present the distributions of the ultimate soil resistance along pile length at scour depth of 0.2D and 0.5D, respectively. The ultimate soil resistance calculated based on Lin's [10] p-y curves is also included in the figure. It should be pointed out that for 3D FE analysis, the soil resistance at deep depth may not be fully mobilized as the lateral displacement was quite small. Thus, those p-y curves were fitted by a hyperbolic function to obtain the ultimate soil resistance [42,43]. As expected, consideration of the stress history of the remaining soil results in a decrease in the ultimate soil resistance when compared with the results obtained when ignoring the stress history effect. A comparison between Figure 12a,b shows that, the difference in the ultimate soil resistance between considering and ignoring the stress history effect increases with increasing the scour depth. To quantitatively evaluate the effect of stress history on ultimate soil resistance, the percentage reduction in ultimate soil resistance between considering stress history and ignoring stress history is plotted in Figure 13. At any given scour depth, the percentage reductions in ultimate soil resistance in Figure 13 are relative to the values of that with ignoring the stress history effects. As presented in the figure, the percentage reduction decreases when the soil depth increases. Based on the computed results, the percentage reduction in ultimate soil resistance near ground surface between considering stress history and ignoring stress history can be up to 30.1% (Sd = 0.2D) and 39.8% (Sd = 0.5D), which is lower than that when using the p-y curves proposed by Lin et al. [10], i.e., the percentage reduction is 36.2% (Sd = 0.2D) and 52.3% (Sd = 0.2D). The comparison demonstrates that Lin's [10] p-y method overestimates the percentage reduction in ultimate soil resistance between considering and ignoring the stress history effect. The overestimation leads to a larger percentage difference in pile-head deflection when compared with that obtained from 3D FE analysis (see Table 3). To quantitatively evaluate the effect of stress history on ultimate soil resistance, the percentage reduction in ultimate soil resistance between considering stress history and ignoring stress history is plotted in Figure 13. At any given scour depth, the percentage reductions in ultimate soil resistance in Figure 13 are relative to the values of that with ignoring the stress history effects. As presented in the figure, the percentage reduction decreases when the soil depth increases. Based on the computed results, the percentage reduction in ultimate soil resistance near ground surface between considering stress history and ignoring stress history can be up to 30.1% (S d = 0.2 D) and 39.8% (S d = 0.5 D), which is lower than that when using the p-y curves proposed by Lin et al. [10], i.e., the percentage reduction is 36.2% (S d = 0.2 D) and 52.3% (S d = 0.2 D). The comparison demonstrates that Lin's [10] p-y method overestimates the percentage reduction in ultimate soil resistance between considering and ignoring the stress history effect. The overestimation leads to a larger percentage difference in pile-head deflection when compared with that obtained from 3D FE analysis (see Table 3). Due to the stress dependency of soil behavior, the loss of the overburden stress caused by scour decreases the stiffness of the remaining soil, and finally results in a lower initial stiffness of the p-y curves. The p-y curves proposed by Lin et al. [10] are based on Matlock p-y curves, and thus adopt the same parabolic curve shape, making the initial stiffness of the p-y curves infinite. Therefore, the p-y curves proposed by Lin et al. [10] cannot reflect the stress history effect on the initial stiffness of the p-y curves which has a significant influence on natural frequency and fatigue life of wind turbine structure [44].

Depth (m)
Since the initial stiffness of the p-y curves proposed by Lin et al. [10] are infinite, only the distributions of the initial stiffness of the computed p-y curves at scour depth of Sd = 0.2D and Sd = 0.5D are presented in Figure 14a,b, respectively. In this study, the secant modulus at a small displacement, i.e., y = D/1000, of the p-y curves were regarded as initial stiffness [44]. It can be clearly seen in the figure that considering stress history effect leads to a decrease in the initial stiffness of the p-y curves when compared with the results that ignoring the stress history. The percentage reduction in initial stiffness between considering and ignoring stress history effect is further examined in Figure  15. At any given scour depth, the percentage reductions in the initial stiffness of p-y curves in Figure  15 are relative to the values of that with ignoring the stress history effects. The figure reveals that the percentage reduction gradually decreases when the soil depth increases. On the other hand, the percentage reduction in the initial stiffness of p-y curves is found to increase with increasing scour depth. When the scour depth is 0.2D and 0.5D, the percentage reduction in the initial stiffness of the p-y curves can be up to 20.7% and 25.8%, respectively. The initial stiffness of p-y curves can affect the natural frequency of an offshore wind turbine directly. It was founded by Wang et al. [45] that 10-50% decrease in the initial stiffness of p-y curves could lead to a 4.6-6.6% drop in the natural frequency of the wind turbine, which has a dramatic effect on the wind turbine fatigue life [46]. Thus, it is recommended that the initial stiffness of p-y curves should be well evaluated.
It is worth noting that the soil properties and thus sour geometries generally vary spatially even within a homogeneous layer [47,48]. It will be the authors' future pursuit to integrate numerical Due to the stress dependency of soil behavior, the loss of the overburden stress caused by scour decreases the stiffness of the remaining soil, and finally results in a lower initial stiffness of the p-y curves. The p-y curves proposed by Lin et al. [10] are based on Matlock p-y curves, and thus adopt the same parabolic curve shape, making the initial stiffness of the p-y curves infinite. Therefore, the p-y curves proposed by Lin et al. [10] cannot reflect the stress history effect on the initial stiffness of the p-y curves which has a significant influence on natural frequency and fatigue life of wind turbine structure [44].
Since the initial stiffness of the p-y curves proposed by Lin et al. [10] are infinite, only the distributions of the initial stiffness of the computed p-y curves at scour depth of S d = 0.2 D and S d = 0.5 D are presented in Figure 14a,b, respectively. In this study, the secant modulus at a small displacement, i.e., y = D/1000, of the p-y curves were regarded as initial stiffness [44]. It can be clearly seen in the figure that considering stress history effect leads to a decrease in the initial stiffness of the p-y curves when compared with the results that ignoring the stress history. The percentage reduction in initial stiffness between considering and ignoring stress history effect is further examined in Figure 15. At any given scour depth, the percentage reductions in the initial stiffness of p-y curves in Figure 15 are relative to the values of that with ignoring the stress history effects. The figure reveals that the percentage reduction gradually decreases when the soil depth increases. On the other hand, the percentage reduction in the initial stiffness of p-y curves is found to increase with increasing scour depth. When the scour depth is 0.2 D and 0.5 D, the percentage reduction in the initial stiffness of the p-y curves can be up to 20.7% and 25.8%, respectively. The initial stiffness of p-y curves can affect the natural frequency of an offshore wind turbine directly. It was founded by Wang et al. [45] that 10-50% decrease in the initial stiffness of p-y curves could lead to a 4.6-6.6% drop in the natural frequency of the wind turbine, which has a dramatic effect on the wind turbine fatigue life [46]. Thus, it is recommended that the initial stiffness of p-y curves should be well evaluated.
It is worth noting that the soil properties and thus sour geometries generally vary spatially even within a homogeneous layer [47,48]. It will be the authors' future pursuit to integrate numerical modelling and spatial variability into the analysis of lateral behavior of piles under irregularly shaped scour.   fourth-order identity tensor N second-order constitutive tensor Y, m degree of non-linearity and tentorial quantity, respectively N position of the isotropic virgin compression line in the ln(1 + e) versus ln (p ' ) plane λ*, κ* slope of the isotropic virgin compression and unloading line in the ln(1 + e) versus ln (p ' ) plane, respectively ϕ' c critical state friction angle ν parameter controlling the shear stiffness δ intergranular strain R size of the elastic rangê δ direction of the intergranular strain u fourth-order tensor m rat path-dependent parameter β r , χ strain-dependent parameters A g , n g stress-dependent parameters G 0 soil initial stiffness F hub wind load acting on the hub ρ a , ρ density of air and sea water, respectively A R rotor swept area C T , C s thrust and shape coefficient, respectively F tower wind load acting on tower A tower wind pressure area on the tower V z average wind speed V hub wind speed at the height of the hub α power law exponent ϕ soil porosity F wave wave load c m , c d inertia and drag coefficient, respectively . x, ..
x wave-induced velocity and acceleration K 0 coefficient of lateral earth pressure k initial stiffness of p-y curves