Nonlinear Predictive Framework of the Undrained Clay Slope Effect on the Initial Stiffness of p-y Curves of Laterally Loaded Piles by FEM

: The hyperbolic p-y curve method is commonly used to design laterally loaded piles, in which the initial stiffness is one of the two key parameters that need to be determined. In this paper, the effect of an undrained clay slope on the initial stiffness of p-y curves of laterally loaded piles was explored, and nonlinear models of the reduction factor ( µ ) were proposed. A series of ﬁnite-element analyses was performed for different pile–slope geometric relationships according to whether the slope geometry had inﬂuence on pile–soil interaction, in which the geometrical parameters were varied. Based on simulation results, a semi-empirical method was mainly used to derive nonlinear formulations for the undrained clay slope effect on the initial stiffness of hyperbolic p-y curves. The wedge failure theory was also used to analyze the effect of the dimensionless slope height on the initial stiffness. Comparing with other researchers’ reduction factor models, the results of the present model are in the reasonable range and can predict more cases. Test cases were used to validate the proposed model, and the results show that theoretical results are in better agreement with test results using the present model.


Introduction
Piles have been widely used in marine geotechnical engineering to support axial and lateral loads for different structures, such as offshore wind turbines, bridges, and offshore drilling platforms constructed at or near the slope [1][2][3][4]. There are several methods to analyze the load-displacement characteristics of laterally loaded piles. The p-y curves method is known to be a relatively accurate and effective method [5,6]. Many investigators have developed p-y curves for sand [7][8][9], clay [10][11][12], and rock masses, [13,14] mainly based on laboratory tests [15,16], field tests [17], and numerical simulation [18][19][20]. A popular hyperbolic p-y curves method is commonly used to predict the load-displacement characteristics of the soil-pile system under lateral load. There are two key parameters in p-y curves, which are initial stiffness (K i ) and ultimate resistance of the soil (p u ) [13,21,22]. The ultimate resistance of the soil (p u ) determines the ultimate state, and the initial stiffness (K i ) determines the initial state for a pile under lateral load. Accurately calculating p u and K i is a crucial step to establish p-y curves, and it is meaningful in the design of laterally loaded piles.
In clayey soils, the ultimate soil resistance (p u ) is found to increase nonlinearly with depth, and a consensus regarding how this term is evaluated is evident in the literature [18,[23][24][25]. However, different formulas of the initial stiffness (K i ) have been proposed based on the results of laboratory tests, in-situ tests, and finite-element methods. Vesic [26] proposed an elastic solution model which was deduced from the theory of a beam on an elastic Winkler foundation. Rajashree and Sitharam [5] suggested that the value of initial stiffness could be taken as twice that of Vesic's model to account for real soil resistance around the circumference of the pile, a suggestion which was widely used by researchers. Then the model was modified to take into account the effect of diameter or depth [13,[27][28][29]. These models considered the effect of different soil-pile parameters and resulted in distinct engineering performances [30]. However, these models just focus on the case of piles on level ground. When the pile is embedded in or near the slopes, the initial stiffness is decreased because of the weakening of the soil in front of the pile [22], which is rather different from that of level ground.
The slope effect on the pile-soil response often refers to the weakening effect of the slope parameters (e.g., slope angle, slope height, near-slope distance and so on) on the pile lateral load capacity. Available information concerning the clay slope effect on the initial stiffness of p-y curves is rather limited and mainly refers to linear analysis. Several investigations have considered the equivalent problem of initial stiffness of p-y curves for the case of piles in or near clay slope. Georgiadis and Georgiadis [18] performed a series of FEA simulations of piles at relatively gentle slopes, and provided a reduction factor (µ) which represented the linear weakening effect of the slope angle. Then Georgiadis and Georgiadis [22] provided a smooth transition from a pile at the slope crest to a pile near the slope crest. However, the equation proved to be relatively accurate only for slope angles lower than 30 • . Yang et al. [16,31] performed laboratory tests and FEA simulations for a pile at the slope crest or in the slope. The results indicated that the critical depth (the depth with the initial stiffness equal to one) increased with the slope angle. A linear model of the reduction factor (µ) was then provided. Jiang et al. [9] provided an elastic stress solution of the soil in a slope under the self-weight, which proved to be accurate only for a narrow range of the slope angle. Considering that the current researches for the undrained clay slope effect on the initial stiffness are insufficient and have some weaknesses, systematical and comprehensive studies should be performed to explore the real slope effect on the initial stiffness of the pile-soil system.
To analyze the undrained clay slope effect on the initial stiffness, this paper established a series of 3D undrained finite-element models. Different pile-slope position relationships were considered according to whether the slope geometry had influence on pile-soil interaction. Based on simulation results, a semi-empirical method was mainly used to derive nonlinear formulations for the slope effect on the initial stiffness of hyperbolic p-y curves. The wedge failure theory was used to analyze the effect of the dimensionless slope height. Then the model performance of the present study was compared with those of other researchers' models. Moreover, the pile head load-displacement curves and p-y curves deduced with the proposed model were validated with other test results. The new model proved to be more appropriate and reasonable to fit real cases and could predict more cases.

Basic Assumptions
Based on the length relationship between the pile length and the slope height, two cases are considered in this paper: the case of "no slope geometry effect" and the case of "existence of the slope geometry effect". The case of "no slope geometry effect" refers to the slope geometry having no influence on the shape of the affected area with load increase. The case of "existence of slope geometry effect" refers to the slope geometry having influence on the shape of the affected area with load increase. The detail of the concept will be discussed in Section 4.2.1. Based on these two cases, considering the distance between the pile and slope crest, pile-slope geometry relationships can be divided into five detailed cases, as in Table 1.

(c) Dimensionless Pile Headcrest Distance /
No slope geometry effect

Variables (d) Dimensionless Slope Height (e) Pile Location in the Slope
Existence of the slope geometry effect Three parameters are focused on in the "no slope geometry effect" case, which are slope angle ( ), dimensionless near-slope distance ( / ) and dimensionless pile headcrest distance (ℎ/ ), as seen in Table 1(a-c). Two parameters are focused on in the "existence of the slope geometry effect" case, which are dimensionless slope height (ℎ ) and pile location in the slope ( ), as seen in Table 1(d,e).
The basic assumptions in this paper are described as follows: (1) The slope is stable and flat without a sliding surface, and the pile head is free.
(2) The pile bends because of the pile head lateral load. Based on the discussions on p-y curves by Reese et al. [21] and Terzaghi [32], the soil resistance is assumed to be nonlinear with pile displacement. The popular hyperbolic p-y curve method is adopted to reflect the nonlinear response of laterally loaded piles as the following: where is the ultimate lateral load per unit length, and is the initial stiffness. (3) The ultimate soil resistance ( ) in front of and behind the pile varies nonlinearly with depth [18,[23][24][25]. The initial stiffness ( ) of undrained clayey soil in the level ground case is independent on depth [26,28,29]. (4) Based on the study of Georgiadis and Georgiadis [18], a reduction factor ( ) is adopted to reflect the undrained clay slope effect on the initial stiffness: where is the initial stiffness in the slope condition, and is the initial stiffness in the level ground condition. For the convenience of readers, this paper provides the description of some used symbols, as shown in Table 2  Three parameters are focused on in the "no slope geometry effect" case, which are slope angle ( ), dimensionless near-slope distance ( / ) and dimensionless pile headcrest distance (ℎ/ ), as seen in Table 1(a-c). Two parameters are focused on in the "existence of the slope geometry effect" case, which are dimensionless slope height (ℎ ) and pile location in the slope ( ), as seen in Table 1(d,e).
The basic assumptions in this paper are described as follows: (1) The slope is stable and flat without a sliding surface, and the pile head is free.
(2) The pile bends because of the pile head lateral load. Based on the discussions on p-y curves by Reese et al. [21] and Terzaghi [32], the soil resistance is assumed to be nonlinear with pile displacement. The popular hyperbolic p-y curve method is adopted to reflect the nonlinear response of laterally loaded piles as the following: where is the ultimate lateral load per unit length, and is the initial stiffness. (3) The ultimate soil resistance ( ) in front of and behind the pile varies nonlinearly with depth [18,[23][24][25]. The initial stiffness ( ) of undrained clayey soil in the level ground case is independent on depth [26,28,29]. (4) Based on the study of Georgiadis and Georgiadis [18], a reduction factor ( ) is adopted to reflect the undrained clay slope effect on the initial stiffness: where is the initial stiffness in the slope condition, and is the initial stiffness in the level ground condition. For the convenience of readers, this paper provides the description of some used symbols, as shown in Table 2  Three parameters are focused on in the "no slope geometry effect" case, which are slope angle ( ), dimensionless near-slope distance ( / ) and dimensionless pile headcrest distance (ℎ/ ), as seen in Table 1(a-c). Two parameters are focused on in the "existence of the slope geometry effect" case, which are dimensionless slope height (ℎ ) and pile location in the slope ( ), as seen in Table 1(d,e).
The basic assumptions in this paper are described as follows: (1) The slope is stable and flat without a sliding surface, and the pile head is free.
(2) The pile bends because of the pile head lateral load. Based on the discussions on p-y curves by Reese et al. [21] and Terzaghi [32], the soil resistance is assumed to be nonlinear with pile displacement. The popular hyperbolic p-y curve method is adopted to reflect the nonlinear response of laterally loaded piles as the following: where is the ultimate lateral load per unit length, and is the initial stiffness. (3) The ultimate soil resistance ( ) in front of and behind the pile varies nonlinearly with depth [18,[23][24][25]. The initial stiffness ( ) of undrained clayey soil in the level ground case is independent on depth [26,28,29]. (4) Based on the study of Georgiadis and Georgiadis [18], a reduction factor ( ) is adopted to reflect the undrained clay slope effect on the initial stiffness: where is the initial stiffness in the slope condition, and is the initial stiffness in the level ground condition. For the convenience of readers, this paper provides the description of some used symbols, as shown in Table 2  Three parameters are focused on in the "no slope geometry effect" case, which are slope angle ( ), dimensionless near-slope distance ( / ) and dimensionless pile headcrest distance (ℎ/ ), as seen in Table 1(a-c). Two parameters are focused on in the "existence of the slope geometry effect" case, which are dimensionless slope height (ℎ ) and pile location in the slope ( ), as seen in Table 1(d,e).
The basic assumptions in this paper are described as follows: (1) The slope is stable and flat without a sliding surface, and the pile head is free.
(2) The pile bends because of the pile head lateral load. Based on the discussions on p-y curves by Reese et al. [21] and Terzaghi [32], the soil resistance is assumed to be nonlinear with pile displacement. The popular hyperbolic p-y curve method is adopted to reflect the nonlinear response of laterally loaded piles as the following: where is the ultimate lateral load per unit length, and is the initial stiffness. (3) The ultimate soil resistance ( ) in front of and behind the pile varies nonlinearly with depth [18,[23][24][25]. The initial stiffness ( ) of undrained clayey soil in the level ground case is independent on depth [26,28,29]. (4) Based on the study of Georgiadis and Georgiadis [18], a reduction factor ( ) is adopted to reflect the undrained clay slope effect on the initial stiffness: where is the initial stiffness in the slope condition, and is the initial stiffness in the level ground condition. For the convenience of readers, this paper provides the description of some used symbols, as shown in Table 2. (5) Poulos' relative stiffness ( ) method is used to judge whether a pile is flexible or rigid, as follows [33]:

Variables (d) Dimensionless Slope Height (e) Pile Location in the Slope
Existence of the slope geometry effect Three parameters are focused on in the "no slope geometry effect" case, which are slope angle ( ), dimensionless near-slope distance ( / ) and dimensionless pile headcrest distance (ℎ/ ), as seen in Table 1(a-c). Two parameters are focused on in the "existence of the slope geometry effect" case, which are dimensionless slope height (ℎ ) and pile location in the slope ( ), as seen in Table 1(d,e).
The basic assumptions in this paper are described as follows: (1) The slope is stable and flat without a sliding surface, and the pile head is free.
(2) The pile bends because of the pile head lateral load. Based on the discussions on p-y curves by Reese et al. [21] and Terzaghi [32], the soil resistance is assumed to be nonlinear with pile displacement. The popular hyperbolic p-y curve method is adopted to reflect the nonlinear response of laterally loaded piles as the following: where is the ultimate lateral load per unit length, and is the initial stiffness. (3) The ultimate soil resistance ( ) in front of and behind the pile varies nonlinearly with depth [18,[23][24][25]. The initial stiffness ( ) of undrained clayey soil in the level ground case is independent on depth [26,28,29]. (4) Based on the study of Georgiadis and Georgiadis [18], a reduction factor ( ) is adopted to reflect the undrained clay slope effect on the initial stiffness: where is the initial stiffness in the slope condition, and is the initial stiffness in the level ground condition. For the convenience of readers, this paper provides the description of some used symbols, as shown in Table 2. (5) Poulos' relative stiffness ( ) method is used to judge whether a pile is flexible or rigid, as follows [33]: Three parameters are focused on in the "no slope geometry effect" case, which are slope angle (θ), dimensionless near-slope distance (b/D) and dimensionless pile head-crest distance (h/L), as seen in Table 1(a-c). Two parameters are focused on in the "existence of the slope geometry effect" case, which are dimensionless slope height (h p ) and pile location in the slope (l p ), as seen in Table 1(d,e).
The basic assumptions in this paper are described as follows: (1) The slope is stable and flat without a sliding surface, and the pile head is free.
(2) The pile bends because of the pile head lateral load. Based on the discussions on p-y curves by Reese et al. [21] and Terzaghi [32], the soil resistance is assumed to be nonlinear with pile displacement. The popular hyperbolic p-y curve method is adopted to reflect the nonlinear response of laterally loaded piles as the following: where p u is the ultimate lateral load per unit length, and K i is the initial stiffness. (3) The ultimate soil resistance (p u ) in front of and behind the pile varies nonlinearly with depth [18,[23][24][25]. The initial stiffness (K i ) of undrained clayey soil in the level ground case is independent on depth [26,28,29]. (4) Based on the study of Georgiadis and Georgiadis [18], a reduction factor (µ) is adopted to reflect the undrained clay slope effect on the initial stiffness: where K iθ is the initial stiffness in the slope condition, and K i0 is the initial stiffness in the level ground condition. For the convenience of readers, this paper provides the description of some used symbols, as shown in Table 2.
(5) Poulos' relative stiffness (K R ) method is used to judge whether a pile is flexible or rigid, as follows [33]: where E p I p is the bending stiffness, E s is the soil Young's modulus, and L is the pile length. When 0.0025 ≤ K R ≤ 0.208, the pile is not rigid or flexible, which is named elastic pile. Table 2. The description of some used symbols.

Symbols Description
Initial stiffness in the level ground condition K iθ Initial stiffness in the slope condition K ii Initial stiffness at the ground surface µ Reduction factor µ θ Reduction factor in the "no slope geometry effect" case with varying (θ) µ b/D Reduction factor in the "no slope geometry effect" case with varying (b/D) µ h/D Reduction factor in the "no slope geometry effect" case with varying (h/D) µ h p Reduction factor in the "existence of the slope geometry effect" case with varying h p µ l p Reduction factor in the "existence of the slope geometry effect" case with varying l p

Finite-Element Analysis
The finite-element analysis program ABAQUS was used to perform analyses in which the geometrical characteristics of the problem (slope angle θ, dimensionless near-slope distance b/D, dimensionless pile head-crest distance h/L, dimensionless slope height h p , pile location in slope l p , pile length L) were varied. Material properties were the same as Georgiadis and Georgiadis' numerical simulation data [18], as shown in Table 3, and variable values are shown in Table 4. The soil was modeled as a linear elastic-perfectly plastic Tresca material [34][35][36] with the undrained Young's modulus E s = 14 MPa, the undrained shear strength c u = 70 kPa, Poisson ratio ν u = 0.49, and the bulk unit weight γ s = 18 kN/m 3 . Since the loading was undrained, the limiting pile-soil adhesion τ = α s c u , where α s is adhesion factor. Considering the lower bound of the relationship between α s and c u , α s was calculated to be 0.5 for the contact between the pile and soil, according to Georgiadis and Georgiadis' study [18]. The breadth and the bottom depth of the FE model were set to be ten times and six times the pile diameter, respectively, which is enough to eliminate the boundary effects [37,38].
All piles were modeled with a diameter (D) of 1 m. The piles were assumed to behave in a linear elastic manner [18,36], with Young's modulus E p = 2.9 × 10 7 kPa, Poisson's ratio υ = 0.1, and the bulk unit weight γ p = 25 kN/m 3 . Piles were modeled with lengths of L = 4 m, 12 m, and 20 m. The relative stiffness (K R ) of the three length piles are calculated in Equation (3) to be 0.3972, 0.0050, and 0.0006, respectively. Thus, the three piles can be regarded as rigid pile, elastic pile, and flexible pile, respectively.
More than 80 models were established in the present study. A typical 3D finiteelement meshing of the 30 • slope soil is shown in Figure 1. Eight-node quadrilateral (C3D8R) interface elements were used for the majority part of both the pile and the soil. There was more detailed meshing around the pile to allow effective pile-soil separating at the back of the pile. In this method, the total element number of models was around 8000 to 20,000. Each analysis step was performed referring to the construction process of drilled piles. The simulation could be divided into five construction parts. First of all, the bottom boundary of the mesh was fixed in all directions, and other boundaries were fixed only in the normal direction. Secondly, the geostatic equilibrium of soil domain under self-weight was considered without the pile. Thirdly, the soil elements located at pile positions were removed and the pile elements were activated. The normal and tangential contact between pile and soil were applied. Fourthly, self-weight was applied on the pile through gravity loading, and the overall model was balanced again. Finally, the lateral load was applied on the top of the pile and the response of the pile was analyzed. Note that each load increment was 100 kN for 12 m and 20 m piles and 20 kN for 4 m piles. For the case of θ = 90 • , all boundary conditions would work at first to balance the crustal stress, and then the boundary near the pile-placed side would be removed to form a 90 • slope. It is noted that large angles like 75 • and 90 • are rare in engineering, and are just used to explore the real slope effect on the initial stiffness in this paper.  Figure 2 shows the comparison between FEA results and the theoretical method of Georgiadis and Georgiadis [18]. The pile length is 20 m, the slope angles are 15°, 45°, and 60°, and the other parameters are the same as shown in Table 1. The results prove that the FEA results are within the reasonable range. To derive p-y curves, the curves of shear force ( ) versus depth ( ) were exported from the nodes of the piles of FEA results. The Q-z curves were curve-fitted in a relatively accurate method of seven-order polynominal function. Then the polynominal functions were differentiated to give curves of lateral load per unit pile length ( ) and depth ( ) for each pile head lateral load. The obtained p-z curves  Figure 2 shows the comparison between FEA results and the theoretical method of Georgiadis and Georgiadis [18]. The pile length is 20 m, the slope angles are 15 • , 45 • , and 60 • , and the other parameters are the same as shown in Table 1. The results prove that the FEA results are within the reasonable range. To derive p-y curves, the curves of shear force (Q) versus depth (z) were exported from the nodes of the piles of FEA results. The Q-z curves were curve-fitted in a relatively accurate method of seven-order polynominal function. Then the polynominal functions were differentiated to give curves of lateral load per unit pile length (p) and depth (z) for each pile head lateral load. The obtained p-z curves were combined with y-z curves to produce p-y curves. Then, the initial slopes of these p-y curves could be obtained as the initial stiffness (K i ). Then, the reduction factors (µ) could be obtained using Equation (2). The whole procedure is shown in Appendix B. The total number of data on the initial stiffness was around 700.  Figure 2 shows the comparison between FEA results and the theoretical method of Georgiadis and Georgiadis [18]. The pile length is 20 m, the slope angles are 15°, 45°, and 60°, and the other parameters are the same as shown in Table 1. The results prove that the FEA results are within the reasonable range. To derive p-y curves, the curves of shear force ( ) versus depth ( ) were exported from the nodes of the piles of FEA results. The Q-z curves were curve-fitted in a relatively accurate method of seven-order polynominal function. Then the polynominal functions were differentiated to give curves of lateral load per unit pile length ( ) and depth ( ) for each pile head lateral load. The obtained p-z curves were combined with y-z curves to produce p-y curves. Then, the initial slopes of these p-y curves could be obtained as the initial stiffness ( ). Then, the reduction factors ( ) could be obtained using Equation (2). The whole procedure is shown in Appendix B. The total number of data on the initial stiffness was around 700. Note that the used polynomial method may have resulted in discontinuity at the ground surface and pile end. Thus, the data on fitting results at ground surface were invalid and not used in the present study. Moreover, the data of the nodes below the first turning point were not considered (e.g., as shown in Figure 3a, though the pile length is 20 m, only the data for the depths from 1 m to 7 m are used, because the pile first turns at around the depth of 8 m). Note that the used polynomial method may have resulted in discontinuity at the ground surface and pile end. Thus, the data on fitting results at ground surface were invalid and not used in the present study. Moreover, the data of the nodes below the first turning point were not considered (e.g., as shown in Figure 3a, though the pile length is 20 m, only the data for the depths from 1 m to 7 m are used, because the pile first turns at around the depth of 8 m).   Dep Depth z / D

Theoretical Analysis for FEA Results
Based on cases with different slope heights, the pile-slope geometry relationships can be divided into two main cases, as Table 1 shows: (1) the "no slope geometry effect" case and (2) the "existence of the slope geometry effect" case. The following sub-sections will discuss the slope effect on the initial stiffness based on these two cases.

The "No Slope Geometry Effect" Case
In this case, three parameters are focused on, which are slope angle ( ), dimensionless near-slope distance ( / ), and dimensionless pile head-crest distance (ℎ/ ), as seen in Table 1(a-c). Theoretical analyses about the effect of these parameters are discussed in the following sub-sections, based on FEA simulation results.

Case 1: Variation of Slope Angle ( )
In this case, the pile is installed at the crest, as seen in Table 1(a). The reduction factor ( ) is calculated as shown in Figure 3. Different from the linear results of Georgiadis and Georgiadis [18] and Yang et al. [30], the reduction factor ( ) increases nonlinearly with the depth in this paper. These nonlinear curves indicate that the value of the initial stiffness ( ) in slope conditions approaches with depth. The empirical exponential equation is provided as following:

Theoretical Analysis for FEA Results
Based on cases with different slope heights, the pile-slope geometry relationships can be divided into two main cases, as Table 1 shows: (1) the "no slope geometry effect" case and (2) the "existence of the slope geometry effect" case. The following sub-sections will discuss the slope effect on the initial stiffness based on these two cases.

The "No Slope Geometry Effect" Case
In this case, three parameters are focused on, which are slope angle (θ), dimensionless near-slope distance (b/D), and dimensionless pile head-crest distance (h/L), as seen in Table 1(a-c). Theoretical analyses about the effect of these parameters are discussed in the following sub-sections, based on FEA simulation results.

Case 1: Variation of Slope Angle (θ)
In this case, the pile is installed at the crest, as seen in Table 1(a). The reduction factor (µ θ ) is calculated as shown in Figure 3. Different from the linear results of Georgiadis and Georgiadis [18] and Yang et al. [30], the reduction factor (µ θ ) increases nonlinearly with the depth in this paper. These nonlinear curves indicate that the value of the initial stiffness (K iθ ) in slope conditions approaches K i0 with depth. The empirical exponential equation is provided as following: where α and β are non-dimensional coefficients, which are determined in Section 4.1.2. Equation (4) shows that the reduction factor (µ θ ) increases from cos a θ at horizontal ground (z = 0) and nonlinearly approaches 1 with the depth ( z → ∞ ).

Case 2: Variation of Dimensionless Near-Slope Distance (b/D)
In this case, the pile is installed at a distance (b) from the slope crest, as seen in Table 1(b). Considering the effect of the near-slope distance, an idealization is illustrated in Figure 4, where K ii is the value of the initial stiffness at depth z = 0 at horizontal ground. It provides a smooth translation from a pile at the slope crest to one near the slope. By replacing depth (z) in Equation (4) with (z + z 1 ), where z 1 = (b − D/2) tan θ, it becomes: Equation (5) can also be expressed as follows: In the ideal case of = 90° at the slope crest, there is no Winkler spring in the outside part of the pile, so = 0 and / = 0. When a pile is not installed at the slope crest, there is still soil constraining the outside of the pile. In this case, Equation (6) becomes: Equation (7) demonstrates that for the ideal case of = 90°, the slope has the same effect on at any depth, and just changes with the distance ( ). Then, the parameter in Equation (7) can be determined through a series of / effect simulations in the case of = 90°. Figure 5 shows results of the slope angle ( ) effect on / by deriving the results from FEA. Averaging the results to obtain the ratio of slope case to level ground case, / becomes 0.32, 0.46, 0.62, 0.74, 0.81, and 0.88, respectively. Then, substituting / values into Equation (7), is obtained as 0.39, 0.42, 0.39, 0.38, 0.37, and 0.38, respectively. Thus, = 0.4 is suggested to be reasonable for use in this study. Equation (5) can also be expressed as follows: In the ideal case of θ = 90 • at the slope crest, there is no Winkler spring in the outside part of the pile, so K iθ = 0 and µ b/D = 0. When a pile is not installed at the slope crest, there is still soil constraining the outside of the pile. In this case, Equation (6) becomes: Equation (7) demonstrates that for the ideal case of θ = 90 • , the slope has the same effect on K iθ at any depth, and K iθ just changes with the distance (b). Then, the parameter β in Equation (7) can be determined through a series of b/D effect simulations in the case of θ = 90 • . Figure 5 shows results of the slope angle (θ) effect on µ b/D by deriving the results from FEA. Averaging the results to obtain the ratio of slope case to level ground case, µ b/D becomes 0.32, 0.46, 0.62, 0.74, 0.81, and 0.88, respectively. Then, substituting µ b/D values into Equation (7), β is obtained as 0.39, 0.42, 0.39, 0.38, 0.37, and 0.38, respectively. Thus, β = 0.4 is suggested to be reasonable for use in this study. Since is determined, can be determined in Equation (6) by fitting curves in Figure 3. Figure 6 shows the distribution of for a total of 95 data from piles. Then, a recommended value of = 1.2 is able to represent the most cases precisely. Since β is determined, α can be determined in Equation (6) by fitting curves in Figure 3. Figure 6 shows the distribution of α for a total of 95 data from piles. Then, a recommended value of α = 1.2. is able to represent the most cases precisely. Since is determined, can be determined in Equation (6) Figure 6 shows the distribution of for a total of 95 data from piles. Then, a recommended value of = 1.2 is able to represent the most cases precisely. Finally, now that = 1.2 and = 0.4 have been obtained, Equation (6) can be provided as: Additional cases of both the dimensionless near-slope distance ( / ) and the slope angle ( ) effects are conducted using FEA, and the results of are plotted in Figure 7, together with curves of Equation (8), which shows good agreement between predictions and simulation results. Finally, now that α = 1.2 and β = 0.4 have been obtained, Equation (6) can be provided as: Additional cases of both the dimensionless near-slope distance (b/D) and the slope angle (θ) effects are conducted using FEA, and the results of µ are plotted in Figure 7, together with curves of Equation (8), which shows good agreement between predictions and simulation results. In this case, the pile is installed on the slope, at a distance ℎ away from the slope crest, as seen in Table 1(c). As shown in Figure 8, the lateral bearing capacity is almost independent of the pile head-crest distance (ℎ). This phenomenon indicates that the soil behind the pile hardly participates in the soil-pile response for undrained clay slopes, which is similar to the results of Sivapriya et al.'s laboratory test [39]. This can be explained by the theory of Coulomb earth pressure, which considers that the active earth pressure is dependent on the soil properties and the height of the structure. Thus, the case of the In this case, the pile is installed on the slope, at a distance h away from the slope crest, as seen in Table 1(c). As shown in Figure 8, the lateral bearing capacity is almost independent of the pile head-crest distance (h). This phenomenon indicates that the soil behind the pile hardly participates in the soil-pile response for undrained clay slopes, which is similar to the results of Sivapriya et al.'s laboratory test [39]. This can be explained by the theory of Coulomb earth pressure, which considers that the active earth pressure is dependent on the soil properties and the height of the structure. Thus, the case of the variation of dimensionless pile head-crest distance (h/L) is nearly the same as that of the pile at the slope crest. The slope effect on the initial stiffness in this case can be expressed as: In this case, the pile is installed on the slope, at a distance ℎ away from the slope crest, as seen in Table 1(c). As shown in Figure 8, the lateral bearing capacity is almost independent of the pile head-crest distance (ℎ). This phenomenon indicates that the soil behind the pile hardly participates in the soil-pile response for undrained clay slopes, which is similar to the results of Sivapriya et al.'s laboratory test [39]. This can be explained by the theory of Coulomb earth pressure, which considers that the active earth pressure is dependent on the soil properties and the height of the structure. Thus, the case of the variation of dimensionless pile head-crest distance (ℎ/ ) is nearly the same as that of the pile at the slope crest. The slope effect on the initial stiffness in this case can be expressed as:  Figure 8. Effect of the slope height on pile head displacement of piles. Figure 8. Effect of the slope height on pile head displacement of piles.

The "Existence of the Slope Geometry Effect" Case
In this case, two parameters are focused on, which are dimensionless slope height (h p ) and pile location in slope (l p ), as seen in Table 1(d,e). Theoretical analyses about these parameters are discussed in the following sub-sections, based on simulation results.

Case 4: Variation of Dimensionless Slope Height (h p )
In this case, the pile is installed at the crest and dimensionless slope height (h p ) varies, as seen in Table 1(d). Figure 9 shows that the lateral bearing capacity decreases with the increase of h p from 0 to 1 for slope angle θ = 45 • . It is obvious that h p hardly affects the soil-pile response when h p is beyond a certain height (e.g., the curves in Figure 9a when h p ≥ 0.30). This indicates that the dimensionless slope height (h p ) affects the soil-pile response in a nonlinear way.
The phenomenon can be explained by a wedge failure theory first proposed by Reese and Welch and Reese et al. [40,41]. The wedge failure model considers that the soil-pile response relates to the lowest value of two possible failure mechanisms: a wedge failure mechanism at shallow depths and a flow failure around the pile at greater depths. Moreover, the passive wedge develops along the depth as the load increases to form a larger wedge failure body, as shown in Figure 10, where φ is the angle between the sliding surface and the vertical plane. Under undrained conditions, the angle (φ) is suggested to be 45 • by Reese et al. [41].
In this case, the pile is installed at the crest and dimensionless slope height (ℎ ) varies, as seen in Table 1(d). Figure 9 shows that the lateral bearing capacity decreases with the increase of ℎ from 0 to 1 for slope angle = 45°. It is obvious that ℎ hardly affects the soil-pile response when ℎ is beyond a certain height (e.g., the curves in Figure 9a when ℎ ≥ 0.30). This indicates that the dimensionless slope height (ℎ ) affects the soil-pile response in a nonlinear way.   The phenomenon can be explained by a wedge failure theory first proposed by Reese and Welch and Reese et al. [40,41]. The wedge failure model considers that the soil-pile response relates to the lowest value of two possible failure mechanisms: a wedge failure mechanism at shallow depths and a flow failure around the pile at greater depths. Moreover, the passive wedge develops along the depth as the load increases to form a larger  Yang et al. analyzed the near-slope failure wedge for three "existence of the slope geometry effect" cases. The study considered that as the load increased, the passive wedge expanded, and the passive wedge shape changed from a triangle to a concave quadrilateral [31]. However, a theoretical largest passive wedge exists because the depth of the bottom of the wedge will not exceed the first bending point of the pile. Similarly, for the changing ℎ case in this paper, the shape of the largest passive failure wedge changes from a concave quadrilateral for a small ℎ value to a triangle for a larger ℎ value, as shown in Figure 11, where is the depth of the first turning point, and is the depth Yang et al. analyzed the near-slope failure wedge for three "existence of the slope geometry effect" cases. The study considered that as the load increased, the passive wedge expanded, and the passive wedge shape changed from a triangle to a concave quadrilateral [31]. However, a theoretical largest passive wedge exists because the depth of the bottom of the wedge will not exceed the first bending point of the pile. Similarly, for the changing h p case in this paper, the shape of the largest passive failure wedge changes from a concave quadrilateral for a small h p value to a triangle for a larger h p value, as shown in Figure 11, where z t is the depth of the first turning point, and z cr is the depth of the intersection of the bottom plane of the largest passive wedge and the slope surface. Yang et al. analyzed the near-slope failure wedge for three "existence of the slope geometry effect" cases. The study considered that as the load increased, the passive wedge expanded, and the passive wedge shape changed from a triangle to a concave quadrilateral [31]. However, a theoretical largest passive wedge exists because the depth of the bottom of the wedge will not exceed the first bending point of the pile. Similarly, for the changing ℎ case in this paper, the shape of the largest passive failure wedge changes from a concave quadrilateral for a small ℎ value to a triangle for a larger ℎ value, as shown in Figure 11, where is the depth of the first turning point, and is the depth of the intersection of the bottom plane of the largest passive wedge and the slope surface. For "existence of the slope geometry effect" cases, when the slope angle remains constant, the shape of the theoretical largest passive wedge changes with ℎ when ℎ < / , but remains constant when ℎ ≥ / , as shown in Figure 11. This indicates that "existence of the slope geometry effect" cases are the same as "no slope geometry effect" cases when ℎ ≥ / . The reduction factor ( ) can be calculated for the dimensionless slope height (ℎ ) variation case as follows: For "existence of the slope geometry effect" cases, when the slope angle remains constant, the shape of the theoretical largest passive wedge changes with h p when h p < z cr /L, but remains constant when h p ≥ z cr /L, as shown in Figure 11. This indicates that "existence of the slope geometry effect" cases are the same as "no slope geometry effect" cases when h p ≥ z cr /L. The reduction factor (µ h p ) can be calculated for the dimensionless slope height (h p ) variation case as follows: where z cr is related to the slope angle, as shown in Figure 12, and can be calculated in geometry as follows: where is related to the slope angle, as shown in Figure 12, and can be calculated in geometry as follows: Calculation Method of Actually, the depth of the first turning point ( ) of the pile is dependent on soil and pile properties, loading conditions, and so on (e.g., for a certain pile, the value of approaches 0 when the soil stiffness is extremely large, but approaches the value of when Calculation Method of z t Actually, the depth of the first turning point (z t ) of the pile is dependent on soil and pile properties, loading conditions, and so on (e.g., for a certain pile, the value of z t approaches 0 when the soil stiffness is extremely large, but approaches the value of L when the soil stiffness is extremely small). Sun et al. [36] also indicated this, providing a range of 0.62 ∼ 0.73L for z t of a large-diameter monopile on sandy level ground. In this paper, the variation of z t is considered to be related to the relative stiffness (K R ). It is obvious that the relative stiffness (K R ) is influenced by the pile length because of the fourth order of the pile length (L) in Equation (3). Recording the depths of the first turning point (z t ) of piles derived from FEA results, as shown in Figure 13, the ratio (z t /L) results of the depth of the first turning point and the pile length are shown in Figure 14. The figure indicates that z t /L decreases as the pile length increases. For each length of the undrained, laterally loaded pile, z t /L remains almost the same value. However, when the pile is flexible, the lower part of the pile hardly bears the load. As a result, does not increase with the pile length. Thus, the value of / is around 0.7 for 4 m and 12 m length piles, but is as small as 0.4 for 20 m length piles, as shown in Figure 14. To solve the problem, a minimum flexible pile length ( ) is proposed in this paper, which can be calculated with Equation (3)   However, when the pile is flexible, the lower part of the pile hardly bears the load. As a result, does not increase with the pile length. Thus, the value of / is around 0.7 for 4 m and 12 m length piles, but is as small as 0.4 for 20 m length piles, as shown in Figure 14. To solve the problem, a minimum flexible pile length ( ) is proposed in this paper, which can be calculated with Equation (3)  However, when the pile is flexible, the lower part of the pile hardly bears the load. As a result, z t does not increase with the pile length. Thus, the value of z t /L is around 0.7 for 4 m and 12 m length piles, but is as small as 0.4 for 20 m length piles, as shown in Figure 14. To solve the problem, a minimum flexible pile length (L flex ) is proposed in this paper, which can be calculated with Equation (3) when K R = 0.0025, as follows: The minimum flexible pile length (L flex ) for the pile and soil in this paper is calculated to be 14.28 m (0.71L). This means that when the pile is longer than 14.28 m, the lower part of the pile will not bear the load, and z t will not increase with depth. Thus, in this case, L flex and z t are considered to be unchanged for the flexible pile. Then, for the flexible pile with a length of 20 m in this paper, z t = 0.4L(0.56L flex ) is finally obtained. The variation of z t can be regarded as 0.8L ∼ 0.6L flex for a pile varied from rigidity to flexibility, as shown in Figure 15. A simple calculation method of is provided in this paper, as seen in Equation (13). Further discussions of can be studied by researchers in the future. Elastic pile 0. 6 Flexible pile (13) Calculation Method of ℎ , As shown in Figure 11, the shape of the passive wedge is a triangle and does not change when ℎ ≥ / . At this time, the reduction factor ( ) remains the same as . Furthermore, the soil-pile system is not affected by the slope for cases with level ground (ℎ = 0) and = 1. Thus, a conclusion can be drawn that decreases from 1 to as ℎ increases from 0 to . To figure out the problem of how varies, a normalized reduction factor is proposed in this paper that can be calculated as follows: A simple calculation method of z t is provided in this paper, as seen in Equation (13). Further discussions of z t can be studied by researchers in the future.
Rigid pile 0.7L Elastic pile 0.6L flex Flexible pile (13) Calculation Method of f h p , µ θ As shown in Figure 11, the shape of the passive wedge is a triangle and does not change when h p ≥ z cr /L. At this time, the reduction factor (µ h p ) remains the same as µ θ . Furthermore, the soil-pile system is not affected by the slope for cases with level ground (h p = 0) and µ h p = 1. Thus, a conclusion can be drawn that µ h p decreases from 1 to µ θ as h p increases from 0 to z cr . To figure out the problem of how µ h p varies, a normalized reduction factor λ norm is proposed in this paper that can be calculated as follows: Equation (14) uses the normalized method to reflect the weight of µ h p for different values of µ θ . Figure 16 shows the nonlinear variation of λ norm with h p for different depths of piles. At deep depths, though curves show differences to shallow-depth curves, these curves are not considered because the µ h p of these curves is close to 1. Thus, the relationship between λ norm and h p is provided in a nonlinear form as follows: Rearranging Equation (15), can be calculated as the following when ℎ < / : where is calculated to be 1.6 m, 4.2 m, and 4.26 m with Equations (11) and (13), for pile lengths of 4 m, 12 m, and 20 m, respectively; is a parameter that controls the shape of the curve, as shown in Figure 16. Equation (16) reflects more characteristics of shallow depths with the increase of . To fit more cases of depths, = 2 is considered to be suitable and reasonable. Additionally, the numbers "1" and " " in Equation (16) can be re- Rearranging Equation (15), µ h p can be calculated as the following when h p < z cr /L: where z cr is calculated to be 1.6 m, 4.2 m, and 4.26 m with Equations (11) and (13), for pile lengths of 4 m, 12 m, and 20 m, respectively; k is a parameter that controls the shape of the curve, as shown in Figure 16. Equation (16) reflects more characteristics of shallow depths with the increase of k. To fit more cases of depths, k = 2 is considered to be suitable and reasonable. Additionally, the numbers "1" and "µ θ " in Equation (16) can be regarded as the reduction factor of the slope bottom (level ground) and that of the slope, respectively. Finally, Equation (10) can be expressed as: Figure 17 shows the curves derived using Equation (17) and the data from FEA, which show good agreement between the prediction and the simulation results.  Figure 17 shows the curves derived using Equation (17) and the data from FEA, which show good agreement between the prediction and the simulation results.

Case 5: Variation of Pile Location on Slope ( )
In this case, the pile is installed on the slope, at a distance ( ) away from the crest, as shown in Table 1(e), where = ℎ/ . Figure 18 shows that the lateral bearing capacity increases with , increasing from zero to one for cases where the slope angle = 45°. It is obvious that, although increases in a linear way, the lateral bearing capacity increases nonlinearly. Considering that the soil-pile response is independent of the soil behind the pile, as Case 3 shows, the case of variation of pile location on slope ( ) (case 5) can be regarded as the combination of Case 3 and Case 4. Replacing ℎ in Equation (17) with ( / − ), the reduction factor ( ) can be calculated using the following equation: In this case, the pile is installed on the slope, at a distance (l p ) away from the crest, as shown in Table 1(e), where l p = h/L. Figure 18 shows that the lateral bearing capacity increases with l p , increasing from zero to one for cases where the slope angle θ = 45 • . It is obvious that, although l p increases in a linear way, the lateral bearing capacity increases nonlinearly. Considering that the soil-pile response is independent of the soil behind the pile, as Case 3 shows, the case of variation of pile location on slope (l p ) (case 5) can be regarded as the combination of Case 3 and Case 4. Replacing h p in Equation (17) with H/L − l p , the reduction factor (µ l p ) can be calculated using the following equation:

Comparison with Other Methods
Since Georgiadis and Georgiadis [18] proposed the concept of the reduction factor ( ) of the initial stiffness for soil-pile response in undrained clay slopes, some researchers have developed the method. Figures 19 and 20 show the comparison of fromthe present study with the methods of Georgiadis and Georgiadis [18,22] and Yang et al. [16]. Figure  19a shows that approaches one with depth for both the present study and the method of Yang et al., but equals one below = 6 for the method of Georgiadis and Georgiadis.
Moreover, the present study shows a relatively smaller value of , which is more conservative for the engineering design. Figure 19b shows that the value of in the present study is similar to the value of the method of Georgiadis and Georgiadis, and the value of the method of Yang et al. is larger.

Comparison with Other Methods
Since Georgiadis and Georgiadis [18] proposed the concept of the reduction factor (µ) of the initial stiffness for soil-pile response in undrained clay slopes, some researchers have developed the method. Figures 19 and 20 show the comparison of µ fromthe present study with the methods of Georgiadis and Georgiadis [18,22] and Yang et al. [16]. Figure 19a shows that µ approaches one with depth for both the present study and the method of Yang et al., but equals one below z = 6D for the method of Georgiadis and Georgiadis. Moreover, the present study shows a relatively smaller value of µ, which is more conservative for the engineering design. Figure 19b shows that the value of µ in the present study is similar to the µ value of the method of Georgiadis Figure 20a shows that all the values of the three methods are nearly the same for small angle slopes. As the slope angle increases, the trend of in the present study is similar to Yang et al.'s method. However, Georgiadis and Georgiadis' method only shows good agreement with real cases for small slope angles, which is also indicated in Georgiadis and Georgiadis' study [18]. Figure 20b shows the similar trend of for both the present study and Yang et al.'s method, and the value of Yang et al.'s method is smaller for large slope angles. In the same way, the value of Georgiadis and Georgiadis' method only shows good agreement for small slope angles, and the value of the present study is better predicted for large slope angles. Additionally, note that both Georgiadis and Georgiadis' method and Yang et al.'s method cannot predict the ideal case of the slope angle = 90°, but it can be obtained by using the method of the present study.

Application to Test Cases
The present study is validated by two pile tests in this section. Usually, the soil-pile response is reflected by monitoring the pile head load ( ) and displacement ( ), which is known to be a − curve. The hyperbolic p-y curves in Equation (1) and the deflection differential equation are used to derive − curves as follows:   Figure 20a shows that all the values of the three methods are nearly the same for small angle slopes. As the slope angle increases, the trend of in the present study is similar to Yang et al.'s method. However, Georgiadis and Georgiadis' method only shows good agreement with real cases for small slope angles, which is also indicated in Georgiadis and Georgiadis' study [18]. Figure 20b shows the similar trend of for both the present study and Yang et al.'s method, and the value of Yang et al.'s method is smaller for large slope angles. In the same way, the value of Georgiadis and Georgiadis' method only shows good agreement for small slope angles, and the value of the present study is better predicted for large slope angles. Additionally, note that both Georgiadis and Georgiadis' method and Yang et al.'s method cannot predict the ideal case of the slope angle = 90°, but it can be obtained by using the method of the present study.

Application to Test Cases
The present study is validated by two pile tests in this section. Usually, the soil-pile response is reflected by monitoring the pile head load ( ) and displacement ( ), which is known to be a − curve. The hyperbolic p-y curves in Equation (1) and the deflection differential equation are used to derive − curves as follows:  Figure 20a shows that all the µ values of the three methods are nearly the same for small angle slopes. As the slope angle increases, the trend of µ in the present study is similar to Yang et al.'s method. However, Georgiadis and Georgiadis' method only shows good agreement with real cases for small slope angles, which is also indicated in Georgiadis and Georgiadis' study [18]. Figure 20b shows the similar trend of µ for both the present study and Yang et al.'s method, and the µ value of Yang et al.'s method is smaller for large slope angles. In the same way, the µ value of Georgiadis and Georgiadis' method only shows good agreement for small slope angles, and the µ value of the present study is better predicted for large slope angles. Additionally, note that both Georgiadis and Georgiadis' method and Yang et al.'s method cannot predict the ideal case of the slope angle θ = 90 • , but it can be obtained by using the method of the present study.

Application to Test Cases
The present study is validated by two pile tests in this section. Usually, the soil-pile response is reflected by monitoring the pile head load (H 0 ) and displacement (y 0 ), which is known to be a H 0 − y 0 curve. The hyperbolic p-y curves in Equation (1) and the deflection differential equation are used to derive H 0 − y 0 curves as follows: The finite difference method is used to transform differential equations into difference equations. The displacement of the laterally loaded pile can be obtained by solving these difference equations. Then the pile head load-displacement curve (H 0 − y 0 ) can also be obtained.
The values of the parameters used to calculate H 0 − y 0 curves are as shown in Table 5, where e is the height of the load applied above the ground level; c u is undrained shear strength; and E s is Young's modulus of the soil.  [42], and Georgiadis and Georgiadis [18]. The case is the same as the case in Section 4.1.1 where θ = 20 • , and the parameter values are shown in Table 5. The calculation method of p u and K i is the same as the method of Georgiadis and Georgiadis [18], as seen in Appendix A, and the method of the present study is used to calculate µ θ = K iθ /K i . The H 0 − y 0 curves are shown in Figure 21, which shows better agreement with the pile test when using the method of µ in the present study.
The finite difference method is used to transform differential equations into difference equations. The displacement of the laterally loaded pile can be obtained by solving these difference equations. Then the pile head load-displacement curve ( − ) can also be obtained.
The values of the parameters used to calculate − curves are as shown in Table  5, where is the height of the load applied above the ground level; is undrained shear strength; and is Young's modulus of the soil.  [42], and Georgiadis and Georgiadis [18]. The case is the same as the case in Section 4.1.1 where = 20°, and the parameter values are shown in Table 5. The calculation method of and is the same as the method of Georgiadis and Georgiadis [18], as seen in Appendix A, and the method of the present study is used to calculate = / . The − curves are shown in Figure 21, which shows better agreement with the pile test when using the method of in the present study.   [16] was used to validate the present study. The case is the same as the case in Section 4.2.1 where h p /L = 0.5, and the parameter values are shown in Table 5. The calculation method of p u and K i is the same as the method of Yang et al. [16], as seen in Appendix A, and the method of the present study is used to calculate µ h p = K ih p /K i . Note that z t is calculated to be 0.19 m in the case of the pile test. The H 0 − y 0 curves and p-y curves at z = 0 m depth are shown in Figure 22. Compared with Yang et al.'s method, the result of the present study proves to have good agreement with the pile test for the large load when the slope angle is small. Moreover, when the slope angle is large, the results of the present study are similar to the result of Yang et al.'s method. Additionally, the lateral load per unit (p) is calculated to be larger than the results of Yang et al.

Conclusions
The paper proposes nonlinear prediction methods for the slope effect ( ) on the initial stiffness ( ) of p-y curves of the undrained laterally loaded pile. Five cases are classified to analyze the slope effect ( ) on the initial stiffness ( ) of a laterally loaded pile, based on the pile-slope position relationship. A series of FEA simulations were conducted, and the data were derived to calculate . Nonlinear models of were established for these five cases, and parameters in the models were calibrated using the results of FEA. The nonlinear models showed good agreement with FEA results. Compared with other researchers' models, the method in this paper is in a reasonable range and can predict more cases. The test cases were used to validate the proposed model, which shows that the p-y curves are in better agreement with test results when using the model of the present study. Additionally, some other conclusions were also obtained: (1) The pile-soil response in undrained clay conditions was nearly independent of the dimensionless pile head-crest distance (ℎ/ ) for the "no slope geometry effect" case, which can be explained by the theory of Coulomb earth pressure. (2) The depth of the first turning point of the pile under undrained loading was explored as 0.8 for the flexible pile, 0.7 for the elastic pile, and 0.6 for the flexible pile. The soil-pile response for the case of ℎ ≥ / was found to be the same as the case of ℎ ≥ 1 ("no slope geometry effect" case).
(3) Compared with the available initial stiffness ( ) calculation methods, the results are more reasonable for large slope angle cases, and can predict some ideal cases, such as the case of = 90°, by using the method of the present study.

Conclusions
The paper proposes nonlinear prediction methods for the slope effect (µ) on the initial stiffness (K i ) of p-y curves of the undrained laterally loaded pile. Five cases are classified to analyze the slope effect (µ) on the initial stiffness (K i ) of a laterally loaded pile, based on the pile-slope position relationship. A series of FEA simulations were conducted, and the data were derived to calculate µ. Nonlinear models of µ were established for these five cases, and parameters in the models were calibrated using the results of FEA. The nonlinear models showed good agreement with FEA results. Compared with other researchers' µ models, the method in this paper is in a reasonable range and can predict more cases. The test cases were used to validate the proposed model, which shows that the p-y curves are in better agreement with test results when using the model of the present study. Additionally, some other conclusions were also obtained: (1) The pile-soil response in undrained clay conditions was nearly independent of the dimensionless pile head-crest distance (h/L) for the "no slope geometry effect" case, which can be explained by the theory of Coulomb earth pressure. (2) The depth of the first turning point of the pile under undrained loading was explored as 0.8L for the flexible pile, 0.7L for the elastic pile, and 0.6L flex for the flexible pile. The soil-pile response for the case of h p ≥ z cr /L was found to be the same as the case of h p ≥ 1 ("no slope geometry effect" case). (3) Compared with the available initial stiffness (µ) calculation methods, the results are more reasonable for large slope angle cases, and can predict some ideal cases, such as the case of θ = 90 • , by using the method of the present study. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

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

Appendix A. Theoretical Method of Georgiadis and Georgiadis [18] and Yang et al. [16]
For undrained laterally loaded piles, the ultimate lateral load per unit pile length (p u ) in Equation (1) is calculated by: where N p was calculated by Georgiadis and Georgiadis [18] for piles at the slope crest: and calculated by Yang et al. [16] for piles in the middle of the slope, as following: where α 1 = 1/(1 + tan θ), α 2 = 1 − sin θ(1 + sin θ)/2; and N pu , N p0 , and λ can be calculated as follows: N pu = π + 2∆ + 2 cos ∆ + 4(cos where ∆ = 1/ sin α and α is the adhesion factor, which is related to c u and can be determined by figuring out α versus c u relationships from Georgiadis and Georgiadis' [18] study. The initial stiffness K i in Equation (1) is calculated by: The linear reduction factor proposed by Georgiadis and Georgiadis [18] is given as: It is given by Yang et al. [31] as:  Figure A1. Process of deriving from simulation results. Figure A1. Process of deriving K iθ from simulation results.