On the External Failure Surface in PSEs: Numerical and Theoretical Methods

e evolution of soil arching, involving internal and external failure surfaces, is of signicance to the load transfer mechanism of pile-supported embankments (PSEs), and external failure surfaces are generally observed in cases of greater embankment height. In this study, the evolution of external failure surfaces was investigated by both numerical and theoretical methods. To begin with, numerical simulations of trapdoor tests were carried out by a two-dimensional discrete element method.e inuences of two key parameters (i.e., embankment height and net spacing between piles) on the development of external failure surfaces were emphasized. e measured coecient of lateral earth pressure around the external failure surface was close to the predicted coecient of active earth pressure. en, a theoretical solution considering inclined external failure surfaces, more realistic compared to those adopted in Terzaghi’s method, was proposed. Compared with Terzaghi’s method, the proposed solution exhibited better consistency with laboratory observations, especially when external failure surfaces were signicantly observed.


Introduction
Pile-supported embankments have been widely used to treat the highly compressible ground in the construction of expressways [1][2][3] and high-speed railways [4,5] owing to their superiority in reducing postconstruction settlement and shortening the construction period. Soil arching is the most fundamental load transfer mechanism of the pile-supported embankment. To analyze the soil arching e ect, numerous analysis models have been proposed, and generally, they can be divided into three groups. e rst group is based on the Marston theory [6] and the Terzaghi theory [7], considering the equilibrium of the soil mass. e second group considers the limit equilibrium theory of soil arch [8]. e third group [9] assumes a xed arch in the embankment, where the embankment loads above the arch are directly transmitted to the pile and the embankment loads under the arch are applied to the soft soil. However, the results of the same project using di erent models vary greatly, indicating that the analytical models heavily rely on the assumed patterns of arching evolution. e mechanism of arching evolution has been commonly investigated by trapdoor tests [10][11][12][13][14], originally designed by Terzaghi [7]. By this means, soil arching is generally recognized to be evolved with various key parameters, such as embankment height (h), pile diameter (b), and pile spacing (s) [15][16][17][18][19]. In some cases, the failure surface may extend beyond the trapdoor to form an external failure surface. Since vertical failure surfaces were assumed for soilpile stress distribution in Terzaghi's [7] were oversimpli ed, a series of trapdoor tests were performed by Ladanyi and Hoyaux [20] to check Terzaghi's method. A great discrepancy between theoretical results and laboratory observations can be observed when the external failure surfaces are formed at a great magnitude of trapdoor displacement. For this reason, di erent failure surfaces in embankments, especially the external failure surface at great trapdoor displacement, can be further investigated.
Several experimental investigations have been carried out on the external failure surface [21][22][23][24]. In particular, Costa et al. [21] reported that the inclination of external failure surfaces to the horizontal would not be less than the repose angle of the soil, which represents the angle of shearing resistance of soil at its loosest state. eoretical models [25][26][27] with sliding surfaces of arbitrary inclinations are also proposed to examine the effect of the sliding surface on vertical stress. Although numerous numerical simulations have been presented on soil arching in PSEs [18,19,28,29], the discrete element method (DEM) has scarcely been adopted to study the microscopic characteristics of the external failure surfaces.
In the present study, trapdoor tests were simulated by the particle flow software PFC 2D , coded by the discrete element method (DEM) [30]. e DEM can provide better insight into the microbehavior of the embankment as a result of the discrete nature of the granular embankment fill materials. On this basis, more attention will be paid to the influences of embankment height (h) and pile net spacing (s-b) on the evolution of external failure surfaces and the stress state near the external failure surfaces. In parallel, a simple theoretical model considering the external failure surfaces would be proposed, and then it would be compared with Terzaghi's method.

Description of Reference Case.
A numerical model by PFC 2D has been established in this study based on the laboratory tests by Jenck et al. [31] and the configuration of the laboratory model was shown in Figure 1. e embankment fill is represented by steel rods of 3, 4, and 5 mm in diameter (1 : 1:1 by volume) and the porosity of the assembly of steel bars is equal to 0.18. e unit weight and internal friction angle of the assembly are 62 kN/m 3 and 24°, respectively. e rigid pile and the soft soil are simulated by rigid elements and foam, respectively.

