#### Abstract

This study investigates the lateral behavior of pile foundations socketed into bedrocks using 3D finite difference analysis. The lateral load-displacement curve, pile deflection, and bending moment distribution were obtained for different bedrock depths between 3 and 20 m. It was discovered that bedrocks that have a depth of 7 m (7*D*) or less influence the lateral behavior of the pile. The *p*-*y* curves were collected at depths of 2.0–4.5 m, and the effect of the bedrock on the curves was evaluated. It was observed that the *p*-*y* curves were significantly affected by the material properties of the bedrock if the rock is located in close proximity (within 3D), but the effect is diminished if the *p*-*y* curves were 3.5 m (3.5*D*) or farther from the bedrock.

#### 1. Introduction

The subgrade reaction method was developed by the Winkler assumption [1], where the soil resistance per unit length, *p*, is assumed to be proportional to the pile deflection, *y*, at a certain point. The linear relationship between *p* and *y* does not represent the nonlinear characteristics of the soil. To address this phenomenon, the nonlinear *p*-*y* method was introduced by McClelland and Focht [2]. Since the introduction of this method, the *p*-*y* curves are typically used to indicate the nonlinear springs representing lateral soil resistance and deflection along the pile. In conventional designs, the soil supporting the pile can be replaced with a set of nonlinear springs, and the pile is assumed to behave as a beam supported by nonlinear springs. The *p*-*y* method is widely accepted for evaluating the lateral behavior of piles due to its simplicity. Numerous researchers have proposed different *p*-*y* curves with considerations for soil types, loading type (static or cyclic), and water conditions (saturated or unsaturated soil). In the early stages of *p*-*y* curve development, the *p*-*y* curve for submerged soft clay was developed by Matlock [3], stiff clay with free water by Reese et al. [4], stiff clay without free water by Reese and Welch [5], and submerged sand by Reese et al. [6]. The *p*-*y* curves for sand developed by Reese et al. [6] were modified into a simpler equation that yielded better estimation by Murchison and O’Neill [7]. The curve of Murchison and O’Neill was later adopted by the American Petroleum Institute [8] for the design of offshore platforms.

There are limited analytical researches available for evaluating the lateral response of piles socketed into rocks. Carter and Kulhawy [9] investigated the lateral behavior of drilled shafts (flexible and rigid type) socketed into rocks using finite element methods. Closed-form equations were proposed considering factors such as loading conditions, material properties, and rock mass stiffness. Zhang et al. [10] proposed a nonlinear continuum method to predict the lateral response of drilled shafts, in which the ultimate resistance of a rock mass was calculated based on the Hoek–Brown failure criterion. To et al. [11] believed the *p*-*y* curves are inappropriate for estimating the behavior of jointed rock masses as the curves consider soil as a continuum. Thus, they proposed a discontinuum model based on rock mass with two or three joint sets. Recently, the elastoplastic continuum method was proposed considering the variation in the flexural rigidity of drilled shafts in multilayers of soil and rock masses [12].

Numerical simulations have been performed extensively to account for anisotropic and discontinuous characteristics of rock medium [13–18]. The three-dimensional distinct element method was adopted using 3DEC to simulate the jointed rock mass under lateral loading, and *p*-*y* curves were proposed for mudstone [13–15]. It was discovered that the number of joints significantly affects the lateral response of piles in rock masses. Liang and Shatnawi [16] performed a series of finite element analyses, presented charts for the initial modulus of subgrade reaction for transversely isotropic rock media, and recommended using five elastic constants for estimating the initial modulus of subgrade reactions. Comparisons were made on the lateral behavior of drilled shafts in the rock medium for isotropic properties following the Hoek–Brown criterion and for transverse isotropic rock properties [17]. Shantnawi [19] expanded the *p*-*y* criteria considering the anisotropic behavior as well as discontinuous joint sets in rock masses.

