A Hybrid Approach Calculating Lateral Spreading Induced by Seismic Liquefaction

Liquefaction-induced lateral spreading has caused severe damages to the infrastructures. To predict the liquefaction-induced lateral spreading, a hybrid approach was proposed based on the Newmark sliding-block model. One-dimensional effective stress analysis based on the borehole investigation of the site was conducted to obtain the triggering time of liquefaction and acceleration time history. Shear wave velocity of the liquefiable soil was used to estimate the residual shear strength of liquefiable soil. +e limit equilibrium analysis was conducted to determine the yield acceleration corresponding with the residual shear strength of liquefied soil. +e liquefaction-induced lateral spreading was calculated based on the Newmark sliding-block model. A case study based on Wildlife Site Array during the 1987 Superstition Hills earthquake was conducted to evaluate the performance of the hybrid approach. +e results showed that the hybrid approach was capable of predicting liquefaction-induced lateral spreading and the calculated lateral spreading was 1.5 times the observed displacement in terms of Wildlife Site Array. Numerical simulations with two other constitutive models of liquefiable sand were conducted in terms of effective stress analyses to reproduce the change of lateral spreading and excess pore water ratio over the dynamic time of Wildlife Site Array. Results of numerical simulations indicated that the lateral spreading varied with the triggering time of liquefaction when different constitutive models were used. +e simulations using PM4sand and UBC3D-PLM constitutive models predicted 5.2 times and 4 times the observed lateral spreading, respectively. To obtain the site response, the motions recorded at and below the ground surface were analyzed using the Hilbert–Huang transform. +e low-frequency content of the motion below the ground surface was amplified at the ground surface, and the liquefaction effect resulted in a shift of the frequency content. By comparing the response spectra of the entire ground surface motion and the ground surface motion from the beginning to the triggering time of liquefaction, the liquefaction effect at the site was confirmed.


Introduction
One of the most significant topics of Geotechnical earthquake engineering is seismic liquefaction. Liquefaction-induced lateral spreading is the horizontal displacement of ground with gentle sloping underlain by liquefiable deposits. e sloping of the ground where lateral spreading occurs is usually less than 5% based on Bartlett and Youd [1], and the ground surface with water table is usually underlain by loose sand or silt deposits which will liquefy during an earthquake. Once liquefaction occurs due to the strong shaking, pore water pressure increases in the liquefied soil and the upper soil begins to move over the liquefied layer. Accompany with the lateral spreading, ground surface fissures, or cracks are observed.
Postearthquakes [2][3][4][5] have witnessed damages caused by lateral spreading. To predict the magnitude of lateral spread displacement and avoid unnecessary damages to the adjacent structures, foundations, and earth structures in the liquefied area, the development of a proper method evaluating liquefaction-induced lateral spreading is important. e prediction methods of lateral spreading mainly include empirical methods [6][7][8], numerical methods [9][10][11], Newmark sliding-block methods [12], and probabilistic methods [13][14][15]. Each method has some limitations to the practical use in engineering. For example, the empirical and probabilistic methods depend on the earthquake loading, site geometry, and physical property of the liquefiable soil, but the soil behavior due to dynamic loading is neglected. e numerical method depends on the advanced constitutive models of the liquefiable soils used in the seismic analysis. If the understanding of dynamic response of the soil is lacking, the analysis is difficult to be conducted.
In the engineering practice, a more practical method needs to be adopted for evaluating liquefaction-induced lateral spreading. e displacement of the overburden soil layer accumulates once the liquefaction is triggered and the liquefied soil retains the residual shear strength before the shaking ceases; meanwhile, the static factor of safety is greater than unit one after the shaking ceases. Based on these characteristics of liquefaction-induced lateral spreading, a Newmark-type displacement [16] can be used to describe the lateral spreading induced by liquefaction if assuming that the overburden soil is intact and slides over the liquefied soil during the shaking.
In the paper, a hybrid approach based on the Newmark sliding-block method and borehole investigation was proposed to evaluate the lateral spreading induced by liquefaction. e dynamic input is one of the most important factors for calculating Newmark displacement. Considering that the lateral spreading occurs after the liquefaction was triggered and the overburden soil slides over the sliding surface, the one-dimensional effective stress analysis is used in this paper to obtain the dynamic input motion beneath the liquefiable soil and triggering time of liquefaction. e yield acceleration, corresponding to the acceleration acting on the sliding mass when the seismic factor of safety is equal to 1.0, is the other important factor when conducting Newmark sliding-block analysis. Since the soil retains residual shear strength during the shaking, the yield acceleration based on the residual shear strength of the liquefied soil is the yield acceleration for calculating liquefaction-induced lateral spreading. e residual shear strength could be estimated from the SPT, CPT value, or shear wave velocity of the liquefiable soil. Several methods [17][18][19][20] were proposed to estimate the residual shear strength of liquefiable soil based on the correlation relationship developed from case histories of flow failure or lateral spreading. ough Olson and Johnson [17] stated, "lateral spreads back analyzed using the Newmark sliding-block procedure exhibit mobilized strength ratios essentially identical to liquefied strength ratios back calculated from flow failures." ere is no consensus that the residual shear strength of liquefiable soil for the lateral spreading case can be estimated using the residual shear strength models developed from cases of flow failure. Based on the research conducted by Pelin Ozener [20], the shear wave velocity and liquefaction of soil are both affected by the same factors. In the paper, the shear wave velocity of the liquefiable soil is used to estimate residual shear strength and reduce the uncertainty of the yield acceleration in the lateral spreading calculation.
To evaluate the hybrid approach, the Wildlife Site Array, widely known for successfully capturing the lateral spreading, pore water pressure generation, and ground surface response during the 1987 Superstition Hills earthquake, was reanalyzed in the paper through the proposed method.
e predicted lateral spreading of the 1987 Superstition Hills earthquake was compared with numerical simulations using two other constitutive models in finite difference analysis and finite-element analysis, respectively. e triggering time of liquefaction, ground surface motion, and calculated lateral spreading were compared based on the simulation results. e mechanism of lateral spreading induced by liquefaction was discussed in the paper.