Basic Kinetics Used in PFC.
e sphere (particle) element and the wall element are the basic components used in the PFC model, and the interactions between them are represented by the contact force. e contact force applied to the particles will produce corresponding particle movement as per the motion equation; however, the one applied to the wall would not activate the motion of the wall.
Without the constraint of deformation coordination applied between particles, the equilibrium equation governs the behavior of the PFC model. If the resultant force or moment of the particles is not zero, the unbalanced force or moment will lead to the movement of particles according to Newton's law. e movement of particles is thus dependent on the contact forces of adjacent particles, and the relationship between contact force and resultant displacement is called the physical equation. e calculation of PFC is based on Newton's second law and force-displacement law, where Newton's second law applied to each particle is used to update the position of the particle and the wall. e forcedisplacement law is used to update the contact force among particles at each contact. Each particle will undergo a stepby-step iteration with updated position and contact force as per the kinetics, and the iterations will be terminated when the force equilibrium is reached.
(1) Force-displacement law (physical equation) e force-displacement law is applied at the contact to update the contact force. e contact force, however, only works when it is less than the activity distance, and the contact gap between particles should be identified first. e contact force F i can be decomposed into the normal component F i n perpendicular to the contact plane and the tangential component force F i s in the contact plane.
(2) Newton's second law (equation of motion) If the resultant force or moment acting on the particles is not zero, translational motion or rotation will be produced accordingly. e translational motion of the particle is described by the linear velocity of the particle center, while the rotation of the particle is described by the angular velocity. e above motions are governed by the following equations:  Advances in Civil Engineering where F i is the unbalanced resultant force acting on the particles; m is the particle mass; € x is the acceleration of the particle; gi is the gravitational acceleration; M i is the unbalanced resultant moment on the particles; I is the moment of inertia of the sphere; and _ ω i is the angular acceleration. Subsequently, the accelerations of particles can be obtained by solving equations (1) and (2) at the time step, Δt, and are then used to update the velocity and angular velocity of the particles of the next time step.

e Contact Constitutive Model.
In PFC simulations, the macroscopic mechanical properties of materials are mainly governed by the contact between microparticles, or the particle contact constitutive model. e linear contact model is used in this study because it is applicable for most analysis scenarios and also by using this model the required microparameters can be easily obtained. e mechanical model of particle contact is shown in Figure 2.
Normal and tangential stiffness are, respectively, defined in PFC, and the stiffness corresponds with the particle stiffnesses acting in series. e normal stiffness, k n, and tangential (shear) stiffness, k s, can thus be estimated as follows: where k [A] n and k [B] n are the normal stiffness of particles A and B, respectively, whereas k [A] s and k [B] s are the tangential stiffness of particles A and B, respectively.

2.2.3.
Damping. In this model, the local damping is always active with the default value of 0.7.

Servo Control.
Because the stress cannot be directly endowed to activate the wall in PFC, confining pressure or loading is usually applied by displacing the wall through the servo control.

Calibration of Embankment Fills.
e linear-based model was adopted in the DEM simulation to represent the contact behaviors among embankment fills. ree major parameters were required for the linear contact model: the normal stiffness of particles, k n ; the shear stiffness of particles, k s ; and the friction coefficient, μ _s . eir magnitudes were first roughly estimated based on the recommendations by Stahl and Konietzky [32] and Han et al. [28], and then they were precisely calibrated by using a numerical biaxial test before modeling trapdoor tests. A calibration chamber with a size of 200 mm × 200 mm was built. According to the observation [26], samples were simulated by disk particles with three different diameters (3, 4, and 5 mm) at a volume ratio of 1 : 1:1 and a porosity of 0.18. e servo-controlled confining stress was applied to the embankment fill in the horizontal direction, and the vertical stress was controlled by displacing the horizontal wall at a constant velocity. As shown in Figure 3(a), the curves of deviatoric stress versus the axial strain are obtained at three different confining stresses (i.e., 20, 30, and 40 kPa). It can be observed that numerical simulations and experimental results were in good agreement. Shear stress versus normal stress is shown in Figure 3(b). e internal friction angle of the embankment fill φ can be calculated as follows: where σ 1 and σ 3 are maximum and minimum principal stresses, respectively. e value of the internal friction angle is 24.4°, which is close to the experimental parameters (i.e., 24°), and the other micromechanical parameters of the particles are listed in Table 1.