The *p*-*y* curves for rocks are very difficult to obtain and validate through experiments because the full-scale tests on rocks are costly and because it is difficult for the rock to fail within the available capacity of the loading system. Indeed, the number of full-scale tests on drilled shafts in rocks was very limited [20, 21]. The complete form of *p*-*y* curves for weak rocks had not been available until that reported by Reese [22]. He proposed *p*-*y* curves based on the results of the two load tests on coral limestone conducted by Nyman [21]. Reese [22] commented that the proposed *p*-*y* curves in weak rock should be used with caution due to the limited number of test results. Gabr et al. [23] proposed hyperbolic *p*-*y* curves for weak rocks based on their field test results on small-diameter drilled shafts. The curves require two parameters: ultimate resistance and initial modulus of the subgrade reaction. However, their curve was criticized as basing their estimation of the initial modulus of subgrade reaction in the study of Vesic [24] that is inappropriate for piles socketed into rocks [25]. Numerous other researches indicate the inappropriateness of conventional *p*-*y* curves for different types of rock medium [26–28]. *p*-*y* curves for the Ohio shale was proposed using the results of field lateral load test on drilled shafts with a diameter of 1.8 m socketed into shale [29]. According to the test results, the *p*-*y* curves of weak rocks [22] were found to underestimate the pile deflection. The same research group developed *p*-*y* curves for rock and intermediate geomaterials using in situ pressure meter tests [30] and *p*-*y* curves for transverse isotropic rock [31, 32]. Cho et al. [26] conducted six field tests on drilled shafts, proposed *p*-*y* curves for soft weathered rock, and concluded that the Reese *p*-*y* curves [22] overestimate the resistance. Conversely, the stiff clay model [5] underestimates the resistance for the same pile deflection. Liang et al. [33, 34] recommended hyperbolic *p*-*y* curves for rock masses using numerical simulations, where the curve was validated through two full-scale load tests on drilled shafts with a diameter of 1.3 m. Recently, centrifuge model tests were employed to simulate large-diameter drilled shafts in rocks [35, 36]. *p*-*y* curves for weakly cemented soils to well cemented soils were proposed by Guo [36]. Full lateral load tests were performed on drilled shafts in weak calcareous sandstones, and centrifuge model tests were performed for preinstalled piles in very weakly cemented sands. Parsons et al. [37] introduced a case history of full-scale cyclic lateral load tests on two drilled shafts rock-socketed in limestone.

The purpose of this study is to investigate the effect of bedrock depth on the lateral behavior of piles socketed into bedrocks. Due to the significant difference in material stiffness and strength between soils and bedrocks, the soil-pile interaction near the boundary between the soil and rock is very complex to understand. Conventional *p*-*y* curves of soils can be created from the material properties of the soils, not from the nearby bedrocks. However, the *p*-*y* curves of soils are likely to be considerably affected by the material properties of the rock if the rock is located in close proximity, implying that the *p*-*y* curves obtained near the bedrock would cause the overestimation of the ultimate soil resistance and the modulus of the subgrade reaction, which is misleading. The three-dimensional finite difference method was adopted to derive the *p*-*y* curves of soils overlying rocks, while the variations in the ultimate soil resistance and the initial modulus of the subgrade reaction at varying distances from the rock were obtained.

#### 2. Numerical Modeling

Figure 1 illustrates the numerical model along with the stratigraphy, pile dimension, socket depth into bedrock, locations where the *p*-*y* curves were collected, variable for the parametric study, 3D model in FLAC3D, and boundary conditions. As shown in Figure 1(a), the pile is socketed 3 m into bedrock that is underlain by sand with a thickness of *H*. The pile head at the ground level was laterally pushed over, and the *p*-*y* curves were collected every 0.5 m from 2.0 to 4.5 m from the ground surface. The distance from the locations of the *p*-*y* curves to the bedrock (*X*_{py}) varied with varying thicknesses of the sand layer (*H*). For example, when *H* is 5 m, the distances, *X*_{py}, are 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 m. For the same thickness, *H*, the *p*-*y* curves of sand at an identical depth would differ due to the influence of varying bedrock depths.

**(a)**

**(b)**

**(c)**

The simulated pile was 1 m in diameter and 20 mm in thickness and embedded 3 m into the bedrock layer. The length of the simulated pile varied from 6 to 23 m, and the pile was created as a so-called “wish-in place” pile. The side boundary was located 10 m from the center of the pile and fixed in the *x* and *y* directions against translation. In addition, the bottom boundary was 7 m below the pile tip and fixed in three directions against translation. The pile head was laterally pushed over with a velocity of 10^{−6} mm per step up to 500 mm in the lateral direction. The applied lateral displacement was far greater than 38 mm that is considered the serviceability limit in lateral [38]. This is because the 38 mm displacement did not induce sufficient pile deflection to generate the *p*-*y* curves. The lateral displacement, lateral load, and moment profiles were obtained and used to calculate the soil resistance per unit length, *p*, through the double integration of the bending moment profile and the pile deflection, *y*, by taking the second derivative of the moment profile.