e Dynamic Motion Beneath the Liquefied Soil and
Triggering Time of Liquefaction. By conducting the effective stress analysis of the site using borehole information, the acceleration time history, the pore water pressure-time history at different depths, and the triggering time of liquefaction can be obtained using effective stress analysis. e one-dimensional effective stress analysis is conducted using the modified Konder-Zelasko model implemented in the computer code D-mod 2000 (GeoMotions [21]). e modified Konder-Zelasko model (the MKZ model) is developed based on the stress-strain relationship by the Konder-Zelasko model (Konder and Zelasko [22]). e Masing's criteria is used to characterize the dynamic behavior of the liquefiable soil during the first cyclic loading in equation (1), where τ is the given stress, c is the given strain, G m0 is the initial shear modulus, and τ m0 is the initial shear strength of the soil: e modified Konder-Zelasko model introduces two more fitting parameters to increase the accuracy of the Konder-Zelasko model. In equation (2), the initial loading curve of the modified Konder-Zelasko model is expressed, where the parameters denoted with * are normalized values, and the parameters β and s are the fitting parameters, respectively: e viscous damping ratio representing the difference between the measured damping value and analytical damping value is used in the model to improve the accuracy of the equivalent damping ratio at the large and moderate strain level. In equation (3), the equivalent viscous damping ratio is expressed, where λ 0 is a function of a given strain c and c c0 is the cyclic shear strain amplitude, respectively: e degraded modified backbone curve is adopted with Masing's criteria during the subsequent cycles. In equation (4), the strain-stress relationship of the soil is expressed, where the parameters denoted with * are the normalized values of the shear modulus and shear stress, respectively. G * mt is the normalized shear modulus at time t, and τ * mt is the cyclic shear stress at time t:

2
Shock and Vibration e residual excess pore pressure is introduced as the degradation parameter. In equations (5) and (6), the normalized shear modulus and shear stress of the soil at time t are expressed, where u * is the normalized residual excess pore water pressure and v is a fitting parameter in equation (6): To solve the equation of dynamic motion in equation (7) In equation (8), the viscous damping is shown, where α R and β R are the coefficients of Rayleigh damping and m and k are the matrix elements of the mass and stiffness, respectively.

