# A Hybrid Approach Calculating Lateral Spreading Induced by Seismic Liquefaction

Authors: Yanxin Yang; Jianguo Liu; Jianlin Ma; Qinke Wang
Journal: Shock and Vibration
 (2020)

---

## Abstract

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. The limit equilibrium analysis was conducted to determine the yield acceleration corresponding with the residual shear strength of liquefied soil. The 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. The 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. The 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. The 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.


---

## Body

## 1. 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. The 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–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. The prediction methods of lateral spreading mainly include empirical methods [6–8], numerical methods [9–11], Newmark sliding-block methods [12], and probabilistic methods [13–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. The 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. The 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. The 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.The 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. The residual shear strength could be estimated from the SPT, CPT value, or shear wave velocity of the liquefiable soil. Several methods [17–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. Though 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.” There 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. The 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. The triggering time of liquefaction, ground surface motion, and calculated lateral spreading were compared based on the simulation results. The mechanism of lateral spreading induced by liquefaction was discussed in the paper.


## 2. The Hybrid Prediction Approach

### 2.1. The 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. The one-dimensional effective stress analysis is conducted using the modified Konder–Zelasko model implemented in the computer code D-mod 2000 (GeoMotions [21]). The modified Konder–Zelasko model (the MKZ model) is developed based on the stress-strain relationship by the Konder–Zelasko model (Konder and Zelasko [22]). The 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, γ is the given strain, Gm0 is the initial shear modulus, and τm0 is the initial shear strength of the soil:(1)τ=fγ=Gm01+Gm0/τm0γ.The 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:(2)τ=fγ=Gm0∗γ1+βGm0∗/τm0∗γs.The 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 γ and γc0 is the cyclic shear strain amplitude, respectively:(3)λ0=4π∫0γc0fγdγfγc0γc0−2π.The 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. Gmt∗ is the normalized shear modulus at time t, and τmt∗ is the cyclic shear stress at time t:(4)τ∗=f∗γ=Gmt∗γ1+βGmt∗/τmt∗γs.The 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):(5)Gmt∗=Gm0∗1−u∗,(6)τmt∗=τm0∗1−u∗ν.To solve the equation of dynamic motion in equation (7), the time-domain Newmark β algorithm is used, where the [M], [C], and [K] are the mass matrix, viscous damping matrix, and nonlinear stiffness matrix, respectively, and u, u˙ and u¨ are the relative displacements, relative velocities, and relative accelerations from mass to the base, respectively:(7)Mu¨+Cu˙+Ku=ft,(8)c=αRm×βRk.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. The Residual Shear Strength and Yield Acceleration