Numerical Filling Technique.
Numerical simulation techniques to generate embankment fills have been improved gradually in recent years [28,29,33]. e embankments primitively generated by the "radius expansion" and "falling rain" methods were unrealistic. e distribution of the fills was nonuniform and the vertical stresses did not vary with the density gradient along the depth. To ensure the natural stress state, Bhandari and Han [33], and Han et al. [28] numerically compacted the embankment fills through the servo mechanism of a moving wall. Lai et al. [29] proposed a multilayer compaction method that can better simulate practical embankment filling processes and can guarantee the natural stress state. A new problem, however, would inherently arise if the embankment fills were compacted to the target porosity. en, the density of the region near the moving wall would be denser than that of the region away from the wall. In this study, overcompaction and release were used to improve the multilayer compaction method. In this study, the initial porosity n i and target n r were 0.59 and 0.18, respectively, and the thickness of each compaction layer of the embankment was 100 mm, which was consistent with the compaction process in experimental tests [31]. e flowchart of the improved multilayer compaction method is shown in Figure 4, and the detailed steps are as follows.
(1) Particle generation. Particles are generated in an area (area � s×200 mm) according to the calculated porosity n i and the cycle has sufficient time for uniform distribution of particles. (2) Overcompression. e particles are compressed by moving down the top platen to a height of less than 100 mm, and the particles in the area are thus compressed to be overlapped within a transitional porosity n g smaller than n r . (3) Normal compression. e platen is moved up to the designed position (h � 100 mm). e particles are then bounced off by the contact force between the overlapped particles.

Advances in Civil Engineering
(4) e second layer is generated by following the steps from 1 to 3. (5) e platen between the layers is removed and it is cycled to the equilibrium.
In this way, the embankment can be generated in lifts without gravity (g � 0) and friction (μ s � 0). e embankment generated by this technique would be uniform, compacted, and can achieve the target porosity. After imposing gravity on the generated embankment, Figure 4 illustrates the accuracy of the embankment stress state generated by the improved multilayer compaction method. By fitting the numerical prediction in Figure 5, the slope of the numerical prediction curve, indicating the unit weight used in the simulation is 64.33, which is close to the magnitude of the unit weight in the experiment (c � 62 kN/m 3 ).   (s) of 320, 450, and 650 mm were considered in reference cases, and the embankment height (h) increased from 100 mm to 600 mm for each pile spacing. A movable wall located at the bottom of the embankment simulates the trapdoor, and its movement was denoted as trapdoor displacement (δ). In the numerical simulation, the moving rate should be slow enough to guarantee the quasistatic equilibrium of the sample during the test [34,35], and a constant vertical velocity of 1 mm/s was applied to the trapdoor. In this model, the default time step for each calculation cycle was set as 2.1 × 10 −5 s and the corresponding moving rate was equivalent to 2.1 × 10 −5 mm/step, so that it takes more than 46,000 steps to move 1 mm. At the same time, a series of measurement circles were set to evaluate the performance of the embankment.

Validation of the Numerical Model.
e efficacy of loads distribution, E, is used to characterize the soil arching effect, and it can be calculated as follows [8]: where F p is the load carried by piles and W is the total weight of the embankment. e results of the DEM simulations were compared with the experimental results for validation. As shown in Figure 7, the efficacy E increases with the embankment height and decreases with the increase of pilenet spacing. It can be found that the results of the DEM simulations match well with the experimental results.