2.2.
e Residual Shear Strength and Yield Acceleration. Ozener [20] proposed equation (9) to correlate the normalized shear wave velocity to the residual shear strength for the soil with shear wave velocity less than 250 m/s: where S u is the residual shear strength of the liquefiable soil in kPa, σ v0 ′ is the effective stress in kPa, and V s is the shear wave velocity of the liquefiable soil in m/s. e residual shear strength of the liquefiable soil in this paper is estimated using equation (9), and the corresponding yield acceleration is obtained with the limit equilibrium method based on the soil profile of the site.

Newmark Sliding-Block Prediction Method.
e lateral spreading can be calculated with the rigorous Newmark sliding-block method following the steps below: Step 1: develop the one-dimensional model from the borehole investigation; Step 2: apply the outcrop motion at the base of the one-dimensional model; Step 3: conduct the effective stress analysis and obtain the acceleration time history beneath the liquefied soil; Step 4: obtain the liquefaction time based on the effective stress analysis; Step 5: estimate the residual shear strength and obtain the yield acceleration using limit equilibrium method; Step 6: conduct the Newmark sliding-block method using the motion beneath the liquefied soil (starting from the time of liquefaction) from Step 3 and 4 and the yield acceleration from Step 5.

e Geotechnical Aspect of Wildlife Site
Array. e Wildlife Site Array (Holzer et al. [23]) was one of the field sites that had observed liquefaction, excess pore water pressure generation, and lateral spreading. e Wildlife Site Array was located to the west bank of the Alamo River in Southern California, US. In Figure 1(a), the liquefaction feature of the Wildlife Site Array (Soroush [24]) is shown, and the soil profile of the site (Makdisi [12]) is shown in Figure 1(b). In Figure 2, the instrumentation at the Wildlife Site Array (Holzer et al. [23]) by the United States Geological Survey is shown. e Wildlife Site Array consists of four layers subsequently from the ground surface to the bottom: the silt, the silty sand, the silty clay, and the silt layer. e water table is at 1.25 m below the ground surface. e silty sand below the silt layer was identified as liquefiable soil during the earthquake. Six piezometers were installed to record the pore water pressure of each location, and 2 strong-motion seismometers were installed to record the strong ground motion. e lateral spreading measured by Holzer et al. [23] was 0.18 m, of which the direction was perpendicular to the flowing direction of the Alamo River.

Dynamic Input and Effective Stress Analysis of the Site.
e dynamic input recorded by the strong-motion station SM1 at the depth of −7 m during the earthquake (Earthquake of November 24, 1987) is shown in Figure 3 and used in this paper. e direction of the motion was parallel with the direction of the lateral spreading of the site. As the motion was exactly recorded below the ground surface, it was directly applied at the bottom of the one-dimensional effective model, and the deconvolution analysis can be conducted automatically in the computer code D-mod 2000. Baseline correction and filtering were conducted to avoid possible drift in the displacement time history, which may induce unexpected large lateral spreading.
e one-dimensional model is developed based on the soil profile by Ching [25], as shown in Figure 4. At the bottom of the column of the site, the dynamic input shown in Figure 3 was applied. e soil property for each soil layer used in the nonlinear stress-strain model is summarized in Table 1. Silty sand-1 and Silty sand-2 in Table 1 represent two different parameters of the silty sand varying with depth. In the MKZ model, the unit weight, shear wave velocity, and thickness for each soil layer are the primary input parameters. e initial shear modulus for each soil layer is calculated based on the unit weight and shear velocity. Assign the dynamic properties (i.e., modulus reduction and damping curves) to each soil layer, and the curve-fitting parameters defined in the MKZ model (i.e., parameter β, s, and v) are calculated by a trial-and-error procedure in D-mod 2000.

