Study on the Dynamic Soil-Pile-Structure Interactive Behavior in Liqueﬁable Sand by 3D Numerical Simulation

: The dynamic behavior of structures in liqueﬁable sand exhibits more complicated characteristics, due to the development of excess pore pressure caused by cyclic loading, than that in dry sand. Therefore, it is crucial to accurately predict the soil–pile structure behavior during liquefaction to prevent damage to the structures. In this study, three-dimensional numerical modeling was performed to predict the dynamic soil–pile behavior during liquefaction. To directly simulate pore pressure generation due to soil shear deformation, the Finn liquefaction model was applied and coupled with the Mohr-Coulomb elasto-plastic model. Soil nonlinearity was considered by applying hysteretic damping, and the interface model was applied to simulate various dynamic phenomena between the soil and pile. Simpliﬁed continuum modeling was introduced to prevent reﬂection wave generation and increase analysis e ﬃ ciency. The applicability of the proposed numerical model was validated using the experimental results. Thereafter, a parametric study was conducted to provide a better understanding of the dynamic behavior of pile foundation during liquefaction. From a series of parametric studies, several important factors that can a ﬀ ect the dynamic pile responses in liqueﬁable sand were identiﬁed. Also, the characteristics of the dynamic soil–pile structure interactive behavior, which are signiﬁcantly di ﬀ erent from each other in liqueﬁed and dry sand, were analyzed qualitatively and quantitatively.


Introduction
Liquefaction, which generally occurs in loose sandy soil, is one of the major causes of damage to facility structures during earthquakes because many structures, such as civil infrastructure systems and plant facilities, are currently constructed on liquefiable soft ground, such as reclaimed area. In particular, in the magnitude 5.4 Pohang earthquake in November 2017, liquefaction was firstly observed in South Korea, raising concerns for the safety of foundations embedded in liquefiable ground. As the pile behavior in liquefiable sand has a more complicated mechanism than that in non-liquefiable soil, it exhibits much more uncertain and unpredictable characteristics. Not only do the inertial force of the upper structure and the kinematic force due to soil deformation apply different dynamic loads on the pile, but the stiffness and shear strength of the soil also decreases over time due to the nonlinear behavior of soil and development of excess pore pressure. Recently, owing to the development of physical model tests, such as the centrifuge test and 1-g shaking table test, and establishment of the numerical analysis method, several studies have evaluated the dynamic behavior of pile foundations in liquefiable sand. In particular, owing to the temporal, economic, and spatial limitations of physical thereby increasing the pore pressure of the soil. This increase in pore pressure decreases the effective stress of the soil, and liquefaction occurs when the generated excess pore pressure equals the initial effective stress of the soil. In this case, soil completely loses its shear strength and exhibits behavior similar to that of water. Soon, it poses a serious threat to the stability of the structure's foundation. Therefore, it is crucial to accurately predict the structural behavior of foundations in the event of liquefaction to prevent damage to the structures. For this purpose, two types of analysis method are generally employed: total stress analysis and effective stress analysis.

Total Stress Analysis
In the event of an earthquake, the seismic wave propagated from the bedrock reaches upper soil particles through the soil layer, as shown in Figure 1. In this case, liquefaction occurs when the underground cyclic shear stress caused by the earthquake exceeds the magnitude of the cyclic shear stress required to cause liquefaction. The zone of liquefaction is determined by comparing the magnitudes of the generated cyclic shear stress and shear strength according to the depth, as shown in Figure 2 (Seed & Idriss, 1971 [18]). Liquefaction is similar to the strain-softening behavior that occurs when sandy soil is subjected to shear stress under undrained conditions. This is the concept whereby the undrained strength gradually decreases to a certain state where the strain exceeds a maximum value; this strength value is referred to as the residual strength of the soil. Using this concept, it is possible to analyze the deformation of the soil-pile system caused by liquefaction through dynamic numerical analysis. This analysis is conducted by evaluating the zone of liquefaction caused by the seismic load, calculating the residual strength after liquefaction, and inputting the strength parameter corresponding to the residual strength into the zone of liquefaction. Liquefaction analysis by this process is referred to as total stress analysis.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 22 decreases the effective stress of the soil, and liquefaction occurs when the generated excess pore pressure equals the initial effective stress of the soil. In this case, soil completely loses its shear strength and exhibits behavior similar to that of water. Soon, it poses a serious threat to the stability of the structure's foundation. Therefore, it is crucial to accurately predict the structural behavior of foundations in the event of liquefaction to prevent damage to the structures. For this purpose, two types of analysis method are generally employed: total stress analysis and effective stress analysis.

Total Stress Analysis
In the event of an earthquake, the seismic wave propagated from the bedrock reaches upper soil particles through the soil layer, as shown in Figure 1. In this case, liquefaction occurs when the underground cyclic shear stress caused by the earthquake exceeds the magnitude of the cyclic shear stress required to cause liquefaction. The zone of liquefaction is determined by comparing the magnitudes of the generated cyclic shear stress and shear strength according to the depth, as shown in Figure 2 (Seed & Idriss, 1971 [18]). Liquefaction is similar to the strain-softening behavior that occurs when sandy soil is subjected to shear stress under undrained conditions. This is the concept whereby the undrained strength gradually decreases to a certain state where the strain exceeds a maximum value; this strength value is referred to as the residual strength of the soil. Using this concept, it is possible to analyze the deformation of the soil-pile system caused by liquefaction through dynamic numerical analysis. This analysis is conducted by evaluating the zone of liquefaction caused by the seismic load, calculating the residual strength after liquefaction, and inputting the strength parameter corresponding to the residual strength into the zone of liquefaction. Liquefaction analysis by this process is referred to as total stress analysis.