Table 1 shows the material input parameters used in the numerical study. The sand and bedrock were modeled using the Mohr–Coulomb model, but the model for the bedrock was reduced to the Tresca model because the internal friction angle was set to zero. The friction angle of the sand was 35.9° while the elastic modulus was 42.0 MPa. A cohesion of 27.6 MPa was assigned to the bedrock with zero friction angle. The elastic modulus of sand was assumed to be constant along the depth for simplification. This may reduce the lateral displacement of the pile near the ground surface. The pile was assumed to have linear elastic behavior with an elastic modulus of 210 GPa and Poisson’s ratio of 0.3. Because the pile was modeled as a solid instead of tabular pile, the elastic modulus of the pile was back-calculated by matching the flexural rigidity (EI). The input value for the elastic modulus was 31.64 GPa, resulting in an identical flexural rigidity of 1,550 MN·m^{4}.

The slip and separate only model was adopted to simulate the pile-soil interaction; the interface element of the model was characterized by the Coulomb friction model in the shear direction and by tensile strength in the normal direction. In the Coulomb friction model, the shear resistance linearly increases up to the maximum side shear with increasing shear displacement, and after the maximum shear, the shear resistance remains constant against further displacement. For the Coulomb input parameters, the friction angle and the cohesion of the surrounding materials were used. The shear and normal stiffnesses of the interface element were set to ten times the equivalent stiffness of the surrounding soil as recommended in the FLAC3D manual [39].

In numerical simulations, the numerical results should not be affected by the mesh size, boundary condition, or interface element of the pile-soil interaction. The load-displacement curves at the pile head were collected to evaluate the effect of these factors. Figure 2 illustrates the effect of mesh size on the numerical results. The lateral load-displacement curves were obtained for different mesh sizes: 0.10, 0.25, 0.50, 1.00, 2.00, and 4.00 m. More lateral load was required to displace the pile head for larger mesh sizes, and the required load drastically dropped as the mesh size decreased. Figure 2(b) shows the variation in the lateral load for a 38 mm pile displacement, indicating that the effect of mesh size appears to be negligible when the mesh size is less than 0.25 m. Therefore, a mesh size of 0.25 m was used in the numerical models.

**(a)**

**(b)**

The boundary should be located sufficiently far from the pile so that the stress transferred from the laterally loaded pile to the soil should not reach the boundary. Figure 3 shows the load-displacement curves for four cases used to evaluate the effect of boundary conditions on the numerical results. The boundaries were located 5 to 20 m away from the center of the pile in the *x* direction. The boundary in the *x* direction plays a key role because the pile was displaced in the *x* direction. It was found that the lateral load-displacement curve did not significantly differ when the boundary was placed 10 m or farther in the *x* direction. Thus, the boundaries were created at 10 times the pile diameter (10*D*) away from the center of the pile in both the *x* and *y* directions.

The stiffness of the interface element (*E*_{i}) simulating the pile-soil interaction was determined by the stiffness of the surrounding soils (*E*_{s}). Figure 4 illustrates the effect of the stiffness of the interface element on the lateral response of the pile. The stiffness was varied from 1 to 20 times that of the surrounding soils, and the obtained load-displacement curves indicate that the stiffness of the interface element does not affect the numerical results if the stiffness was set to 5 times the stiffness of the surrounding soil or higher. Thus, the value of 10 times the stiffness of the surrounding soil recommended in the FLAC3D manual was used [39].

#### 3. Results and Discussions

A total of six numerical simulations were conducted with varying bedrock depths of 3, 5, 7, 10, 15, and 20 m. Figure 5 presents the lateral load-displacement curves for the six simulations, showing that a bedrock depth of 3 m results in dramatic increase in the lateral capacity of the pile. It should be noted that the pile was assumed to be linear elastic with no failure, but in reality, the calculated lateral capacity could cause material failure. The assumption was made because the study focuses on the soil response, not on the material failure. The effect of the bedrock gradually fades away and appears to be nonexistent when the bedrock depth is 10 m or greater.

The same observation was made in the pile deflection and bending moment profiles, as shown in Figure 6. The pile deflection was strongly confined by the surrounding bedrock at depths of 7 m or shallower, but no meaningful difference was noticed at greater depths. The maximum bending moment occurred near the bedrock when it was at a shallow depth, but the bending moment profiles for bedrock depths of 10 m or greater appeared to be identical. From the deflection and bending moment profiles, it could be concluded that the bedrock had no effect at depths of 10 m or greater.

**(a)**

**(b)**

Figure 7 shows the variation in the required lateral load and the maximum bending moment when the pile head was displaced by 38 mm. The lateral load was about 4,800 kN for a bedrock depth of 3 m, dropping to about 820 kN at a depth of 10 m, and then remaining constant at greater depths. The maximum bending moment was also greatest for a bedrock depth of 3 m, drastically declined to a depth of 10 m, and remained constant at depths below 10 m.