Lateral Spreading Predicted by the Hybrid Approach.
e shear wave velocity of the Silty sand-2 shown in Table 1 was used to estimate the value of residual shear Shock and Vibration 3 strength. e effective stress over the liquefiable soil was estimated as 61.65 kPa, the shear wave velocity of the liquefiable soil was 120 m/s, and thus, the residual shear strength of the liquefiable soil was 4.62 kPa. By conducting the limit equilibrium analysis using the soil profile shown in Figure 1(b) (developed from Makdisi [12]), the yield acceleration k y of the site was 0.01 g. Based on the triggering time of the liquefaction of 22.9 s for the Silty sand-2 layer, the motion beneath Silty sand-2 considering the liquefaction time was used in the Newmark sliding-block analysis. e normal-direction displacement and inversedirection were 0.33 m and 0.21 m, respectively, so the average calculated lateral spreading was 0.27 m. Comparing with the reported lateral spreading of 0.18 m, the predicted displacement was 1.5 times the reported lateral spreading.

Numerical Simulations to Lateral Spreading
e effective stress analysis uses constitutive models of liquefiable soil under different loading and stress conditions within geotechnical engineering. In the constitutive models, the soil parameters of liquefiable soil were related to the effective shear strength parameters, and the pore water pressure generation is incorporated. Different research studies have been conducted to evaluate the seismic liquefaction or liquefaction-induced lateral spreading. Gu [26] analyzed postearthquake deformation of the lower San Fernando dam with an incremental finite-element method and concluded that the flow failures of the dam were induced by stress redistribution initiated by strain-softening of the liquefied materials. Gu [27] conducted a static finite-element deformation analysis on postearthquake deformation of Wildlife Site Array in the 1987 Superstition Hills earthquake and showed that the lateral spreading was induced by stress redistribution of liquefied soil. Uzuoka et al. [28] proposed a numerical method based on fluid dynamics for the prediction of lateral spreading. Hadush et al. [29] proposed a cubic interpolated pseudoparticle based on a numerical approach and incorporated the method into a fluid dynamics framework to reproduce the dynamic response of liquefaction-induced ground flow from experiments. Elgamal et al. [30] proposed a plasticity-based constitutive model for cyclic mobility analysis and calibrated the model by using laboratory data and centrifuge experiments. e soil response was impacted by the dominant excitation frequency. Kanibir et al. [31] evaluated the lateral spreading on Lake Sapanca shore in the 1999 Kocaeli earthquake with the finite-element method. e finite-element method was accurate if the soil properties were well defined. Valsamis et al. [32] predicted lateral spreading cases with a numerical method and empirical method. It was concluded that the ground deformation was affected by the presence of nonliquefiable soil or interlayered liquefiable soil. Seid-Karbasi and Byrne [9] studied the low-permeability effects induced by sily layers on the pore water pressure within sloping ground using numerical method. Mayoral et al. [33] developed a hybrid numerical procedure for the evaluation of lateral spreading, which considered wave propagation, displacement distribution, and pore water pressure generation. Phillips et al. [34] developed a three-dimensional numerical model and calibrated the model with displacement, acceleration, and pore water pressure records in a centrifuge test of free-field lateral spreading. Kamai and Boulanger [35] simulated a centrifuge test and analyzed dissipation pattern, lateral spreading, and shear strain localization of the model to verify the proposed numerical modeling approach. Montassar and De Buhan [36] developed an identification procedure to evaluate the undrained shear strength and the Bingham viscosity coefficient of the viscoplastic Bingham media simulating liquefied soil of lateral spreading in two centrifuge tests. Howell et al. [37] conducted numerical simulation on centrifuge tests to predict the site response and lateral spreading and verified the improvement by prefabricated vertical drains. Munter et al. [38] predicted the potential lateral spreading at a site in Turkey during the 1999 Kocaeli earthquake by two-dimensional nonlinear deformation analyses with the magnitude of centimeters. Maza [39] conducted numerical modeling of the lateral spreading induced by liquefaction in the 2010 Maule earthquake. Good agreement between the simulation and field investigations was observed. Ghasemi-Fare and Pak [40,41] simulated a centrifuge test of the gently sloping  Shock and Vibration 5 liquefiable ground using a fully coupled numerical analysis and proposed an equation predicting the maximum lateral spreading based on numerical results. Ramirez et al. [42] validated the two numerical constitutive models using a centrifuge test and concluded applicability of the two models. Manzari et al. [43] simulated the centrifuge test of LEAP (the liquefaction experiment and analysis project), and the simulations showed good agreement with the measured response under different amplitudes of base excitations. Ziotopoulou [44] simulated the LEAP centrifuge tests (the liquefaction experiment and analysis project) using the PM4sand model to capture the dynamic response of sloping ground. Tropeano et al. [45] proposed a numerical model and reproduced the pore water generation and dissipation. As a summary, the numerical simulation is a promising method to reproduce the characteristics of the liquefaction and site response under different loading and stress conditions, but the calculation depends on the pore water pressure generation, stress-strain relationship, stiffness degradation, and the other factors. Among these constitutive models for liquefiable sand, PM4sand (Boulanger [46]) and UBC3D-PLM (Petalas and Galavi [11]) are commonly used in the effective stress analysis for liquefaction problems. e plasticity constitutive model of PM4sand is a stress-ration controlled, critical state compatible, and bounding-surface model for sand, the model is a twodimensional constitutive model, and the out-of-plane direction (i.e., third dimension) is considered as elastic evolution. e UBC3D-PLM model is an elastoplastic model developed based on UBCSAND, which is a two-dimensional plasticity model with a hyperbolic strain-hardening rule based on Duncan-Chang criterion. In this paper, PM4sand and UBC3D-PLM constitutive models are used to analyze the liquefaction-induced lateral spreading of the Wildlife Site Array.