Effective Stress Analysis
Effective stress analysis can identify the liquefaction behavior at a desired position in real time by integrating the three steps of the total stress analysis into one. Although the equations and calculations for this analysis are much more complicated and require more time than that of the total stress analysis, highly reliable results can be expected when an appropriate model is used because this approach is more accurate and straightforward. Therefore, this method requires a liquefaction model that can accurately implement the liquefaction behavior of the target soil in the event of an earthquake. In this study, effective stress analysis was applied and the dynamic behavior of the soilpile-structure system during liquefaction was evaluated using the Finn model, a built-in liquefaction model in FLAC3D. The applicability of the Finn model was verified through comparisons with various field observations and model tests in several studies conducted by Soroush & Koohi (2004) [19], Liu et al. (2011) [20], Moradi et al. (2011) [21], and Ghorbani & Jahanpour (2012) [22]. Moreover, by incorporating the Finn model into Mohr-Coulomb criteria, it demonstrated good compatibility for different types of analysis such as dry (Kwon & Yoo, 2019 [23]) and liquefiable conditions.

Soil Liquefaction Model
If shear deformation is applied to the saturated soil element under undrained conditions, excess pore pressure is generated owing to slippage between the soil particles. This decreases the effective stress and simultaneously expands the volume of the soil. In each cyclic shear deformation step, the compatibility conditions for this volume change can be expressed as shown in Equation (1) (Martin et al., 1975 [24]). In other words, the change in the void ratio under undrained conditions (water compression amount) is equal to the value obtained by subtracting the volume expansion of the soil element caused by pore pressure from the volume change of the soil element under drained conditions.
where ∆u is the pore pressure increment, Δεvd is the volume reduction, Er is the tangential coefficient of the unloading curve, ne is the porosity, and Kw is the bulk modulus of water. The volume strain increment that induces a change in the pore pressure in Equation (1) can be expressed using Equation (2).

Zone of liquefaction
Cyclic shear stress required to cause liquefaction Equivalent cyclic shear stress induced by earthquake Shear stress Depth Figure 2. Determination of zone of liquefaction (Seed & Idriss, 1971 [18]).

Effective Stress Analysis
Effective stress analysis can identify the liquefaction behavior at a desired position in real time by integrating the three steps of the total stress analysis into one. Although the equations and calculations for this analysis are much more complicated and require more time than that of the total stress analysis, highly reliable results can be expected when an appropriate model is used because this approach is more accurate and straightforward. Therefore, this method requires a liquefaction model that can accurately implement the liquefaction behavior of the target soil in the event of an earthquake. In this study, effective stress analysis was applied and the dynamic behavior of the soil-pile-structure system during liquefaction was evaluated using the Finn model, a built-in liquefaction model in FLAC3D. The applicability of the Finn model was verified through comparisons with various field observations and model tests in several studies conducted by Soroush & Koohi (2004) [19], Liu et al. (2011) [20], Moradi et al. (2011) [21], and Ghorbani & Jahanpour (2012) [22]. Moreover, by incorporating the Finn model into Mohr-Coulomb criteria, it demonstrated good compatibility for different types of analysis such as dry (Kwon & Yoo, 2019 [23]) and liquefiable conditions.

Soil Liquefaction Model
If shear deformation is applied to the saturated soil element under undrained conditions, excess pore pressure is generated owing to slippage between the soil particles. This decreases the effective stress and simultaneously expands the volume of the soil. In each cyclic shear deformation step, the compatibility conditions for this volume change can be expressed as shown in Equation (1) (Martin et al., 1975 [24]). In other words, the change in the void ratio under undrained conditions (water compression amount) is equal to the value obtained by subtracting the volume expansion of the soil element caused by pore pressure from the volume change of the soil element under drained conditions.
(∆u · n e )/K w = ∆ε vd -(∆u/E r ) (1) where ∆u is the pore pressure increment, ∆ε vd is the volume reduction, E r is the tangential coefficient of the unloading curve, n e is the porosity, and K w is the bulk modulus of water. The volume strain increment that induces a change in the pore pressure in Equation (1) can be expressed using Equation (2).
where ε vd is the cumulative volume strain up to N cycles, and γ is the shear strain action in cycle N + 1.
In Equation (2), C 1 , C 2 , C 3 , and C 4 are the volume change constants, which vary depending on soil conditions. The volume change constants, which are important variables for liquefaction analysis, can be calculated through indoor experiments. For practical purposes, Byrne (1991) [25] proposed a simplified equation (Equation (3)) that can replace Equation (2). Equation (3) represents an empirical method. Byrne (1991) [25] also proposed a calculation equation (Equation (4)) for the volume change constants C 1 and C 2 using the relative density of the soil.  [27]). By incorporating both constitutive model, basic failure criteria and a related physical mechanism was modeled by Mohr-Coulomb criteria, and additional complex behaviors due to liquefaction, such as dynamic pore pressure generation, reduction of effective stress, etc. were captured by the Finn liquefaction model.