Evolution of External Failure
Surfaces. e evolution of soil arching can be described by the development of failure surfaces, which are generally explored by the displacement contour [18,19,21,36], and a typical case is illustrated in Figure 7 (s � 320 mm, h � 400 mm). Figure 8(a) demonstrates the development of soil displacement during the displacement of the trapdoor. When the trapdoor reaches the designated displacement of 17.25 mm, the embankment is divided into three regions according to the displacement contour, which includes the active region, transition region, and stable region. According to Costa et al. [21], the boundary between the transition region and active region (curve OA) represents the internal surface and the boundary between the transition region and stable region (curve OB) represents the external surface. e active region will flow downward with the displacement of the trapdoor, making the transition region move obliquely downward to the active region accordingly.  Advances in Civil Engineering 5 In this study, the failure surfaces are outlined in Figure 9. To facilitate the description of the failure surface, simplified straight failure surfaces were adopted to approximately represent the actual failure surfaces, and they share the same starting and ending points. e inclination of the simplified failure surface, θ, is the angle between the line of the failure surface and the horizontal plane, as shown in Figure 9. e external failure surfaces are only observed when the embankment height, h, is larger than the pile-net spacing (s-b). For example, for cases with a net spacing of 220 mm, the external failure surface would occur when the embankment heights are 300, 400, 500, and 600 mm. It is also true for the cases with the net spacing of 350 mm when embankment heights are 400, 500, and 600 mm, and the cases with the net spacing of 550 mm when the embankment height is 600 mm. It should be noted that the inclination of the failure surface (θ) decreases gradually with the embankment height.

Distribution of Vertical Stress.
e contact force between particles can be monitored by measurement circles and then can be converted into continuum stress based on the following equation: where N p1 is the number of particles whose centroids are contained within the measurement circle; N c is the number of contacts between particles; n is the porosity within the measurement circle; S (p) is the area of a particle; s i (p) and s i (c) are the locations of a particle centroid and its contact, respectively; n i (c,p) is the unit normal vector directed from a particle centroid to its contact location; F j (c) is the force acting on the contact; and indices i and j are in the set of (x, z).
e vertical stress, σ v , can be calculated based on the following equation for i � z and j � z: where σ zz is the average stress in the direction of the z-axis. Figure 10 shows the distribution of vertical stress in the embankment before and after displacing the trapdoor. e stress redistribution results in the prominent increase of the loads acting on piles increasing from 25 kPa to 60 kPa, but the reduction of loads on the soil decreasing from 25 kPa to      Advances in Civil Engineering 9.5 kPa with the trapdoor displacement.
e areas of increased vertical stress and corresponding equipotential lines of the contour roughly coincide with the stable regions and external failure surfaces in Figure 8, respectively. It implies that more loads of surrounding soil are transferred to the stable regions through the external failure surface, thus increasing the vertical stress on the pile top. e ratio of the horizontal stress to vertical stress, defined as the coefficient of lateral earth pressure (k), is an important parameter and it is calculated before and after displacing the trapdoor as shown in Figure 11. e coefficient of lateral earth pressure around the external failure surface is 0.48, roughly approaching the coefficient of active earth pressure, k a � tan 2 (45°-φ/2) � 0.42. On this basis, the coefficient of active earth pressure can be used to represent the coefficient of lateral earth pressure around the external failure surface in the proposed method.