Effective Stress Analysis Using PM4sand
Model. One of the advantages of the PM4sand model is that the model can be defined using only three main parameters. e other parameters in the model are recommended to be set as default values by the manual. e three primary parameters include: the apparent relative density, D r , the shear modulus coefficient, G 0 , and the contraction rate parameter, h p0 . D r can be estimated based on the SPT value of the liquefiable soil, which controls the dilatancy and stress-strain relationship, as expressed in the following equation: G 0 controls the shear modulus at small strains and restrains the deviatoric and volumetric increments. It can be expressed in equation (11) in terms of the SPT value of the liquefiable soil: e contraction rate parameter adjusts the plastic volumetric strain during the contraction. It can be obtained when calibrating the cyclic resistance ratio to a target value.
In Table 2, the thickness of soil layers and the soil parameters used in the analysis are presented. e dry density, shear wave velocity, and shear modulus are based on the same parameters used in Table 1. In Table 3, based on the SPT value of the liquefiable soil, the PM4sand parameters for the liquefiable soil are summarized, the calibration of contraction rate parameter is based on the target cyclic resistance ratio for 15 cycles of loading for an earthquake of moment magnitude 7.5 based on the SPT value and the corresponding CRR (cyclic resistance ratio) value. Silty sand-1 and Silty sand-2 in Table 3 represent two different PM4sand parameters of the silty sand.
In Figure 5, the mesh of the Wildlife Site Array generated by FLAC 2D [47] is shown. e dimension of the model is large enough to avoid the influence on the boundary condition. e motion recorded from the SM1 strong-motion station was directly applied at the bottom. A Rayleigh damping ratio of 0.5% was assigned in the model to absorb the high-frequency content of the input motion. Regarding the static and dynamic boundary conditions, the x and y directions were fixed at the bottom of the model, and the xdirection for each lateral boundary was fixed. A free-field boundary condition was applied to avoid the reflection of motions at the boundaries and the reflection internally. A quiet boundary (viscous boundary) was applied at the bottom of the model to absorb the increment of stress induced by dynamic loading. Before the dynamic analysis, the static analysis based on the Mohr-Coulomb criterion was conducted to achieve the initial equilibrium of the model. e groundwater flow was calculated to establish the phreatic groundwater surface and pore water pressure distribution.