**(a)**

**(b)**

The lateral soil resistance per unit length, *p*, can be calculated by dividing the sum of the lateral forces applied to the entire interface nodes per unit length along the pile, while the lateral displacement, *y*, can be obtained from the nodal displacements. The available FISH language in FLAC3D was used to calculate the *p*-*y* curves, but the lateral soil resistance per unit length, *p*, was calculated differently. Figure 8 shows the stresses at the pile-soil interface. The *x*-directional (horizontal) component of the normal stress (*σ*_{x,i}) and shear stress (*τ*_{x,i}) at interface node *i* was calculated using Equation (1), and the lateral soil resistance, *p*, was calculated by summing all *p*_{i} at the desired elevation. The detailed process for calculating the *p*-*y* curves can be found in previous researches [40, 41]where is the lateral soil resistance at the interface node , is the normal stress at the interface node , is the shear stress at the interface node , is the representative area of the interface node , and is the pile length between two neighboring interface nodes along the pile.

Figure 9 shows the *p*-*y* curves for sandy soil at six elevations from the ground surface: 2.0, 2.5, 3.0, 3.5, 4.0, and 4.5 m. The relatively shallow depths ranging from 2.0 to 4.5 m were selected because the soil response near the ground surface is the most important factor affecting the lateral behavior of piles. At depths of 2.0, 2.5, 3.0, and 3.5 m, the *p*-*y* curves at the same elevation agree well with each other irrespective of the bedrock depth (*H*). In other words, the lateral soil response at a bedrock depth of 5 m is not different from that at a bedrock depth of 20 m. Conversely, as the location of the *p*-*y* curves gets closer to the bedrock, the *p*-*y* curves at the similar depths become significantly different. The *p*-*y* curves at a depth of 4.0 m displayed much stiffer moduli of subgrade reaction when the bedrock depths were 5.0 and 7.0 m (i.e., 1.0 and 3.0 m to bedrock) than those of other depths (i.e., distance to bedrock greater than 6.0 m). Similar behavior was observed at a depth of 4.5 m. This finding implies that the *p*-*y* curves obtained near the bedrock (as close as 3*D* (3.0 m)) do not properly represent the soil response. In fact, the obtained curves in such area overpredict the stiffness (i.e., modulus of the subgrade reaction) of the *p*-*y* curves compared to those not affected by the bedrock. The stiffer behavior of the soil near the bedrock is likely due to the higher strength and stiffness of the bedrock; the elastic modulus of the bedrock is 500 times that of the soil or higher. The contribution of the bedrock stiffness becomes greater as the obtained *p*-*y* curves gets closer to the bedrock. The soil resistance did not reach the ultimate value for those *p*-*y* curves (Figures 9(e) and 9(f)), but it is likely to increase near the bedrock. Figure 9 also indicates that the *p*-*y* curves obtained at bedrock depths of 10 m or greater can be assumed to have the similar lateral pile response regardless of the presence of the bedrock. This finding is consistent with the facts that can be inferred from Figures 5–7.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

#### 4. Conclusions

In this study, the effect of bedrock depth on the lateral behavior of piles was investigated using 3D finite difference analysis. The variations in pile deflection profiles as well as bending moment distributions for various bedrock depths were discussed. The ultimate soil resistance and stiffness in the obtained *p*-*y* curves were compared for various bedrock depths. The following conclusions can be drawn:(1)The effect of the bedrock gradually diminished with increasing bedrock depth and eventually disappeared at bedrock depths of 10 m (10*D*) or greater. A lateral load of 821 kN was required to displace the pile head by 38 mm for piles not affected by the bedrock, while a load of 4,803 kN was required at a bedrock depth of 3 m (3*D*).(2)The maximum bending moment occurred near the bedrock appearance for a bedrock depth of 7 m (7*D*) or less but occurred at a depth of about 4.5 m (4.5*D*) for greater bedrock depths. In the simulations, bedrock depths of 10 m (10*D*) or greater showed similar pile deflection and bending moment profiles along the pile, indicating that the effect of bedrock depth becomes negligible at greater depths.(3)It was found that the *p*-*y* curves near the bedrock were strongly influenced by the material properties of the bedrock; there were significant increases in the ultimate soil resistance per unit length and in the stiffness (modulus of subgrade reaction). In other words, the *p*-*y* curves obtained near the bedrock (within 3*D*) do not properly represent the soil response, and the modulus of the subgrade reaction and the ultimate soil resistance may be overestimated.

#### Data Availability

The numerical simulation data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was partially supported by the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning (NRF-2016R1C1B2013478) and by 2018 Hongik University Research Fund.