eoretical Model Establishment.
It is essential to establish a theoretical model for calculating the efficacy of pilesupported embankments with external failure surfaces, and the schematic of the theoretical model is shown in Figure 12(a). e coordination system is established with the origin located at the center of the pile top plane, and the right direction and upward direction are denoted as positive, respectively. Considering the elastic-plastic behavior of the soil commonly described by the Mohr-Coulomb model [37], two simplified straight failure surfaces AB and A′B′ are assumed to deviate with the displacement direction (i.e., zaxis) at an angle of α (α � π/4 + φ/2), as shown in Figure 12(a). e pile-supported embankment is a threedimensional (3D) problem; however, the 2D plane-strain model is much easier to solve and can also obtain reasonable accuracy for engineering practices [38,39]. In this study, the sliding soil over the trapdoor is sliced into segments along the depth, and the coefficient of active earth pressure is used in the proposed method based on the embankment behavior in the numerical simulation.
Considering the vertical force equilibrium of the segment with the thickness of dz, shown in Figure 12(b), the following equation can be derived: where P s (z) is the vertical stress at the depth of z, F n (z) is the stress perpendicular to the failure surface at a depth of z, D is the lower width of the soil unit, α is the angle of the failure surface with a z-axis, φ is the internal friction of soil, and c is the unit weight of soil. For the horizontal force equilibrium of the differential segment dz, the following equation can be obtained: where k a is the coefficient of active earth pressure, equal to tan 2 (45°-φ/2). By rearranging equation (10), the following relationship can be obtained: where a 1 � cosα-tanφsinα. According to geometrical relationships, we can obtain the following formulas: dD � 2 tan α dz.  Substituting equations (11)-(13) into (9), the following equation stands: e boundary condition should satisfy P s (z) can then be solved by combining equations (14) and (15).
For the equilibrium of vertical stress at z � 0, we obtain the following equation: where m is the area replacement ratio (m � b/s), P s|z � 0 is the embankment load suffered by the soil, and P p|z � 0 is the embankment load suffered by the pile, which can be solved by using equation (16). Expression of efficacy E can be obtained from equation (17) and is displayed as

Comparative
Analyses. e parameters used in the comparative analyses are consistent with the experimental tests [31]: the pile net spacings (s-b) are 350 mm and 550 mm, the internal friction of soil φ is 24°, and the weight of soil c is 62 kN/m 3 . e calculated results of the efficiency E from the proposed method and the comparison between the experimental results and those from Terzaghi's method [7] are shown in Figure 12. e method for predicting the efficiency E suggested by Terzaghi [7] can be briefly clarified by the following formulas: When h ≤ 2s When h > 2s where k, equals 1, as suggested by Terzaghi. As shown in Figure 13, the results calculated by Terzaghi's method show good agreement with the experimental results when the embankment height, h, is less than the pilenet spacing, s-b. As concluded from the numerical simulation, the embankment failure surface is vertical when the embankment height is less than that of the pile-net spacing and the shape of the failure surface is consistent with the assumption proposed in Terzaghi's method. Nevertheless, the significant differences from the assumption can be observed when the embankment height is larger than the net spacing because the failure surface no more extends vertically along the edge of the trapdoor but it extends beyond the trapdoor. As a result, the results calculated by the proposed method would be more accurate for high embankments with external failure surfaces compared to those by Terzaghi's method with the simple vertical failure surface. In conclusion, Terzaghi's method is appropriate when the embankment height is less than the net spacing, whereas the proposed method is appropriate when the embankment height is larger than the net spacing.

Conclusion
e evolution of soil arching in pile-supported embankments plays an important role in the selection of arching models, and it evolves with embankment height, pile diameter, and pile spacing. In addition to internal failure surfaces, there are also external failure surfaces in some cases of the pile-supported embankment. Numerical models were developed in this study to analyze the external failure surfaces. To accurately replicate the initial state of the embankment, an improved multilayer compaction method was adopted so that the porosity of the generated particles in the PFC approaches the realistic embankment fill. e micromechanical parameters of embankment fills were calibrated by numerical biaxial tests at various confining pressures, and the numerical model was then validated by comparing the  efficacy with the case observation. It is concluded that the external failure surface will only occur when the embankment height is greater than the pile-net spacing and the embankment load is redistributed towards the pile head through the external failure surface. e coefficient of lateral earth pressure around the external failure surface is 0.48, and it is close to the coefficient of active earth pressure. e calculation results indicate that the proposed solution demonstrates better consistency with experimental results than that of Terzaghi's method for the case with external failure surfaces.