Effective Stress Analysis Using UBC3D-PLM Model.
e UBC3D-PLM constitutive model was used in the effective stress analysis to model liquefiable soil. e UBC3D-PLM model is the three-dimensional formulation of UBCSAND. It contains the Mohr-Coulomb yield surface in the three-dimension principal stress space for the primary loading and a yield surface with kinematic hardening rule for the secondary loading. A nonassociated plastic potential function is developed based on Drucker-Prager's criterion. e flow rule is based on the modified stress-dilatancy theory developed by Rowe [48]. e initial calculation was created to generate the initial stress in soil layers using the K0 procedures, and then the dynamic analysis was conducted. In dynamic analysis, as the development of the hardening soil small model overcame the disadvantage that the hardening soil model used high stiffness at small strain levels (i.e., the strain was smaller than 10e − 5), the nonliquefiable soil was modeled using hardening soil small constitutive model [49].
In Table 4, based on the relative density of each soil layer, the hardening soil small model parameters are summarized. In equations (12)  is the reference shear stiffness at small strain, c 0.7 is the shear strain at which G reduced to 72.2%, m is the rate of stress-level dependency in stiffness behavior, K ref nc is the horizontal over vertical stress ratio in primary 1 − D compression, R f is the failure ratio, and D r is the relative density of the soil.
In the UBC3D-PLM constitutive model, it requires 15 primary parameters in PLAXIS 3D (Brinkgreve et al. [49]) as input parameters. e calibration of the input parameters depends on the SPT value of the liquefiable soil. Among the 15 parameters, the stiffness parameters k * e B , k * e G , and k * P G , in equations (20)- (22), are correlating with the normalized SPT value, (N 1 ) 60 of the liquefiable soil: where k * e B is the elastic bulk modulus, k * e G is the elastic shear modulus, and k * P G is the plastic shear modulus. e constant volume friction angle, φ cv , the peak friction angle, φ p , and the cohesion, c, can be derived from the drained triaxial test or direct shear simple test. e peak friction angle is expressed in equation (23). In the analysis, the cohesion was equal to zero: ree modulus indexes are recommended to be default values: the elastic bulk modulus index, me � 0.5; the elastic  shear modulus index, ne � 0.5; and the plastic shear modulus index, np � 0.5.
Based on the constitutive model manual, the tension cutoff σ t is equal to zero. e atmospheric pressure, p ref , is recommended to be default value as 100 kP by the constitutive model manual. e other 4 parameters include the failure ratio R f , the densification factor f dens , the SPT value of liquefiable soil, (N 1 ) 60 and the postliquefaction factor, f Epost .
In Table 5, the UBC3D-PLM model parameters and the other input parameters [50] for liquefiable soil of Wildlife Site Array are shown. Silty sand-1 and Silty sand-2 in Table 5 represent different parameters of the silty sand.
e Rayleigh damping coefficients of soil layers were determined based on the two target frequencies, f 1 and f 2 , respectively. e target frequency f 1 was calculated based on the average shear velocity of the soil layers and the height of soil layers, and the target frequency f 2 was obtained based on the maximum Fourier amplitude of the input motion and the target frequency f 1 . In equations (24) and (25), the target frequencies f 1 and f 2 are expressed, respectively: where V s,average is the average shear wave velocity of the soil profile and H is the thickness of the soil profile. Once f 1 is obtained, the target frequency f 2 can be calculated. In equation (25), f fud is the maximum Fourier amplitude of the input motion, which can be obtained through the Fourier transform of the input motion. Regarding the input motion, the maximum Fourier amplitude is 0.005 at the frequency of 1.77 Hz. e target frequency f 1 � 3.85 Hz and f 2 � 2.175 Hz. e input motion in Figure 3 was applied as a prescribed displacement of 1.0 m multiplied by dynamic multipliers, at the bottom of the autogenerated model mesh in Figure 6. e numerical model was a pseudo-3D one in essence, as the model at y direction was extended from a 2-D model. In the static analysis, the displacement of x and y-direction was fixed along the lateral boundaries and the bottom boundary was fixed at three directions (i.e., x, y, and z-direction). In the dynamic analysis, the free-field boundary was assigned at the x-direction and y-direction boundary, and the compliant base was assigned at the bottom (z-direction boundary).