Özener [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:(9)Suσv0′=0.0218exp0.0103Vs,where Su is the residual shear strength of the liquefiable soil in kPa, σv0′ is the effective stress in kPa, and Vs is the shear wave velocity of the liquefiable soil in m/s.The 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.


### 2.3. Newmark Sliding-Block Prediction Method

The 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.


## 2.1. The 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. The one-dimensional effective stress analysis is conducted using the modified Konder–Zelasko model implemented in the computer code D-mod 2000 (GeoMotions [21]). The modified Konder–Zelasko model (the MKZ model) is developed based on the stress-strain relationship by the Konder–Zelasko model (Konder and Zelasko [22]). The 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, γ is the given strain, Gm0 is the initial shear modulus, and τm0 is the initial shear strength of the soil:(1)τ=fγ=Gm01+Gm0/τm0γ.The 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:(2)τ=fγ=Gm0∗γ1+βGm0∗/τm0∗γs.The 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 γ and γc0 is the cyclic shear strain amplitude, respectively:(3)λ0=4π∫0γc0fγdγfγc0γc0−2π.The 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. Gmt∗ is the normalized shear modulus at time t, and τmt∗ is the cyclic shear stress at time t:(4)τ∗=f∗γ=Gmt∗γ1+βGmt∗/τmt∗γs.The 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):(5)Gmt∗=Gm0∗1−u∗,(6)τmt∗=τm0∗1−u∗ν.To solve the equation of dynamic motion in equation (7), the time-domain Newmark β algorithm is used, where the [M], [C], and [K] are the mass matrix, viscous damping matrix, and nonlinear stiffness matrix, respectively, and u, u˙ and u¨ are the relative displacements, relative velocities, and relative accelerations from mass to the base, respectively:(7)Mu¨+Cu˙+Ku=ft,(8)c=αRm×βRk.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. The Residual Shear Strength and Yield Acceleration

Özener [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:(9)Suσv0′=0.0218exp0.0103Vs,where Su is the residual shear strength of the liquefiable soil in kPa, σv0′ is the effective stress in kPa, and Vs is the shear wave velocity of the liquefiable soil in m/s.The 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.


## 2.3. Newmark Sliding-Block Prediction Method

The 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.


## 3. Lateral Spreading Predicted by the Hybrid Approach

### 3.1. The Geotechnical Aspect of Wildlife Site Array

The 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. The 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. The 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. The water table is at 1.25 m below the ground surface. The 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. The 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.Figure 1
The liquefaction feature and soil profile of Wildlife Site Array. (a) Liquefaction feature of the Wildlife Site Array (after Soroush [24]). (b) The soil profile of Wildlife Site Array (after Makdisi [12]).
(a)(b)Figure 2
The cross section and instrumentation at the Wildlife Site Array (after Holzer [23]).


### 3.2. Dynamic Input and Effective Stress Analysis of the Site

The 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 Figure3 and used in this paper. The 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.Figure 3
The dynamic input (at the depth of −7 m).The 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. The 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. The 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.Figure 4
The one-dimensional model of Wildlife Site Array (after Ching [25]).Table 1
The thickness of soil layers and the soil parameters.
Soil layersThickness (m)Saturated unit weight (kN/m3)Shear wave velocity (m/s)Initial shear modulus (kPa)Shear stress at reference strain (kPa)Curve-fitting parameter-βCurve-fitting parameter-sCurve-fitting parameter-v controls cyclic degradation of saturated sandSilt2.519.69016171.07.921.990.811.0Silty sand-11.019.611024173.830.944.340.811.0Silty sand-22.719.612024183.639.425.280.811.0Silty clay1.319.618064756.1103.611.170.810.0


### 3.3. Lateral Spreading Predicted by the Hybrid Approach

The shear wave velocity of the Silty sand-2 shown in Table1 was used to estimate the value of residual shear strength. The 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 ky 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. The normal-direction displacement and inverse-direction 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.


## 3.1. The Geotechnical Aspect of Wildlife Site Array

The 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. The 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. The 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. The water table is at 1.25 m below the ground surface. The 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. The 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.Figure 1
The liquefaction feature and soil profile of Wildlife Site Array. (a) Liquefaction feature of the Wildlife Site Array (after Soroush [24]). (b) The soil profile of Wildlife Site Array (after Makdisi [12]).
(a)(b)Figure 2
The cross section and instrumentation at the Wildlife Site Array (after Holzer [23]).


## 3.2. Dynamic Input and Effective Stress Analysis of the Site

The 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 Figure3 and used in this paper. The 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.Figure 3
The dynamic input (at the depth of −7 m).The 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. The 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. The 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.Figure 4
The one-dimensional model of Wildlife Site Array (after Ching [25]).Table 1
The thickness of soil layers and the soil parameters.
Soil layersThickness (m)Saturated unit weight (kN/m3)Shear wave velocity (m/s)Initial shear modulus (kPa)Shear stress at reference strain (kPa)Curve-fitting parameter-βCurve-fitting parameter-sCurve-fitting parameter-v controls cyclic degradation of saturated sandSilt2.519.69016171.07.921.990.811.0Silty sand-11.019.611024173.830.944.340.811.0Silty sand-22.719.612024183.639.425.280.811.0Silty clay1.319.618064756.1103.611.170.810.0


## 3.3. Lateral Spreading Predicted by the Hybrid Approach

The shear wave velocity of the Silty sand-2 shown in Table1 was used to estimate the value of residual shear strength. The 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 ky 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. The normal-direction displacement and inverse-direction 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.


## 4. Numerical Simulations to Lateral Spreading

The 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. The 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. The 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 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. The plasticity constitutive model of PM4sand is a stress-ration controlled, critical state compatible, and bounding-surface model for sand, the model is a two-dimensional constitutive model, and the out-of-plane direction (i.e., third dimension) is considered as elastic evolution. The 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.


### 4.1. 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. The other parameters in the model are recommended to be set as default values by the manual. The three primary parameters include: the apparent relative density,Dr, the shear modulus coefficient, G0, and the contraction rate parameter, hp0. Dr 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:(10)Dr=N16046.G0 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:(11)G0=167N160+2.5.The 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 Table2, the thickness of soil layers and the soil parameters used in the analysis are presented. The 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.Table 2
The thickness of soil layers and the soil parameters.
Soil layersThickness (m)Dry density (kg/m3)Shear wave velocity (m/s)Cohesion (kPa)Internal friction angle (degrees)Bulk modulus (kPa)Shear modulus (kPa)Silt2.51440909.533.64.04e41.86e4Silty sand-11.014681102.033.05.25e42.42e4Silty sand-22.714781202.033.05.25e42.42e4Silty clay4.8152318021.530.01.20e55.52e4Table 3
The parameters for the PM4sand model [10].
Soil layerSPTN160Dr (%)G0hp0Silty sand-16.64383501.90Silty sand-215.47582660.86In Figure5, the mesh of the Wildlife Site Array generated by FLAC 2D [47] is shown. The dimension of the model is large enough to avoid the influence on the boundary condition. The 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 x-direction 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. The groundwater flow was calculated to establish the phreatic groundwater surface and pore water pressure distribution.Figure 5
The mesh of Wildlife Site Array using the PM4sand model.


### 4.2. Effective Stress Analysis Using UBC3D-PLM Model

The UBC3D-PLM constitutive model was used in the effective stress analysis to model liquefiable soil. The 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. The flow rule is based on the modified stress-dilatancy theory developed by Rowe [48].The 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 Table4, based on the relative density of each soil layer, the hardening soil small model parameters are summarized. In equations (12)–(19), the 8 parameters of the hardening soil small model are presented: secant stiffness from the triaxial test at a reference pressure, E50ref; tangent stiffness from the oedometer test at a reference pressure, Eoderef; unloading/reloading stiffness, Eurref; reference shear stiffness at small strain, G0ref; shear strain at which G reduced to 72.2%, γ0.7; the rate of stress-level dependency in stiffness behavior, m; the horizontal over vertical stress ratio in primary 1 − D compression, Kncref, and failure ratio, Rf; and the relative density of the soil, Dr:(12)E50ref=60∗Dr,(13)E50ref=Eoedref=60∗Dr,(14)Eurref=2Eoedref,(15)G0ref=60+68∗Dr100,(16)γ0.7=2−Dr100∗10−4,(17)m=0.7−Dr320,(18)Kncref=1−sinφ′,(19)Rf=1−Dr800,where E50ref is secant stiffness from the triaxial test at a reference pressure, Eoderef is tangent stiffness from the oedometer test at a reference pressure, Eurref is unloading/reloading stiffness, G0ref is the reference shear stiffness at small strain, γ0.7 is the shear strain at which G reduced to 72.2%, m is the rate of stress-level dependency in stiffness behavior, Kncref is the horizontal over vertical stress ratio in primary 1 − D compression, Rf is the failure ratio, and Dr is the relative density of the soil.Table 4
Soil parameters of the hardening soil small model.
Soil layersSymbolSilt-1Silt-2Silty clayThickness (m)1.21.34.8Relative density (%)Dr36.7437.547.43Secant stiffness (kPa)E50ref220442250028458Tangent stiffness (kPa)Eoderef220442250028458Unloading/reloading stiffness (kPa)Eurref440884500056916Reference shear stiffness (kPa)G0ref849838550092252Shear strain at whichG has reduced to 72.2%γ0.70.0001630.0001630.000153Rate of stress-level dependencym0.5851880.5828130.551781Horizontal over vertical stress ratioKncref0.4613390.4599430.441838Failure ratioRf0.9540.9530.941In the UBC3D-PLM constitutive model, it requires 15 primary parameters in PLAXIS 3D (Brinkgreve et al. [49]) as input parameters. The calibration of the input parameters depends on the SPT value of the liquefiable soil. Among the 15 parameters, the stiffness parameters kB∗e, kG∗e, and kG∗P, in equations (20)–(22), are correlating with the normalized SPT value, N160 of the liquefiable soil:(20)kB∗e=0.7∗kG∗e,(21)kG∗e=21.7∗20∗N1600.3333,(22)kG∗P=kG∗e∗N1602∗0.003+100,where kB∗e is the elastic bulk modulus, kG∗e is the elastic shear modulus, and kG∗P is the plastic shear modulus.The 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. The peak friction angle is expressed in equation (23). In the analysis, the cohesion was equal to zero:(23)φp=φcv+N16010+max0;N160−155.Three 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. The atmospheric pressure, pref, is recommended to be default value as 100 kP by the constitutive model manual. The other 4 parameters include the failure ratio Rf, the densification factor fdens, the SPT value of liquefiable soil, N160 and the postliquefaction factor, fEpost.In Table5, 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.Table 5
Soil parameters of the UBC3D-PLM model [50].
Soil layersSymbolSilty sand-1Silty sand-2Thickness (m)1.02.7Young’s modulus (kN/m2)Eref7.28e47.28e4Poisson’s ratiov0.30.3Unsaturated unit weight (kN/m3)γusta19.719.7Saturated unit weight (kN/m3)γsat21.821.8Void ratioe0.740.74Constant volume friction angle (degrees)φcv22.022.0Peak friction angle (degrees)φp22.823.1Cohesion (kPa)C00Elastic shear moduluskG∗e864.6954.1Plastic shear moduluskG∗P250.0424.7Elastic bulk moduluskB∗e598.2667.9Elastic shear modulus indexne0.50.5Plastic shear modulus indexnp0.50.5Elastic bulk modulus indexme0.50.5Failure ratioRf0.8110.771Atmospheric pressure (kPa)pref100100Tension cutoff (kPa)σt00Densification factorfdens0.450.45SPT valueN1607.6510.65Post liquefaction factorfEpost0.20.2Permeability (m/s)K5.0e − 72.0e − 8Tangent stiffness for oedometer (kPa)Eoedref9800098000Dilatancy angleΨ19.018.0The Rayleigh damping coefficients of soil layers were determined based on the two target frequencies,f1 and f2, respectively. The target frequency f1 was calculated based on the average shear velocity of the soil layers and the height of soil layers, and the target frequency f2 was obtained based on the maximum Fourier amplitude of the input motion and the target frequency f1. In equations (24) and (25), the target frequencies f1 and f2 are expressed, respectively:(24)f1=Vs,average4H,(25)f2=ffudf1,where Vs,average is the average shear wave velocity of the soil profile and H is the thickness of the soil profile. Once f1 is obtained, the target frequency f2 can be calculated. In equation (25), ffud 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. The target frequency f1=3.85Hz and f2=2.175Hz.The input motion in Figure3 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. The 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).Figure 6
The three-dimensional mesh of Wildlife Site Array.


## 4.1. 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. The other parameters in the model are recommended to be set as default values by the manual. The three primary parameters include: the apparent relative density,Dr, the shear modulus coefficient, G0, and the contraction rate parameter, hp0. Dr 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:(10)Dr=N16046.G0 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:(11)G0=167N160+2.5.The 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 Table2, the thickness of soil layers and the soil parameters used in the analysis are presented. The 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.Table 2
The thickness of soil layers and the soil parameters.
Soil layersThickness (m)Dry density (kg/m3)Shear wave velocity (m/s)Cohesion (kPa)Internal friction angle (degrees)Bulk modulus (kPa)Shear modulus (kPa)Silt2.51440909.533.64.04e41.86e4Silty sand-11.014681102.033.05.25e42.42e4Silty sand-22.714781202.033.05.25e42.42e4Silty clay4.8152318021.530.01.20e55.52e4Table 3
The parameters for the PM4sand model [10].
Soil layerSPTN160Dr (%)G0hp0Silty sand-16.64383501.90Silty sand-215.47582660.86In Figure5, the mesh of the Wildlife Site Array generated by FLAC 2D [47] is shown. The dimension of the model is large enough to avoid the influence on the boundary condition. The 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 x-direction 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. The groundwater flow was calculated to establish the phreatic groundwater surface and pore water pressure distribution.Figure 5
The mesh of Wildlife Site Array using the PM4sand model.


## 4.2. Effective Stress Analysis Using UBC3D-PLM Model

The UBC3D-PLM constitutive model was used in the effective stress analysis to model liquefiable soil. The 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. The flow rule is based on the modified stress-dilatancy theory developed by Rowe [48].The 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 Table4, based on the relative density of each soil layer, the hardening soil small model parameters are summarized. In equations (12)–(19), the 8 parameters of the hardening soil small model are presented: secant stiffness from the triaxial test at a reference pressure, E50ref; tangent stiffness from the oedometer test at a reference pressure, Eoderef; unloading/reloading stiffness, Eurref; reference shear stiffness at small strain, G0ref; shear strain at which G reduced to 72.2%, γ0.7; the rate of stress-level dependency in stiffness behavior, m; the horizontal over vertical stress ratio in primary 1 − D compression, Kncref, and failure ratio, Rf; and the relative density of the soil, Dr:(12)E50ref=60∗Dr,(13)E50ref=Eoedref=60∗Dr,(14)Eurref=2Eoedref,(15)G0ref=60+68∗Dr100,(16)γ0.7=2−Dr100∗10−4,(17)m=0.7−Dr320,(18)Kncref=1−sinφ′,(19)Rf=1−Dr800,where E50ref is secant stiffness from the triaxial test at a reference pressure, Eoderef is tangent stiffness from the oedometer test at a reference pressure, Eurref is unloading/reloading stiffness, G0ref is the reference shear stiffness at small strain, γ0.7 is the shear strain at which G reduced to 72.2%, m is the rate of stress-level dependency in stiffness behavior, Kncref is the horizontal over vertical stress ratio in primary 1 − D compression, Rf is the failure ratio, and Dr is the relative density of the soil.Table 4
Soil parameters of the hardening soil small model.
Soil layersSymbolSilt-1Silt-2Silty clayThickness (m)1.21.34.8Relative density (%)Dr36.7437.547.43Secant stiffness (kPa)E50ref220442250028458Tangent stiffness (kPa)Eoderef220442250028458Unloading/reloading stiffness (kPa)Eurref440884500056916Reference shear stiffness (kPa)G0ref849838550092252Shear strain at whichG has reduced to 72.2%γ0.70.0001630.0001630.000153Rate of stress-level dependencym0.5851880.5828130.551781Horizontal over vertical stress ratioKncref0.4613390.4599430.441838Failure ratioRf0.9540.9530.941In the UBC3D-PLM constitutive model, it requires 15 primary parameters in PLAXIS 3D (Brinkgreve et al. [49]) as input parameters. The calibration of the input parameters depends on the SPT value of the liquefiable soil. Among the 15 parameters, the stiffness parameters kB∗e, kG∗e, and kG∗P, in equations (20)–(22), are correlating with the normalized SPT value, N160 of the liquefiable soil:(20)kB∗e=0.7∗kG∗e,(21)kG∗e=21.7∗20∗N1600.3333,(22)kG∗P=kG∗e∗N1602∗0.003+100,where kB∗e is the elastic bulk modulus, kG∗e is the elastic shear modulus, and kG∗P is the plastic shear modulus.The 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. The peak friction angle is expressed in equation (23). In the analysis, the cohesion was equal to zero:(23)φp=φcv+N16010+max0;N160−155.Three 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. The atmospheric pressure, pref, is recommended to be default value as 100 kP by the constitutive model manual. The other 4 parameters include the failure ratio Rf, the densification factor fdens, the SPT value of liquefiable soil, N160 and the postliquefaction factor, fEpost.In Table5, 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.Table 5
Soil parameters of the UBC3D-PLM model [50].
Soil layersSymbolSilty sand-1Silty sand-2Thickness (m)1.02.7Young’s modulus (kN/m2)Eref7.28e47.28e4Poisson’s ratiov0.30.3Unsaturated unit weight (kN/m3)γusta19.719.7Saturated unit weight (kN/m3)γsat21.821.8Void ratioe0.740.74Constant volume friction angle (degrees)φcv22.022.0Peak friction angle (degrees)φp22.823.1Cohesion (kPa)C00Elastic shear moduluskG∗e864.6954.1Plastic shear moduluskG∗P250.0424.7Elastic bulk moduluskB∗e598.2667.9Elastic shear modulus indexne0.50.5Plastic shear modulus indexnp0.50.5Elastic bulk modulus indexme0.50.5Failure ratioRf0.8110.771Atmospheric pressure (kPa)pref100100Tension cutoff (kPa)σt00Densification factorfdens0.450.45SPT valueN1607.6510.65Post liquefaction factorfEpost0.20.2Permeability (m/s)K5.0e − 72.0e − 8Tangent stiffness for oedometer (kPa)Eoedref9800098000Dilatancy angleΨ19.018.0The Rayleigh damping coefficients of soil layers were determined based on the two target frequencies,f1 and f2, respectively. The target frequency f1 was calculated based on the average shear velocity of the soil layers and the height of soil layers, and the target frequency f2 was obtained based on the maximum Fourier amplitude of the input motion and the target frequency f1. In equations (24) and (25), the target frequencies f1 and f2 are expressed, respectively:(24)f1=Vs,average4H,(25)f2=ffudf1,where Vs,average is the average shear wave velocity of the soil profile and H is the thickness of the soil profile. Once f1 is obtained, the target frequency f2 can be calculated. In equation (25), ffud 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. The target frequency f1=3.85Hz and f2=2.175Hz.The input motion in Figure3 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. The 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).Figure 6
The three-dimensional mesh of Wildlife Site Array.