Soil Dynamic Damping and Nonlinear Model
Soil nonlinearity significantly influences the dynamic behavior of a pile and superstructure. If a simple elastic or elasto-plastic model as Mohr-Coulomb criteria is employed, implementation of an additional damping is necessary to properly represent the hysteretic characteristics (Han & Hart (2006) [28], Miglio et al. (2008) [29], Kim et al. (2012) [30], Cheng et al. (2013) [31], Kwon & Yoo (2019) [23]). In this study, the hysteretic damping model was applied to consider both the nonlinearly varying soil modulus and energy dissipation. The hysteretic damping model expresses the tangential elastic modulus of the G/G max -γ curve as a function of shear strain, as shown in Equation (5) (Itasca Consulting Group, 2006 [17]). L 1 and L 2 are coefficients that determine the rate and onset of the decrease in the G/G max value in the G/G max -γ curve. In this study, the values of L 1 and L 2 that simulated the experimental results most accurately were determined to be −3.65 and 0.5 by comparing the G/G max -γ curve of Nevada sand derived through indoor experiments (Seed & Idriss, 1970 [32] and Papadimitriou et al., 2001 [33]) with those derived from the changes in L 1 and L 2 through iterations (Kwon et al., 2016 [34], Figure 3).
where M t is the tangential elastic modulus, s = (L 2 − L)/(L 2 -L 1 ), L = log 10 γ, L 1 and L 2 are arbitrary constants, γ is the shear strain, and e is the void ratio. Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 22 The initial shear modulus of the soil depends on the magnitude of the confining pressure according to the soil depth. In this study, the initial shear modulus was calculated using Equation (6) proposed by Hardin & Drnevich (1972) [35].
where F(e) = 1/(0.3 + 0.7e 2 ), e is the void ratio, k is the overconsolidation ratio index (0 for sand), σ'm is the average principal stress (σ'm = (σ'1 + σ'2 + σ'3)/3, Pa is the atmospheric pressure; empirical constants A and n are the empirical coefficients determined to be 247.73 and 0.567, respectively, through bender element tests under various confining pressures (Yang, 2009 [36]). The distribution of the initial shear elastic modulus based on depth was input by applying the above relationship to every element using the FISH function, which is a built-in function of FLAC3D.

Soil-Pile Foundation Interface Model
In the event of a strong earthquake, slippage and separation phenomena may occur at the contact between the pile foundation and soil according to the magnitude of the applied load. As such phenomena directly affect the seismic behavior of the pile foundation, it is essential to utilize appropriate soil-pile interface elements for accurate modeling. Therefore, in this study, a soil-pile interface model that considered full contact, slippage, and separation between the pile and soil was applied. The interface model calculated the stiffness of the soil-pile boundary elements using vertical and shear springs at each position, and the value was determined through Equation (7) (Itasca Consulting Group, 2006 [17]). As the shear elastic modulus (G) introduced in Equation (6) considers the nonlinearity of soil obtained from the hysteretic damping model, the applied boundary element model also simulates the nonlinear behavior of soil; it was continuously applied and entered in the vertical and shear directions based on depth through the FISH function.
where KE is the equivalent stiffness of the interface spring, K is the bulk modulus, G is the shear elastic modulus, and Δzmin is the smallest width of the soil element adjacent to the pile in the vertical direction. Further, the yield criterion for the pile-soil boundary was determined using Equation (8), and the three boundary conditions between the pile and soil (full contact, slippage, and separation) were considered depending on the magnitude of the external force.  The initial shear modulus of the soil depends on the magnitude of the confining pressure according to the soil depth. In this study, the initial shear modulus was calculated using Equation (6) proposed by Hardin & Drnevich (1972) [35].
where F(e) = 1/(0.3 + 0.7e 2 ), e is the void ratio, k is the overconsolidation ratio index (0 for sand), σ' m is the average principal stress (σ' m = (σ' 1 + σ' 2 + σ' 3 )/3, P a is the atmospheric pressure; empirical constants A and n are the empirical coefficients determined to be 247.73 and 0.567, respectively, through bender element tests under various confining pressures (Yang, 2009 [36]). The distribution of the initial shear elastic modulus based on depth was input by applying the above relationship to every element using the FISH function, which is a built-in function of FLAC3D.

Soil-Pile Foundation Interface Model
In the event of a strong earthquake, slippage and separation phenomena may occur at the contact between the pile foundation and soil according to the magnitude of the applied load. As such phenomena directly affect the seismic behavior of the pile foundation, it is essential to utilize appropriate soil-pile interface elements for accurate modeling. Therefore, in this study, a soil-pile interface model that considered full contact, slippage, and separation between the pile and soil was applied. The interface model calculated the stiffness of the soil-pile boundary elements using vertical and shear springs at each position, and the value was determined through Equation (7) (Itasca Consulting Group, 2006 [17]). As the shear elastic modulus (G) introduced in Equation (6) considers the nonlinearity of soil obtained from the hysteretic damping model, the applied boundary element model also simulates the nonlinear behavior of soil; it was continuously applied and entered in the vertical and shear directions based on depth through the FISH function.
where K E is the equivalent stiffness of the interface spring, K is the bulk modulus, G is the shear elastic modulus, and ∆z min is the smallest width of the soil element adjacent to the pile in the vertical direction. Further, the yield criterion for the pile-soil boundary was determined using Equation (8), and the three boundary conditions between the pile and soil (full contact, slippage, and separation) were considered depending on the magnitude of the external force. where c is the adhesion at the boundary, Φ is the friction angle at the boundary, p is the pore pressure, and A is the representative area of the boundary. For the internal friction angle at the soil-pile interface, it is generally appropriate to input a value slightly less than the maximum value of far-field soil internal friction angle because dynamic behavior of the pile foundation doesn't have effect on far-field soil. Kraft (1990) [37] suggested 70% of far-filed value as the interface value, in addition, Reddy et al. (2000) [38] suggested 60% of far-field value for the interface. In this study, calculation basis of the internal friction angle at the soil-pile interface proposed by Beringen et al. (1979) [39] and Randolph et al. (1994) [40] was applied as Equation (9).
where δ is the internal friction angle of the interface, and Φ max is the maximum internal friction angle of the far-field soil.

Boundary Conditions
Three-dimensional numerical modeling is essential for appropriate and accurate simulation of soil-pile-structure dynamic interaction; however, 3D full-scale analysis is often time consuming, particularly when applied to dynamic problems. If an infinite number of soil elements are created to simulate the semi-infinite boundary of the actual site, the numerical analysis of various conditions becomes difficult. Further, the analysis efficiency decreases significantly because the analysis time increases sharply. However, if boundary conditions are not properly set in dynamic numerical modeling, the input wave propagating from the center of the model can generate a reflection wave at the boundary, complicating the simulation of the actual behavior. Therefore, in this study, the simplified continuum modeling method shown in Figure 4 was applied to reduce the analysis time and to properly the simulate dynamic behavior at the model boundary. Figure 4 shows the schematic of the simplified continuum modeling method adopted in the present study. The entire region is divided into two parts: near-field and far-field zones. The soil motion in the near-field zone is influenced by the soil-pile interaction and modeled in finite zones. The soil motion in the far-field zone not influenced by the soil-pile interaction is replaced as the boundary motion of the near-field zone. This modeling method can reduce analysis time compared to full modeling approach of two zones. Nogami [46] applied similar approaches. However, their approaches replaced the far-field zone with nonlinear spring elements that simulate the interaction stress between the pile and nearby soil. In addition, these approaches were limited to 1D or 2D modeling. Kim (2011) [47], Kim et al. (2012) [30], and Kwon et al. (2013) [48] applied the simplified method for 3D modeling to model a boundary condition of the soil-pile system in dry sand and validate applicability of the boundary model which input proper level of acceleration-time histories at the horizontal boundary of the model. In the present study, the method has been extended to simulate the dynamic behavior of a pile in liquefiable sand. The soil motion in the far-field zone is obtained from the independent site response analysis, which can simulate the liquefaction behavior of sandy soil. The horizontal soil acceleration-time histories along the depth in the far-field zone are applied at the boundary nodes of the near-field zone.
In the simplified continuum modeling method, it is necessary to determine the boundary extent of the near-field zone where the soil motion becomes identical to the far-field movement. In order to determine the optimal boundary extent of the near-field, the acceleration amplification ratio, which is the ratio of the amplitude of the maximum acceleration at the surface to the base input, was analyzed. The outline procedure of simplified continuum modeling is as follows. (1) Target system is modeled and analyzed by exerting earthquake load. (2) Obtain extent of near-field elements by plotting amplification ratio according to the distance from the pile center. (3) Only near-field is modeled and acceleration-time histories in far-field for each depth are input as boundary condition. Figure 5 shows the amplification ratio according to the distance from the center of the pile. The dynamic excitation of the pile increases the surface acceleration. Therefore, as the distance from the pile increases, the ratio decreases. The optimal distance with a constant ratio was determined to be 8 D (where D is the pile diameter). Finally, the boundary extent of the near-field was applied as 10 D from the center of the pile for the conservative modeling.
Appl. Sci. 2020, 10, x FOR PEER REVIEW  8 of 22 pile increases, the ratio decreases. The optimal distance with a constant ratio was determined to be 8 D (where D is the pile diameter). Finally, the boundary extent of the near-field was applied as 10 D from the center of the pile for the conservative modeling.   pile increases, the ratio decreases. The optimal distance with a constant ratio was determined to be 8 D (where D is the pile diameter). Finally, the boundary extent of the near-field was applied as 10 D from the center of the pile for the conservative modeling.  2.6. Numerical Modeling Figure 6 shows typical mesh and boundary extents adopted in this study. By applying the simplified modeling, the number of elements was reduced from approximately 6000 to 4000 and the analysis time was remarkably reduced by about 50% in comparison with that of the full modeling of near-and far-field zones. By applying symmetric conditions, half of the entire system was modeled. The vertical and horizontal displacements at the bottom boundary and the normal displacements at the lateral boundaries were constrained. The size of the soil elements gradually increased from the pile foundation to the domain boundary. In addition, the optimum mesh sizes were determined to minimize the effect of the mesh size on the results from the preliminary analyses. The 3D elements of the hexahedron type were selected to model the soil domain and the pile foundation. As a result of modeling the pile by applying elastic model, it was found that the stress on pile foundation induced in the analysis process was even smaller than the yielding stress of material itself. The weight of the superstructure was correctly simulated by controlling the unit weight of the corresponding elements. Numerical analyses were performed following three steps; (1) simulation of geostatic condition, (2) pile installation, and (3) dynamic loading. At the first step, an in-situ geostatic condition was simulated by achieving the stress equilibrium under gravitational load. In the second step, the soil elements for pile area were eliminated and pile elements were created with the interface elements. The stress equilibrium condition was then gained for increasing stress due to pile weight once again. Finally, dynamic loading was applied at the bottom boundary and the left and right boundaries corresponding to loading direction. The time history of input seismic acceleration was applied to the bottom boundary nodes, and time histories of far-filed horizontal acceleration were applied to the left and right boundary nodes.

Numerical Modeling
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 22 Figure 6 shows typical mesh and boundary extents adopted in this study. By applying the simplified modeling, the number of elements was reduced from approximately 6000 to 4000 and the analysis time was remarkably reduced by about 50% in comparison with that of the full modeling of near-and far-field zones. By applying symmetric conditions, half of the entire system was modeled. The vertical and horizontal displacements at the bottom boundary and the normal displacements at the lateral boundaries were constrained. The size of the soil elements gradually increased from the pile foundation to the domain boundary. In addition, the optimum mesh sizes were determined to minimize the effect of the mesh size on the results from the preliminary analyses. The 3D elements of the hexahedron type were selected to model the soil domain and the pile foundation. As a result of modeling the pile by applying elastic model, it was found that the stress on pile foundation induced in the analysis process was even smaller than the yielding stress of material itself. The weight of the superstructure was correctly simulated by controlling the unit weight of the corresponding elements. Numerical analyses were performed following three steps; (1) simulation of geostatic condition, (2) pile installation, and (3) dynamic loading. At the first step, an in-situ geostatic condition was simulated by achieving the stress equilibrium under gravitational load. In the second step, the soil elements for pile area were eliminated and pile elements were created with the interface elements. The stress equilibrium condition was then gained for increasing stress due to pile weight once again. Finally, dynamic loading was applied at the bottom boundary and the left and right boundaries corresponding to loading direction. The time history of input seismic acceleration was applied to the bottom boundary nodes, and time histories of far-filed horizontal acceleration were applied to the left and right boundary nodes.

Centrifuge Tests and Input Properties
Verification of the proposed numerical model was carried out using the results from the centrifuge test conducted by Wilson (1998) [16]. The centrifuge test was conducted under centrifugal acceleration of 30 g, which corresponds to the length scaling factor of 30 using the National Geotechnical Centrifuge tester from UC DAVIS, and the model soil was constructed with two layers

Centrifuge Tests and Input Properties
Verification of the proposed numerical model was carried out using the results from the centrifuge test conducted by Wilson (1998) [16]. The centrifuge test was conducted under centrifugal acceleration of 30 g, which corresponds to the length scaling factor of 30 using the National Geotechnical Centrifuge tester from UC DAVIS, and the model soil was constructed with two layers using uniform Nevada sand. The dense bottom layer (D r = 80%) was 11.4 m thick and the loose top layer (D r = 35% or 55%) was 9.1 m thick. All the soil layers were completely saturated (Figure 7).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 22 using uniform Nevada sand. The dense bottom layer (Dr = 80%) was 11.4 m thick and the loose top layer (Dr = 35% or 55%) was 9.1 m thick. All the soil layers were completely saturated (Figure 7).
The model pile was fabricated with aluminum. It measured 0.67 m in diameter and was 0.072 m thick. The pile head was located 3.8 m from the ground surface, and a superstructure load of 480 kN was installed. The penetration depth of the pile was 16.8 m, and the pile end simulated the friction pile located approximately 3.7 m from the bottom of the soil container. The soil container was filled with a fluid of a viscosity approximately ten times greater than water, which was fabricated by mixing a viscous fluid (Hydroxyl-propylmethyl-cellulose) with water to generate saturated soil conditions and to satisfy similitude law for the permeability (Stewart et al., 1998 [49]). For the input seismic wave, the Kobe (1995) seismic wave (Port Island 83-m depth, NS direction) was used, for which the maximum acceleration value was scaled to 0.22 g. Figure 8 shows the input accelerationtime history used in experiments and numerical analysis.
For the numerical simulation of the experimental model, the modeling method described in Section 2 was applied. The soil properties applied in the numerical model are summarized in Table  1 (Popescu & Prevost, 1993 [50]). Variation of internal friction angle and Poisson's ratio according to the relative density was determined by Beringen et al. (1979) [39] and Kumar & Madhusudhan (2010) [51], respectively. The pile was modeled as a cylinder and had an equivalent flexible stiffness with the model pile of the centrifuge test. The pile properties are summarized in Table 2.   (Wilson, 1998 [16]).
The model pile was fabricated with aluminum. It measured 0.67 m in diameter and was 0.072 m thick. The pile head was located 3.8 m from the ground surface, and a superstructure load of 480 kN was installed. The penetration depth of the pile was 16.8 m, and the pile end simulated the friction pile located approximately 3.7 m from the bottom of the soil container. The soil container was filled with a fluid of a viscosity approximately ten times greater than water, which was fabricated by mixing a viscous fluid (Hydroxyl-propylmethyl-cellulose) with water to generate saturated soil conditions and to satisfy similitude law for the permeability (Stewart et al., 1998 [49]). For the input seismic wave, the Kobe (1995) seismic wave (Port Island 83-m depth, NS direction) was used, for which the maximum acceleration value was scaled to 0.22 g. Figure 8 shows the input acceleration-time history used in experiments and numerical analysis.
using uniform Nevada sand. The dense bottom layer (Dr = 80%) was 11.4 m thick and the loose top layer (Dr = 35% or 55%) was 9.1 m thick. All the soil layers were completely saturated (Figure 7). The model pile was fabricated with aluminum. It measured 0.67 m in diameter and was 0.072 m thick. The pile head was located 3.8 m from the ground surface, and a superstructure load of 480 kN was installed. The penetration depth of the pile was 16.8 m, and the pile end simulated the friction pile located approximately 3.7 m from the bottom of the soil container. The soil container was filled with a fluid of a viscosity approximately ten times greater than water, which was fabricated by mixing a viscous fluid (Hydroxyl-propylmethyl-cellulose) with water to generate saturated soil conditions and to satisfy similitude law for the permeability (Stewart et al., 1998 [49]). For the input seismic wave, the Kobe (1995) seismic wave (Port Island 83-m depth, NS direction) was used, for which the maximum acceleration value was scaled to 0.22 g. Figure 8 shows the input accelerationtime history used in experiments and numerical analysis.
For the numerical simulation of the experimental model, the modeling method described in Section 2 was applied. The soil properties applied in the numerical model are summarized in Table  1 (Popescu & Prevost, 1993 [50]). Variation of internal friction angle and Poisson's ratio according to the relative density was determined by Beringen et al. (1979) [39] and Kumar & Madhusudhan (2010) [51], respectively. The pile was modeled as a cylinder and had an equivalent flexible stiffness with the model pile of the centrifuge test. The pile properties are summarized in Table 2.   For the numerical simulation of the experimental model, the modeling method described in Section 2 was applied. The soil properties applied in the numerical model are summarized in Table 1 (Popescu & Prevost, 1993 [50]). Variation of internal friction angle and Poisson's ratio according to the relative density was determined by Beringen et al. (1979) [39] and Kumar & Madhusudhan (2010) [51], respectively. The pile was modeled as a cylinder and had an equivalent flexible stiffness with the model pile of the centrifuge test. The pile properties are summarized in Table 2.

Comparison of Dynamic Responses
The applicability of the proposed modeling method was verified by comparing various dynamic responses under earthquake loading using two experimental cases. Figure 9 shows time histories of the excess pore water pressure according to the relative density of upper soil layer and depth. In this figure, the excess pore pressure ratio is calculated by dividing the excess pore water pressure by the initial vertical effective stress in soil. The peak value of the excess pore pressure ratio in Case 1 (D r = 55%) reached 1.0 (means liquefaction) at 4.5 m and decreased to 0.5 at 7 m because the soil at such depths had higher confining pressure. However, the peak ratio in Case 2 (D r = 35%) was 1.0 in the entire region of the upper layer. It appears that liquefaction also occurs at greater depth in Case 2, unlike that in Case 1, because the relative density of the upper soil is less than that of Case 1; thus, the confining pressure is relatively lower. The time histories of the ratios computed from the proposed model agreed well with the experimental results in both cases. Thus, it is demonstrated that the proposed model can properly simulate the excess pore pressure response of the liquefiable sand with various relative densities.   Figure 10 shows time histories of the pile bending moment according to the depth and the soil relative density. It shows that the amplitude of the bending moment decreased after 3 s when the excess pore pressure reached the peak values. In addition, the peak values of Case 2 (D r = 35%) were smaller than those of Case 1 (D r = 55%). The proposed model can accurately predict the peak amplitude and the decreasing trend of the bending moment during liquefaction. Figure 11 shows time history of the horizontal displacement at the pile head. As shown in the figure, the displacement calculated from the proposed model accurately predicted the displacement amplitude and the decrease of the amplitude during liquefaction. Despite the slight discrepancy in phase, the calculated maximum moment and displacement reasonably agreed with the measured values of the experiments with discrepancy less than 10%. Finally, it can be concluded that the proposed model can simulate the dynamic soil-pile-structure interaction in liquefiable sand under various conditions. In the verification results, on the other hand, the maximum bending moment at 1 m depth near the ground surface and the maximum pile displacement occurs at approximately 3.5 s. This almost coincides with the time when the excess pore pressure ratio sharply increases in Figure 9. This appears to be because of the influence of kinematic force caused by softening and deformation of the soil around the pile as well as the inertial force generated by the pile head.

Parametric Study
A parametric study was performed using the proposed model to investigate various dynamic behavior of a pile embedded in liquefiable sand. The selected parameters were thickness of liquefiable soil layer, soil relative density, and weight of a superstructure ( Table 3). The input earthquake motion at the base was the sinusoidal wave with amplitude of 0.13 g or 0.25 g, which are similar to the level of design ground acceleration in South Korea. Sinusoidal waves were adopted as the input earthquake motion for parametric study to obtain the dynamic performance of the soil-pile system under various seismic input conditions. The modeling of the pile and the soil was basically identical with the verification study. The relative density of the upper liquefiable layer was assumed to be 30% or 50%, and the soil property corresponding to each relative density was determined by applying the suggested methods, as shown in Table 4. The basic input values of the parameters were as follows: input frequency of 1 Hz, upper liquefiable layer of 9-m thickness and D r = 30%, and superstructure weight of 300 kN. The value of each target parameter was changed after fixing the values of the other parameters. Table 3. Selected parameters in a parametric study.

Effect of Thickness of Liquefiable Soil Layer
Dynamic pile responses obtained from repetitive analysis for varying thickness of liquefiable soil layer (i.e., 3, 6, and 9 m) were discussed. Figure 12 shows the peak envelopes of the bending moment and the pile lateral displacement according to the thickness of the liquefiable layer when input acceleration of 0.13 g is applied. Pile foundation with superstructure is usually exerted by two kinds of dynamic forces representatively under earthquake. One is inertial force induced by superstructure, the other is kinematic force induced by ground movement. These forces induce complex dynamic interactions to the pile, and combination mechanism of the two forces acting on the pile depends on the site condition. The way of combination mechanism can govern overall performance of the pile foundation under earthquake. Tokimatsu et al. (2004) [52] reported that, if the natural period of the superstructure is less than that of the ground, the ground displacement tends to be in phase with the inertial force from the superstructure. In contrast, if natural period of the superstructure is greater than that of the ground, the ground displacement tends to be out of phase with the inertial force, restraining the pile stress from increasing. In Figure 12, when the thickness of liquefiable layer is increased, larger soil deformation is induced. The maximum displacement decreased by approximately 34% as the thickness of liquefiable layer decreased from 9 m to 3 m. Peak bending moment value also increased as thickness of liquefiable layer increased. If soil layer is liquefied, natural period of the ground increased because soil stiffness decreased (natural period and stiffness are inverse proportion). Therefore, the natural period of the liquefied soil is always greater than that of the superstructure, resulting in kinematic force due to ground movement in phase with the inertial force in this case. In this situation, as the thickness of liquefiable layer increased, larger kinematic force is exerted with inertial force to the pile structure. On the other hand, the bending moment profile reached the maximum at the interface between the liquefied and non-liquefied layers regardless of analysis case. The pile with the upper liquefied sand behaved like a cantilever beam embedded below the non-liquefied layer. Finn & Fujita (2002) [11] reported that many pile damages were observed at the interface between the liquefied and non-liquefied layers from the investigations of past earthquakes. It can be demonstrated that this kind of destruction pattern was caused by sharp increase of bending moment at the interface and captured by proposed numerical model. Also, it is confirmed that critical seismic design section with the bending moment is the interface between liquefied and non-liquefied regions by numerical simulation.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 22 than that of the superstructure, resulting in kinematic force due to ground movement in phase with the inertial force in this case. In this situation, as the thickness of liquefiable layer increased, larger kinematic force is exerted with inertial force to the pile structure. On the other hand, the bending moment profile reached the maximum at the interface between the liquefied and non-liquefied layers regardless of analysis case. The pile with the upper liquefied sand behaved like a cantilever beam embedded below the non-liquefied layer. Finn & Fujita (2002) [11] reported that many pile damages were observed at the interface between the liquefied and non-liquefied layers from the investigations of past earthquakes. It can be demonstrated that this kind of destruction pattern was caused by sharp increase of bending moment at the interface and captured by proposed numerical model. Also, it is confirmed that critical seismic design section with the bending moment is the interface between liquefied and non-liquefied regions by numerical simulation.  Figure 13 shows the peak envelopes of the bending moment and the pile lateral displacement along the depth both for Dr = 30% and Dr = 50% when input acceleration of 0.13 g is applied. As the relative density was increased from 30% to 50%, the peak value of the bending moment and the pile lateral displacement was decreased by approximately 45% and 40%, respectively. It would be due to the increases of liquefiable soil stiffness and relatively smaller values of kinematic force induced by soil movement. Usually, liquefaction hardly occurs in deeper depth because of the high confining pressure. The fully liquefied depth was approximately 3 m at Dr = 50% and approximately 5 m at Dr = 30%. Therefore, as the liquefiable layer became denser, the liquefied depth became shallower and the bending moment and the pile lateral displacement significantly decreased. These results show that soil compaction would be very effective for stabilization of a pile in liquefiable soil. On the other hand, this kind of behavior is wholly different from the dry sand case. Figure 14 shows peak  Figure 13 shows the peak envelopes of the bending moment and the pile lateral displacement along the depth both for D r = 30% and D r = 50% when input acceleration of 0.13 g is applied. As the relative density was increased from 30% to 50%, the peak value of the bending moment and the pile lateral displacement was decreased by approximately 45% and 40%, respectively. It would be due to the increases of liquefiable soil stiffness and relatively smaller values of kinematic force induced by soil movement. Usually, liquefaction hardly occurs in deeper depth because of the high confining pressure. The fully liquefied depth was approximately 3 m at D r = 50% and approximately 5 m at D r = 30%. Therefore, as the liquefiable layer became denser, the liquefied depth became shallower and the bending moment and the pile lateral displacement significantly decreased. These results show that soil compaction would be very effective for stabilization of a pile in liquefiable soil. On the other hand, this kind of behavior is wholly different from the dry sand case. Figure 14 shows peak envelopes of the bending moment and pile lateral displacement for different relative densities in dry sand (Kwon & Yoo, 2019 [23]). Strict quantitative comparison between two analysis cases would be improper because analysis condition of dry sand case was slightly different from liquefiable case, but qualitative comparison would be worthwhile to investigate critical distinction in trend of dynamic pile behavior according to the soil condition. It can be seen that dynamic pile responses show almost identical values for different relative densities in dry sand. It is significantly different pattern with liquefiable condition. In Figure 15, this phenomenon can be clearly identified by comparing peak values of dynamic bending moment obtained in different soil conditions. In the figure, black line graph with grey dots shows percentage of gap between peak value of dynamic bending moment for D r = 30% case and D r = 50% case for different soil conditions. The gap in liquefied soil condition shows much larger value than the gap in dry soil condition. It can be demonstrated that kinematic force is much more dominant in liquefiable sand than in dry sand from the results. Moreover, peak value of dynamic bending moment in dry soil condition for D r = 50% was larger than that for D r = 30%, which is different pattern with liquefied soil condition. Ground displacement tends to be out of phase with the inertial force, restraining the pile stress from increasing in dry soil condition, so it can be identified that D r = 50% case exerted larger reaction forces on the pile foundation than D r = 30% case in the figure. Based on a series of analyses, important dynamic characteristics as kinematic and inertial interaction for different soil condition were clarified by full-scale numerical simulation.

Effect of Relative Density of Liquefiable Soil
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 22 envelopes of the bending moment and pile lateral displacement for different relative densities in dry sand (Kwon & Yoo, 2019 [23]). Strict quantitative comparison between two analysis cases would be improper because analysis condition of dry sand case was slightly different from liquefiable case, but qualitative comparison would be worthwhile to investigate critical distinction in trend of dynamic pile behavior according to the soil condition. It can be seen that dynamic pile responses show almost identical values for different relative densities in dry sand. It is significantly different pattern with liquefiable condition. In Figure 15, this phenomenon can be clearly identified by comparing peak values of dynamic bending moment obtained in different soil conditions. In the figure, black line graph with grey dots shows percentage of gap between peak value of dynamic bending moment for Dr = 30% case and Dr = 50% case for different soil conditions. The gap in liquefied soil condition shows much larger value than the gap in dry soil condition. It can be demonstrated that kinematic force is much more dominant in liquefiable sand than in dry sand from the results. Moreover, peak value of dynamic bending moment in dry soil condition for Dr = 50% was larger than that for Dr = 30%, which is different pattern with liquefied soil condition. Ground displacement tends to be out of phase with the inertial force, restraining the pile stress from increasing in dry soil condition, so it can be identified that Dr = 50% case exerted larger reaction forces on the pile foundation than Dr = 30% case in the figure. Based on a series of analyses, important dynamic characteristics as kinematic and inertial interaction for different soil condition were clarified by full-scale numerical simulation.

Effect of Superstructure Weight
The effect of the weight of superstructure on dynamic pile response was analyzed by varying the weight of the superstructure at 0, 300, 600, and 900 kN. Figure 16 shows peak envelopes of the pile lateral displacement according to the weight of a superstructure when input frequency of 0.5 Hz was applied. Although the superstructure weight increased from 0 to 900 kN, the lateral pile displacements were very similar regardless of the magnitude of input acceleration. The maximum

Effect of Superstructure Weight
The effect of the weight of superstructure on dynamic pile response was analyzed by varying the weight of the superstructure at 0, 300, 600, and 900 kN. Figure 16 shows peak envelopes of the pile lateral displacement according to the weight of a superstructure when input frequency of 0.5 Hz was applied. Although the superstructure weight increased from 0 to 900 kN, the lateral pile displacements were very similar regardless of the magnitude of input acceleration. The maximum

Effect of Superstructure Weight
The effect of the weight of superstructure on dynamic pile response was analyzed by varying the weight of the superstructure at 0, 300, 600, and 900 kN. Figure 16 shows peak envelopes of the pile lateral displacement according to the weight of a superstructure when input frequency of 0.5 Hz was applied. Although the superstructure weight increased from 0 to 900 kN, the lateral pile displacements were very similar regardless of the magnitude of input acceleration. The maximum value of the pile displacement decreased only by 26% when the superstructure weight decreased from 900 kN to 0 kN under liquefiable conditions. Similar results were reported by Han (2006) [1]. Han (2006) [1] analyzed the dynamic earth pressure on piles in liquefiable soil by 1 g shaking table tests and found that the amplitude of the dynamic earth pressure acting on the pile and the pile responses are similar irrespective of the superstructure weight. Figure 17 shows maximum pile lateral displacement envelopes for four different weight of the superstructure embedded in dry sand (Kwon & Yoo, 2019 [23]). It was noted that lateral pile displacements were remarkably varied for different weights of the superstructure in dry condition, in contrast to liquefiable condition. It was identified that the maximum pile displacement in dry condition decreased to 92% when the superstructure weight decreased form 900 kN to 0 kN, whereas those in liquefiable condition decreased only 26%. Figure 18 shows the effect of inertial force in different soil condition under earthquake by comparing peak values of lateral pile displacement. The gap in liquefied soil condition shows much smaller value than the gap in dry soil condition. The difference between results from different types of soil condition was strikingly distinctive. It can be demonstrated that inertial force is much more dominant in dry soil conditions than in liquefiable soil conditions. The similar observation was reported through field investigations by Ishihara (1997) [53]. He concluded that the inertial force is the predominant force in dry condition and responsible for the maximum bending moment near the pile head. These results showed that the effect of inertia induced by the superstructure was not remarkable under liquefiable condition compared with that under dry conditions. Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 22 value of the pile displacement decreased only by 26% when the superstructure weight decreased from 900 kN to 0 kN under liquefiable conditions. Similar results were reported by Han (2006) [1]. Han (2006) [1] analyzed the dynamic earth pressure on piles in liquefiable soil by 1 g shaking table tests and found that the amplitude of the dynamic earth pressure acting on the pile and the pile responses are similar irrespective of the superstructure weight. Figure 17 shows maximum pile lateral displacement envelopes for four different weight of the superstructure embedded in dry sand (Kwon & Yoo, 2019 [23]). It was noted that lateral pile displacements were remarkably varied for different weights of the superstructure in dry condition, in contrast to liquefiable condition. It was identified that the maximum pile displacement in dry condition decreased to 92% when the superstructure weight decreased form 900 kN to 0 kN, whereas those in liquefiable condition decreased only 26%. Figure 18 shows the effect of inertial force in different soil condition under earthquake by comparing peak values of lateral pile displacement. The gap in liquefied soil condition shows much smaller value than the gap in dry soil condition. The difference between results from different types of soil condition was strikingly distinctive. It can be demonstrated that inertial force is much more dominant in dry soil conditions than in liquefiable soil conditions. The similar observation was reported through field investigations by Ishihara (1997) [53]. He concluded that the inertial force is the predominant force in dry condition and responsible for the maximum bending moment near the pile head. These results showed that the effect of inertia induced by the superstructure was not remarkable under liquefiable condition compared with that under dry conditions.

Concluding Remark
A numerical modeling methodology was proposed to analyze the dynamic soil-pile-structure interactive behavior embedded in liquefiable sand. The applicability of the proposed method was verified by comparing with the dynamic responses from the centrifuge tests performed by Wilson

Concluding Remark
A numerical modeling methodology was proposed to analyze the dynamic soil-pile-structure interactive behavior embedded in liquefiable sand. The applicability of the proposed method was verified by comparing with the dynamic responses from the centrifuge tests performed by Wilson

Concluding Remark
A numerical modeling methodology was proposed to analyze the dynamic soil-pile-structure interactive behavior embedded in liquefiable sand. The applicability of the proposed method was verified by comparing with the dynamic responses from the centrifuge tests performed by Wilson (1998) [16]. A parametric study was conducted to provide better discernment into the dynamic behavior of pile foundation under liquefaction. Detailed conclusions are as follows: (1) The proposed methodology includes the liquefaction model, which can simulate pore pressure generation and soil strength reduction under strong earthquake in real time. Simplified continuum modeling was adopted firstly for liquefaction related soil-pile-structure system to accurately simulate boundary condition and to improve analysis efficiency. Damping and interface models were proposed to properly consider soil nonlinearity, and various dynamic soil properties such as initial shear modulus, interface friction angle, etc. were carefully determined. (2) Based on the verification, it was confirmed that time histories of the excess pore pressure ratio, pile bending moment, and pile head lateral displacement by depth derived using the proposed modeling method effectively simulated those observed in centrifuge test. (3) In a parametric study, the bending moment reached its maximum at the interface between the liquefied and non-liquefied layers, and the lateral pile displacement remarkably decreased as the thickness of liquefiable layer was decreased. As the liquefiable layer became denser, the liquefied depth became shallower and dynamic pile responses significantly decreased. The kinematic force induced by soil deformation was dominant in liquefiable sand, whereas the inertial force induced by superstructure was relatively insignificant. These kinds of important dynamic characteristics for different soil conditions were clarified by full-scale numerical simulation. (4) Based on the results from a series of parametric studies using the proposed numerical model, several recommendations can be drawn in design stage of pile foundation. Additional careful consideration should be made when constructing piles in multi-layered ground condition because unexpected unusual sharp increase of dynamic response can be occurred at the interface between soil layers. Even under similarly liquefied conditions, significant differences in the dynamic responses of pile can occur due to differences in ground relative density because kinematic force is dominant in liquefied condition. Therefore, when the liquefiable layer is distributed to a deep depth, a more rational design is possible if an effective soil improvement method such as dynamic compaction is used in combination with a pile, rather than massive reinforcement or the addition of a pile.
Author Contributions: S.Y.K. organized the paperwork, made a analysis plan, performed numerical analysis and suggested the numerical analysis method; M.Y. supported the numerical model verification by comparing test results helped the data analysis; all authors contributed to the writing of paper. All authors have read and agreed to the published version of the manuscript.