Excess Pore Water Pressure Generation.
e variation of excess pore water pressure ratios (i.e., the ratio of excess pore pressure to the effective vertical stress) and triggering time of liquefaction of Silty sand-1 at the depth of −3.5 m by the three analyses are shown in Figure 7. Pore water pressure recorded by P2 is also plotted in Figure 7, indicating that the liquefaction occurred at the time of 40 s (note: the record was shifted forward by 15 seconds corresponding to the initiation of strong motion). When applying the MKZ constitutive model, the pore water pressure increased gradually. e ratio was equal to 0.95 at the time of 84 s. e back curve degraded as a function of the increasing of the pore water pressure, which was a strain-based model developed from the result of undrained, strain-controlled cyclic shear tests; thus, the  When applying the PM4sand constitutive model, the pore water pressure ratio increased when the dynamic time reached 18.50 s. When applying the UBC3D-PLM constitutive model, the variation of the excess pore water pressure ratio showed that the site liquefied at the time of 8.19 s. Different from the MKZ constitutive model, the pore water pressure generations in PM4sand and UBC3D-PLM constitutive models were depending on the volumetric strain change of the soil rather than threshold strain used by the MKZ model. Furthermore, when the UBC3D-PLM model was used in the three-dimensional analysis, the pore water pressure was generated using the two yield surfaces in the hardening process and the soil densification effect was considered.

Lateral Spreading Induced by Liquefaction
e displacement of the point at the edge of the free surface was recorded and regarded as liquefaction-induced lateral spreading in the simulations. In Figure 8, the variations of lateral spreading versus dynamic time for PM4sand and UBC3D-PLM constitutive models show that the two constitutive models predict 0.94 m (the UBC3D-PLM constitutive model) and 0.72 m (the PM4sand constitutive model), respectively. e two different predictions are induced by the postliquefaction effect. e displacement became greater due Shock and Vibration to the soil softening, and the triggering time of liquefaction will affect the accumulation of the permanent displacement. erefore, the effective stress analysis using the UBC3D-PLM constitutive model predicted a more conservative value comparing to the PM4sand constitutive model.

Hilbert-Huang Transform.
e Hilbert-Huang transform [51] is a signal-processing technique for the nonlinear and nonstationary signals. e signal is firstly decomposed into different simple intrinsic modes of oscillation using empirical mode decomposition (EMD), and then each mode is analyzed by performing Hilbert transform to obtain the Hilbert spectrum. In each intrinsic mode, an intrinsic function (IMF) is defined by the following two conditions: (1) the number of extreme and the number of zero-crossings are equal or at most differ by one; (2) the mean value of the envelope defined by the local maxima or minima is zero. Any signal X(t) can be expressed in equation (26), where c j is the jth IMF component and r n is the nth residual: e Hilbert transform of each intrinsic function, c j , is expressed in equation (27), where P is the Cauchy principal value. e analytical signal, Z j , can be formed by combing c j and d j , as expressed in equation (28): where a j (t) � , a j (t) is the time-dependent instantaneous Hilbert amplitude, and θ j (t) is the instantaneous phase of c j (t).
e instantaneous frequency is expressed in equation (29). By applying the Hilbert transform to the IMFs in equation (29), the signal X(t) is expressed in equation (30), where RP is the real part of the value to be calculated: e frequency-time distribution of the amplitude is representing the Hilbert amplitude spectrum, as expressed in equation (31), where H j is the jth Hilbert spectrum: e marginal spectrum, h(w), is expressed in equation (32), which is the integral of the Hilbert spectrum with respect to time, where T is the time duration of the signal. h j (w) is the jth marginal spectrum: e Hilbert-Huang transform was conducted for analyzing the nonlinear response of the site in this paper.
e Hilbert-Huang spectra of the recorded ground surface motion and the motion recorded at the depth of −7.5 m were obtained. Figure 9 shows the frequency-time-amplitude distributions of the two acceleration time histories. e maximum amplitude of the Hilbert-Huang spectrum is at 8.42 Hz corresponding with the time of 15.63 s for the silt layer. For the silty clay layer, the maximum amplitude of the Hilbert-Huang spectrum is at 1.0 Hz corresponding with the time of 18.59 s. e differences between the two Hilbert-Huang spectra show that the  amplification due to the soft soil is observed from the depth of 7.5 m below the ground surface to the ground surface. e lowfrequency motion is transmitted to the ground surface and results in a shift of the frequency content due to the liquefaction. e high-frequency content of the motion is identified at the ground surface after the liquefaction occurs. e marginal amplitude spectra of the two motions are obtained and shown in Figure 10 to examine the liquefaction effects. e two motions have the same peak marginal amplitudes. Both of them have the same corresponding frequency of 0.3 Hz. For the motion recorded at the depth of 7.5 m below the ground surface, two peaks are identified, and the second peak amplitude is at the frequency of 1.25 Hz. is downshift phenomenon of the frequency for motion recorded below the ground surface may be attributed to the soil dilatancy induced by the liquefaction effect. is is consistent with the observations of Hilbert-Huang transforms to the two motions.