## 5. Interpretation of Analysis Results

### 5.1. Excess Pore Water Pressure Generation

The 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 Figure7. 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. The ratio was equal to 0.95 at the time of 84 s. The 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 threshold strain of 0.01% for developing the liquefaction was reached at the corresponding time of 84 s for Silty sand-1.Figure 7
Pore water pressure ratio versus time.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.


## 5.1. Excess Pore Water Pressure Generation

The 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 Figure7. 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. The ratio was equal to 0.95 at the time of 84 s. The 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 threshold strain of 0.01% for developing the liquefaction was reached at the corresponding time of 84 s for Silty sand-1.Figure 7
Pore water pressure ratio versus time.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.


## 6. Lateral Spreading Induced by Liquefaction

The displacement of the point at the edge of the free surface was recorded and regarded as liquefaction-induced lateral spreading in the simulations. In Figure8, 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. The two different predictions are induced by the postliquefaction effect. The displacement became greater due to the soil softening, and the triggering time of liquefaction will affect the accumulation of the permanent displacement. Therefore, the effective stress analysis using the UBC3D-PLM constitutive model predicted a more conservative value comparing to the PM4sand constitutive model.Figure 8
Lateral spreading versus dynamic time.


## 7. Hilbert–Huang Transform to Site Response Analysis

### 7.1. Hilbert–Huang Transform

The Hilbert–Huang transform [51] is a signal-processing technique for the nonlinear and nonstationary signals. The 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 Xt can be expressed in equation (26), where cj is the jth IMF component and rn is the nth residual:(26)Xt=∑j=1ncjt+rnt.The Hilbert transform of each intrinsic function,cj, is expressed in equation (27), where P is the Cauchy principal value. The analytical signal, Zj, can be formed by combing cj and dj, as expressed in equation (28):(27)djt=1πP∫cjt′t−t′dt′,(28)Zjt=cjt+djti=ajteiθjt,where ajt=cj2t+dj2t, θjt=arctandjt/cjt, ajt is the time-dependent instantaneous Hilbert amplitude, and θjt is the instantaneous phase of cjt.The instantaneous frequency is expressed in equation (29). By applying the Hilbert transform to the IMFs in equation (29), the signal Xt is expressed in equation (30), where RP is the real part of the value to be calculated:(29)ωjt=dθjtdt,(30)Xt=RP∑j=1najtei∫wjtdt.The 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:(31)Hω,t=∑j=1nH˜jω,t≡∑j=1najt.The marginal spectrum,hw, 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˜jw is the jth marginal spectrum:(32)hω=∑j=1nh˜ω≡∑j=1n∫0Tajtdt.The Hilbert–Huang transform was conducted for analyzing the nonlinear response of the site in this paper. The Hilbert–Huang spectra of the recorded ground surface motion and the motion recorded at the depth of −7.5 m were obtained. Figure9 shows the frequency-time-amplitude distributions of the two acceleration time histories. The 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. The 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. The low-frequency motion is transmitted to the ground surface and results in a shift of the frequency content due to the liquefaction. The high-frequency content of the motion is identified at the ground surface after the liquefaction occurs.Figure 9
Hilbert–Huang spectrum for two earthquake time histories. (a) Hilbert–Huang spectrum for silt layer (Vs = 90 m/s). (b) Hilbert–Huang spectrum for silty clay layer (Vs = 180 m/s).
(a)(b)The marginal amplitude spectra of the two motions are obtained and shown in Figure10 to examine the liquefaction effects. The 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. This downshift phenomenon of the frequency for motion recorded below the ground surface may be attributed to the soil dilatancy induced by the liquefaction effect. This is consistent with the observations of Hilbert–Huang transforms to the two motions.Figure 10
Marginal spectra for two earthquake time histories.