e Liquefaction Effects on Lateral
Spreading. e ground surface motion from the analysis using the PM4sand constitutive model was analyzed to obtain response spectra. Two different types of motions were used to investigate the liquefaction effect. e first motion is part of the ground surface motion, but it is only corresponding to the time from the beginning of the motion to the triggering time of liquefaction. e second motion is the entire ground surface motion. Figure 11 shows the two response spectra. e red response spectrum is for the entire motion of the ground surface motion and the black one is the spectrum corresponding with the motion from the beginning to the triggering time of liquefaction. It can be seen from Figure 11 that, from the short-period 0.01 s to 0.6 s, the two spectra coincide with each other. When the period is greater than 0.6 s, the spectral amplitude of the entire motion is greater than the other one. e increase of the spectral amplitude indicates that the liquefaction amplifies the low-frequency content of the motion. is causes that the transmitting of the high-frequency content to the ground surface affects the site response and then affects the displacement of lateral spreading. erefore, the effects on the site response by the triggering time of liquefaction are consistent with the accumulation of lateral spreading and the observation from the Hilbert-Huang transform to the motions.

Conclusions
It is important to predict the magnitude of lateral spreading, which is one of the seismic phenomena induced by liquefaction. A hybrid approach based on the Newmark sliding-block method was proposed in the paper to predict the liquefactioninduced lateral spreading. e lateral spreading of Wildlife Site Array during the 1987 Superstition earthquake was reanalyzed to evaluate the proposed method. e same case history was used to investigate the capability of the constitutive models simulating the soil behaviors and liquefaction-induced lateral spreading under the cyclic loading. Based on the analyses results, the following conclusions were made. e hybrid approach used the Newmark sliding-block model and considered the residual shear strength of the liquefiable soil and the displacement characteristics of lateral spreading. By analyzing the borehole profile using one-dimensional effective stress analysis, the lateral spreading was calculated. It predicted 1.5 times the reported lateral spreading for the Wildlife Site Array. e pore water pressure ratio and lateral spreading versus dynamic time were obtained by conducting  numerical simulations using PM4sand and UBC-PLM 3D constitutive models. Additionally, the pore water pressure ratio time history from one-dimensional effective stress analysis is compared with the ratios from numerical simulations. It indicated that the lateral spreading increased with the time and depended on the liquefaction time, which varied with the constitutive models used in the analysis. e analyses using PM4sand and UBC-PLM 3D constitutive models predicted 5.2 and 4 times the reported lateral spreading for the Wildlife Site Array, respectively. e Hilbert-Huang transform to site response was conducted using the ground surface motion and the motion recorded below the ground surface. e results showed that the liquefaction amplified the low-frequency content of the motion and caused a shift of the frequency content. is phenomenon was too observed by comparing the spectral amplitudes of the entire ground surface motion with the motion from the beginning to the triggering time of liquefaction in the numerical simulation using PM4sand constitutive model.
When conducting analysis using the hybrid method, the direction of the input motion was parallel with the direction of lateral spreading, so the input motion is exactly the motion that induces lateral spreading during the earthquake. It is the only one used to calculate the liquefaction-induced lateral spreading in this paper. For the condition that there is no available downhole motion, the deconvolution analysis using the free-surface motion is required to obtain the outcrop motion. e analyses conducted in the paper were based on one case, more analyses were recommended to be carried to evaluate the proposed method and understand the performances of different constitutive models of liquefiable soil at the large strain level.

Data Availability
e data used to support the findings of the study are included in the paper.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.