### 7.2. The Liquefaction Effects on Lateral Spreading

The 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. The 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. The second motion is the entire ground surface motion.Figure11 shows the two response spectra. The 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. The increase of the spectral amplitude indicates that the liquefaction amplifies the low-frequency content of the motion. This 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. Therefore, 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.Figure 11
Ground surface response spectra by two-dimensional analysis.


## 7.1. Hilbert–Huang Transform

The Hilbert–Huang transform [51] is a signal-processing technique for the nonlinear and nonstationary signals. The 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 Xt can be expressed in equation (26), where cj is the jth IMF component and rn is the nth residual:(26)Xt=∑j=1ncjt+rnt.The Hilbert transform of each intrinsic function,cj, is expressed in equation (27), where P is the Cauchy principal value. The analytical signal, Zj, can be formed by combing cj and dj, as expressed in equation (28):(27)djt=1πP∫cjt′t−t′dt′,(28)Zjt=cjt+djti=ajteiθjt,where ajt=cj2t+dj2t, θjt=arctandjt/cjt, ajt is the time-dependent instantaneous Hilbert amplitude, and θjt is the instantaneous phase of cjt.The instantaneous frequency is expressed in equation (29). By applying the Hilbert transform to the IMFs in equation (29), the signal Xt is expressed in equation (30), where RP is the real part of the value to be calculated:(29)ωjt=dθjtdt,(30)Xt=RP∑j=1najtei∫wjtdt.The 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:(31)Hω,t=∑j=1nH˜jω,t≡∑j=1najt.The marginal spectrum,hw, 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˜jw is the jth marginal spectrum:(32)hω=∑j=1nh˜ω≡∑j=1n∫0Tajtdt.The Hilbert–Huang transform was conducted for analyzing the nonlinear response of the site in this paper. The Hilbert–Huang spectra of the recorded ground surface motion and the motion recorded at the depth of −7.5 m were obtained. Figure9 shows the frequency-time-amplitude distributions of the two acceleration time histories. The 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. The 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. The low-frequency motion is transmitted to the ground surface and results in a shift of the frequency content due to the liquefaction. The high-frequency content of the motion is identified at the ground surface after the liquefaction occurs.Figure 9
Hilbert–Huang spectrum for two earthquake time histories. (a) Hilbert–Huang spectrum for silt layer (Vs = 90 m/s). (b) Hilbert–Huang spectrum for silty clay layer (Vs = 180 m/s).
(a)(b)The marginal amplitude spectra of the two motions are obtained and shown in Figure10 to examine the liquefaction effects. The 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. This downshift phenomenon of the frequency for motion recorded below the ground surface may be attributed to the soil dilatancy induced by the liquefaction effect. This is consistent with the observations of Hilbert–Huang transforms to the two motions.Figure 10
Marginal spectra for two earthquake time histories.


## 7.2. The Liquefaction Effects on Lateral Spreading

The 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. The 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. The second motion is the entire ground surface motion.Figure11 shows the two response spectra. The 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. The increase of the spectral amplitude indicates that the liquefaction amplifies the low-frequency content of the motion. This 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. Therefore, 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.Figure 11
Ground surface response spectra by two-dimensional analysis.


## 8. 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 liquefaction-induced lateral spreading. The lateral spreading of Wildlife Site Array during the 1987 Superstition earthquake was reanalyzed to evaluate the proposed method. The 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.The 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.The 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. The 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.The Hilbert–Huang transform to site response was conducted using the ground surface motion and the motion recorded below the ground surface. The results showed that the liquefaction amplified the low-frequency content of the motion and caused a shift of the frequency content. This 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. The 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.



---

*Source: 3769206-2020-09-10.xml*