Factors That Affect Liquefaction-Induced Lateral Spreading in Large Subduction Earthquakes

Liquefaction-induced lateral spreading can induce significant deformations and damage in existing structures, such as ports, bridges, and pipes. Past earthquakes have caused this phenomenon in coastal areas and rivers in many parts of the world. Current lateral spreading prediction models tend to either overestimate or underestimate the actual displacements by a factor of two or more when applied to large subduction earthquake events. The purpose of this study was to identify ground motion intensity measures and soil parameters that better correlate with observed lateral spreading under large-magnitude (Mw ≥ 7.5) subduction earthquakes that have occurred in countries like Chile, Japan, and Peru. A numerical approach was first validated against centrifuge and historical cases and then used to generate parametric models on which statistical analysis was applied. Our results show that cumulative absolute velocity (CAV), Housner intensity (HI), and sustained maximum velocity (SMV) have a reasonably good correlation with lateral spreading for the analyzed cases.


Introduction
Several models, either analytical, empirical, or computational, have been formulated to predict the behavior of liquefiable soils and to anticipate the amount of lateral displacements that can be generated during earthquakes. The purpose has been to provide recommendations for the design of, e.g., foundations and embankments to mitigate losses in the case of future earthquakes (Valsamis et al. (2010) [1]). An evaluation of currently used empirical models (Bartlett et al. (1995) [2], Faris et al. (2006) [3], Zhang et al. (2012) [4], Youd et al. (2002) [5], Rauch et al. (2006) [6]) was made by Tryon (2014) [7] and Williams (2015) [8] using historical cases of large-magnitude subduction earthquakes. They found weaknesses in those empirical models. Firstly, there was a lack of historical case data with earthquakes of moment magnitudes greater than 8 (the only historical case was the 1964 Alaska Mw 9.2 earthquake). Secondly, the term referring to the distance from the site to the seismic source (R) is more challenging to determine when dealing with large magnitude earthquakes.
The main aim of this paper is to investigate, using an appropriate numerical model (Elgamal et al. (2002) [9]), the effects of different geotechnical and seismic parameters on the amount of lateral spreading in free-field conditions during large-magnitude subduction earthquakes.

Background of Lateral Spreading in Subduction Earthquakes
Current techniques used to predict liquefaction-induced lateral spreading are mostly empirical (Bray et al. (2010) [10]). Observations from recent earthquakes have shown that these models become inaccurate when extrapolated beyond their limits, such as large-magnitude events or different fault types. In this section, the phenomena of liquefaction and lateral spreading and the issues with current prediction models, when applied to subduction earthquakes, are described.   (Youd (2018) [12]). Liquefaction is a relevant soil phenomenon for geotechnical design as it may cause local or global failures of foundations and even the collapse of complete structures (Jia (2017) [13]). Among the potential consequences of soil liquefaction, one of the most dangerous ones is lateral spreading. Youd (2018) [12] defined this phenomenon as the horizontal displacement of a soil layer riding on liquefied soil either down a gentle slope or toward a free face like a river channel ( Figure 2). When the underlying soil layer liquefies, the non-liquefied upper soil crust continues moving down until it reaches a new equilibrium position.

Liquefaction and Lateral Spreading
The word liquefaction was first used after the 1964 Niigata Mw7.6 earthquake (Kawata et al. 2018 [11]). This phenomenon is defined as a change of the soil phase from a solid to a liquid state due to pore water pressure increment, and the corresponding loss of effective stress, during an earthquake ( Figure 1). As Youd (2018) [12] indicated, when an earthquake occurs, waves propagate through the soil, shear strains increase, pore water pressure goes up, and the intergranular forces get reduced. As pore water pressures reach a critical level, and the intergranular stresses approach zero, the soil behavior goes from a solid to a viscous liquid state.
Liquefaction is a relevant soil phenomenon for geotechnical design as it may cause local or global failures of foundations and even the collapse of complete structures (Jia (2017) [13]). Among the potential consequences of soil liquefaction, one of the most dangerous ones is lateral spreading. Youd (2018) [12] defined this phenomenon as the horizontal displacement of a soil layer riding on liquefied soil either down a gentle slope or toward a free face like a river channel ( Figure 2). When the underlying soil layer liquefies, the non-liquefied upper soil crust continues moving down until it reaches a new equilibrium position. Prediction of lateral spreading is essential because it can cause damage to the overlying and subsurface infrastructure, and the amount of displacement may influence the design of the infrastructure concerning the decision to, for instance, perform soil improvement in the area affected by this phenomenon (Bray et al. (2017) [10]).

Lateral Spreading Prediction Models
Most of the lateral spreading prediction models are empirical. They use regression procedures to fit equations with field case histories (Hamada et al. (1987) [16], Bartlett and Youd (1995) [2], Youd et al. (2002) [5]). These models take different algebraic forms, and they rely on parameters such as the liquefiable soil's thickness, density, and fines content; earthquake magnitude; and site-to-source distance. Semi-empirical models (Zhang et al. (2004) [17]; Faris et al. (2006) [3]) use other variables, like shear strain ratios and earthquake intensity measures, such as peak surface acceleration. Table  1 shows a list of existing lateral spreading prediction models and their main variables.

Author(s)
Variables Hamada et al. (1987) [16] Ground slope, thickness of the liquefiable layer Bartlett and Youd (1995) [2] Ground slope, thickness of the liquefiable layer, fines content, average grain size, earthquake magnitude, horizontal distance from the site to the seismic energy source. Youd et al. (2002) [5] Ground slope, thickness of the liquefiable layer, fines content, average grain size, earthquake magnitude, horizontal distance from the site to the seismic energy source. Zhang et al. (2004) [17] Ground slope, thickness of the liquefiable layer, shear strain, earthquake magnitude, depth to the liquefiable layer. Faris et al. (2006) [3] Seismic coefficient, earthquake magnitude, horizontal distance from the site to the seismic energy source. Olson and Jhonson (2008) [18] Ground slope, thickness of the liquefiable layer, fines content, average grain size, earthquake magnitude, horizontal distance from the site to the seismic energy source, post-liquefaction undrained shear strength. Zhang et al. (2012) [4] Ground slope, thickness of the liquefiable layer, fines content, average grain size, pseudo-spectral displacement. Gillins and Bartlett (2014) [19] Ground slope, thickness of the liquefiable layer, fines content, average grain size, earthquake magnitude, horizontal distance from the site to the seismic energy source. Pirhadi et al. (2019) [20] Ground slope, thickness of the liquefiable layer, fines content, average grain size, earthquake magnitude, cumulative absolute velocity, peak ground acceleration. Prediction of lateral spreading is essential because it can cause damage to the overlying and subsurface infrastructure, and the amount of displacement may influence the design of the infrastructure concerning the decision to, for instance, perform soil improvement in the area affected by this phenomenon (Bray et al. (2017) [10]).

Lateral Spreading Prediction Models
Most of the lateral spreading prediction models are empirical. They use regression procedures to fit equations with field case histories (Hamada et al. (1987) [16], Bartlett and Youd (1995) [2], Youd et al. (2002) [5]). These models take different algebraic forms, and they rely on parameters such as the liquefiable soil's thickness, density, and fines content; earthquake magnitude; and site-to-source distance. Semi-empirical models (Zhang et al. (2004) [17]; Faris et al. (2006) [3]) use other variables, like shear strain ratios and earthquake intensity measures, such as peak surface acceleration. Table 1 shows a list of existing lateral spreading prediction models and their main variables.

Author(s) Variables
Gillins and Bartlett (2014) [19] Ground slope, thickness of the liquefiable layer, fines content, average grain size, earthquake magnitude, horizontal distance from the site to the seismic energy source. Pirhadi et al. (2019) [20] Ground slope, thickness of the liquefiable layer, fines content, average grain size, earthquake magnitude, cumulative absolute velocity, peak ground acceleration.

Current Models and Large-Magnitude Subduction Earthquakes
Tryon (2014) [7] evaluated six empirical models used in practice (Youd et al. (2002) [5]; Bartlett and Youd (1995) [2], Faris et al. (2006) [3], Zhang et al. (2012) [4], and Zhang et al. (2004) [17]) with three case-histories from the 2010 Maule Mw 8.8 subduction earthquake. He found that site-to-source distances are difficult to define accurately for large subduction zone earthquakes. They can vary significantly between seismic regions, making it difficult to recommend a method for calculating such an "R" value. Figure 4 shows a summary of different distance terms that can be considered: D1 = hypocentral distance, D2 = epicentral distance, D3 = closest distance to high-stress zone, D4 = closet distance to the edge of the fault rupture, D5 = closest distance to the surface projection of the rupture (Joyner Boore distance). In large subduction earthquakes, although there is a small area where the earthquake begins (hypocenter), there are multiple zones on the contact between plates ("patches") where energy is released at different times and with different intensities. Hence, although distances D1, D2, and D3 could be defined, they do not necessarily have a reasonable correlation with the intensity of the ground motion at the site of interest. Additionally, for seismically active countries, such as Chile and Peru, D4 and D5 are very small or even zero. From a design point of view, estimating these distances before an earthquake occurs is very difficult.  [17]) with three case-histories from the 2010 Maule Mw 8.8 subduction earthquake. He found that site-to-source distances are difficult to define accurately for large subduction zone earthquakes. They can vary significantly between seismic regions, making it difficult to recommend a method for calculating such an "R" value. Figure 4 shows a summary of different distance terms that can be considered: D1 = hypocentral distance, D2 = epicentral distance, D3 = closest distance to high-stress zone, D4 = closet distance to the edge of the fault rupture, D5 = closest distance to the surface projection of the rupture (Joyner Boore distance). In large subduction earthquakes, although there is a small area where the earthquake begins (hypocenter), there are multiple zones on the contact between plates ("patches") where energy is released at different times and with different intensities. Hence, although distances D1, D2, and D3 could be defined, they do not necessarily have a reasonable correlation with the intensity of the ground motion at the site of interest. Additionally, for seismically active countries, such as Chile and Peru, D4 and D5 are very small or even zero. From a design point of view, estimating these distances before an earthquake occurs is very difficult. Similarly, Williams (2015) [8] used two case-histories from the 2010 Maule Mw 8.8 earthquake to evaluate the empirical methods developed by Youd et al. (2002) [5] and by Bartlett and Youd (1995) [2], concluding that they are extremely sensitive to the distance term, R, and that the current definition of R for these two methods (the Joyner Boore distance) resulted in predictions that were more than two times the measured values. The semi-empirical models by Zhang et al. (2004) [17] and Faris et al. (2006) [3] also over predicted the displacement but in these cases due to the depth weighting factor of their models. In particular, the empirical model of Zhang et al. (2004) [4] predicted displacements roughly six to eight times larger than the measured displacements. On the other hand, De la Maza et al. (2017) [22] studied one case history (Caleta Lo Rojas) from the 2010 Maule Mw8.8 earthquake. They used the Youd et al. (2002) [5] methodology with different distances, finding that the distance to the zone that bounds 10% of the largest slips resulted in satisfactory values when compared against in-situ post-earthquake measurements.
In this study, we analyzed 13 lateral spread cases from six sites affected by the 2010 Maule Mw 8.8 earthquake, where lateral spreading took place ( Figure 5). Figure 6 shows a comparison between observed and calculated lateral spreading using Youd et al.'s (2002) [5] methodology with three R- Similarly, Williams (2015) [8] used two case-histories from the 2010 Maule Mw 8.8 earthquake to evaluate the empirical methods developed by Youd et al. (2002) [5] and by Bartlett and Youd (1995) [2], concluding that they are extremely sensitive to the distance term, R, and that the current definition of R for these two methods (the Joyner Boore distance) resulted in predictions that were more than two times the measured values. The semi-empirical models by Zhang et al. (2004) [17] and Faris et al. (2006) [3] also over predicted the displacement but in these cases due to the depth weighting factor of their models. In particular, the empirical model of Zhang et al. (2004) [4] predicted displacements roughly six to eight times larger than the measured displacements. On the other hand, De la Maza et al. (2017) [22] studied one case history (Caleta Lo Rojas) from the 2010 Maule Mw8.8 earthquake. They used the Youd et al. (2002) [5] methodology with different distances, finding that the distance to the zone that bounds 10% of the largest slips resulted in satisfactory values when compared against in-situ post-earthquake measurements.
In this study, we analyzed 13 lateral spread cases from six sites affected by the 2010 Maule Mw 8.8 earthquake, where lateral spreading took place ( Figure 5). Figure 6 shows a comparison between observed and calculated lateral spreading using Youd et al.'s (2002) [5] methodology with three R-value definitions. The first one is the original R from Youd et al.'s (2002) [5] methodology, the second one is the distance to the maximum observed coastal uplift, and the third one is the distance used by De la Maza et al. (2017) [22], which is defined as the distance to the zone that bounds 10% of the largest slips. The measured lateral displacements at the selected sites were between 1 and 2 m.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 21 In all cases, the conclusion was similar to those of Tryon (2014) [7], Williams (2015) [8], and De la Maza et al. (2017) [22], namely in that the Youd et al. (2002) [5] model, for large-magnitude subduction earthquakes, overestimates the liquefaction-induced lateral displacements by a factor of more than two. Figure 6c shows, however, that there are a few sites where the predictions were close to the measurements. Those sites were those where the R-value was that of De la Maza et al. (2017) [22] and where the average fines content in the cumulative thickness of the saturated granular layer was less than 5%. This is only an initial observation, and much more case-histories need to be studied before generalizing, or not generalizing, this conclusion.  In all cases, the conclusion was similar to those of Tryon (2014) [7], Williams (2015) [8], and De la Maza et al. (2017) [22], namely in that the Youd et al. (2002) [5] model, for large-magnitude subduction earthquakes, overestimates the liquefaction-induced lateral displacements by a factor of more than two. Figure 6c shows, however, that there are a few sites where the predictions were close to the measurements. Those sites were those where the R-value was that of De la Maza et al. (2017) [22] and where the average fines content in the cumulative thickness of the saturated granular layer was less than 5%. This is only an initial observation, and much more case-histories need to be studied before generalizing, or not generalizing, this conclusion.   In all cases, the conclusion was similar to those of Tryon (2014) [7], Williams (2015) [8], and De la Maza et al. (2017) [22], namely in that the Youd et al. (2002) [5] model, for large-magnitude subduction earthquakes, overestimates the liquefaction-induced lateral displacements by a factor of more than two. Figure 6c shows, however, that there are a few sites where the predictions were close to the measurements. Those sites were those where the R-value was that of De la Maza et al. (2017) [22] and where the average fines content in the cumulative thickness of the saturated granular layer was less than 5%. This is only an initial observation, and much more case-histories need to be studied before generalizing, or not generalizing, this conclusion.

Validation of the Numerical Methodology
The simulations in this study were performed using Cyclic1D, a finite-element program for one-dimensional dynamic site-response analyses (Elgamal et al. (2002) [9]). Cyclic1D uses a multi-yield-surface plasticity constitutive model to simulate the cyclic mobility response mechanism. The constitutive model uses a non-associative flow rule for simulating volumetrically contractive or dilative response due to shear loading.

Validation with Centrifuge Tests
To evaluate the accuracy of the selected numerical methodology, several centrifuge tests were simulated. Table 2 gives a list of the centrifuge experiments that were used, where i = surface inclination, Dr = relative density of the liquefiable layer, H = thickness of the liquefiable layer, a max = maximum horizontal acceleration of the input ground motion, and Dh = residual lateral displacement at the surface. Table 3 shows a list of input parameters of the constitutive model, a suggested range of values recommended in Cyclic1D user's manual (Elgamal et al. (2015) [26]) for saturated granular soil, and the calibrated model parameters for Nevada and Ottawa sand used in this study.  [29] The numerical methodology was previously validated against centrifuge experiments by other researchers (Elgamal et al. (2002) [9]). In this study, we reproduced the experimental results from the projects VELACS (Arulanandan and Scott (1993) [32], Taboada and Dobry (1998) [27]) and LEAP (Kutter et al. (2018) [33], Ziotopoulou (2018) [29]), in addition to the centrifuge tests from Sharp et al. (2003) [28].

General Description of the Centrifuge Tests
The project Verification of Liquefaction Analysis by Centrifuge Studies (VELACS) was a cooperative research effort involving eight universities to study soil liquefaction problems (Arulanandan and Scott (1993) [32]). In that project, a series of dynamic centrifuge tests was performed on a variety of different saturated soil models (Arulanandan and Scott (1993) [32]). In this section, we present the simulation of the centrifuge model 2 of the VELACS project conducted at Rennselaer Polytechnic Institute (RPI) by Taboada and Dobry (1998) [27] and numerically validated by Elgamal et al. (2002) [9]. In this test, the soil profile consisted of submerged 20 cm high (physical model) uniform Nevada sand, of Dr = 40%-45%, inclined by 2 • with respect to the horizontal (more details in Table 4). The experiment was conducted at 50 g centrifugal acceleration. A sketch of the laminar box and the instrumentation is shown in Figure 7. The lateral input shaking applied to the base of the model and its corresponding 5% damped pseudo acceleration response spectra are shown in Figure 8.  Similarly, the Liquefaction Experiments and Analysis Project (LEAP) was a cooperative effort among several universities and research institutes to investigate liquefaction and its effects on geostructures (Kutter et al. (2018) [33]). The data are available on the Network for Earthquake Engineering Simulations (NEES) website (Carey et al. (2017) [34]). In this section, we show the simulation of the centrifuge LEAP project conducted at Rennselaer Polytechnic Institute (RPI) (Kutter et al. (2018) [33]). In this test, the soil layer was 4 m high at the center of the model (prototype dimension), and a uniform medium dense sand (Ottawa F-65) was used. The soil had a relative density of Dr = 65% and a 5° slope (more details in Table 5). A sketch of the laminar box and the instrumentation that was used is shown in Figure 9. The lateral input shaking (motion 2) applied to the base of the model and its corresponding 5% damped pseudo acceleration response spectra are shown in Figure 10.   [33]). In this test, the soil layer was 4 m high at the center of the model (prototype dimension), and a uniform medium dense sand (Ottawa F-65) was used. The soil had a relative density of Dr = 65% and a 5 • slope (more details in Table 5). A sketch of the laminar box and the instrumentation that was used is shown in Figure 9. The lateral input shaking (motion 2) applied to the base of the model and its corresponding 5% damped pseudo acceleration response spectra are shown in Figure 10.  [33]). In this test, the soil layer was 4 m high at the center of the model (prototype dimension), and a uniform medium dense sand (Ottawa F-65) was used. The soil had a relative density of Dr = 65% and a 5° slope (more details in Table 5). A sketch of the laminar box and the instrumentation that was used is shown in Figure 9. The lateral input shaking (motion 2) applied to the base of the model and its corresponding 5% damped pseudo acceleration response spectra are shown in Figure 10.

Model Input Parameters
The multi-yield-surface plasticity constitutive model has 14 parameters that must be calibrated to reproduce the liquefaction phenomenon and the lateral spreading. Table 3 shows the parameters we used for Nevada sand (Dr = 40%) and Ottawa F65 sand (Dr = 65%). The calibration was realized by trial and error to obtain a good fit to the measured response of the centrifuge tests. Peak shear strain, friction angle, and phase angle were estimated from the triaxial tests of VELACS (Arulmoli et al. (1992) [36]) and LEAP (Carey et al. (2017) [34]). The coefficient of lateral pressure was estimated using Jaky's relation (Jaky (1944) [37]). The shear wave velocity was estimated using correlations by Seed and Idriss (1970) [38]. Default values were used for the contraction, dilation, and liquefaction parameters (c1, c2, d1, d2, liq), and the number of yield surfaces (Elgamal et al. (2015) [26]). Finally, additional Rayleigh-type damping was set using the modulus reduction curves and the damping curves proposed by Darandelli (2001) [39].

Comparison of Results
The modeling approach was verified by comparing various dynamic responses under earthquakes loading using two experimental cases: M2-2 and RPI-02, from the VELACS and LEAP

Model Input Parameters
The multi-yield-surface plasticity constitutive model has 14 parameters that must be calibrated to reproduce the liquefaction phenomenon and the lateral spreading. Table 3 shows the parameters we used for Nevada sand (Dr = 40%) and Ottawa F65 sand (Dr = 65%). The calibration was realized by trial and error to obtain a good fit to the measured response of the centrifuge tests. Peak shear strain, friction angle, and phase angle were estimated from the triaxial tests of VELACS (Arulmoli et al. (1992) [36]) and LEAP (Carey et al. (2017) [34]). The coefficient of lateral pressure was estimated using Jaky's relation (Jaky (1944) [37]). The shear wave velocity was estimated using correlations by Seed and Idriss (1970) [38]. Default values were used for the contraction, dilation, and liquefaction parameters (c 1 , c 2 , d 1 , d 2 , liq), and the number of yield surfaces (Elgamal et al. (2015) [26]). Finally, additional Rayleigh-type damping was set using the modulus reduction curves and the damping curves proposed by Darandelli (2001) [39].

Comparison of Results
The modeling approach was verified by comparing various dynamic responses under earthquakes loading using two experimental cases: M2-2 and RPI-02, from the VELACS and LEAP projects, respectively. Figures 11-13 show the good fit between predicted and measured excess pore water pressure (EPP), horizontal accelerations, and lateral displacements from these two tests. In the case of RIP-02 lateral displacement, no measurements were made, so we used the numerical results of Ziotopoulou (2018) [29], Case B, for comparison.   Figure 14 shows the good fit between the predicted and the measured lateral displacements from centrifuge tests of Table 4.  Figure 14 shows the good fit between the predicted and the measured lateral displacements from centrifuge tests of Table 4. Figure 14. Estimated versus measured lateral displacements using Cyclic1D from selected centrifuge tests.

Validation with Historical Cases
Validation using field case histories is more challenging due to the inherent variability of soil and earthquake properties. We simulated the response of two case-histories of lateral spreading

Validation with Historical Cases
Validation using field case histories is more challenging due to the inherent variability of soil and earthquake properties. We simulated the response of two case-histories of lateral spreading during large-magnitude subduction earthquakes: the Lo Rojas Port, in the 2010 Mw8.8 Maule Chile earthquake, and the Matanuska Bridge in the 1964 Mw 9.2 Alaska earthquake. As site effects are considered explicitly in the numerical model approach, only strong motions recorded in rock stations are adequate for this study. For the Chile case, we used the records from the Rapel Station (34.0 • S 71.6 • W) from both horizontal components (PGA NS = 0.20 g and PGA EW = 0.19 g), where PGA = Peak Ground Acceleration. For the Alaska case, we used synthetic records estimated by Mavroeidis et al. (2008) [40] for the city of Anchorage in both horizontal components (PGA NS = 0.25 g and PGA EW = 0.23 g).

Description of Field Conditions
For the Lo Rojas site, there is reliable information on layer stratification, in situ testing, and laboratory tests documented in De la Maza et al. (2017) [22] and Barrueto et al. (2017) [41]. Likewise, the papers of Bartlett and Youd (1992) [42] and Gillins and Bartlett (2014) [19] show test field data from the Matanuska site.
For the Lo Rojas site, the same modeling section selected by De la Maza et al. (2017) [22] was used for validation. This geotechnical model was developed according to the bathymetry and field test information provided by the Ports Department of the Ministry of Public Works. According to Barrueto et al. (2017) [41], the soil profile was composed of four soil units, from top to bottom: poorly graded sand (~10 m thick), clayey sand (~9 m thick), high plasticity clay (~5 m thick), and low plasticity clay (down to 70 m deep before a highly cemented soil). Several laboratory tests were conducted to obtain the mechanical parameters for the soil layers: monotonic triaxial, cyclic triaxial, and column shear test (details in De la Maza et al. (2017) [22] and Barrueto et al. (2017) [41]). Table 6 shows the calibrated parameters used for the Lo Rojas model in this study. Using the pore water pressure-based criteria by Wu et al. 2004 [43], the numerical results show that the upper poorly graded sand liquefied as the excess pore pressure ratio (r u ) reached 1.0 after 20 seconds during the seismic event. For the Matanuska site, the modeling section was developed considering the boreholes taken at the Railroad Bridge Mile Post near the Matanuska River ( Bartlett and Youd (1992) [42]). According to Gillins and Bartlett (2014) [19] the selected soil profile was composed, from top to bottom, of gravel sand (~6 m thick), well-graded gravel (~2 m thick), poorly graded sand (~5 m thick), clayed sand (~9 m thick), and low plasticity clay (down to 70 m deep before a highly cemented soil). Table 7 shows the parameters used in the Matanuska model. Figure 15 shows the soil geotechnical layout at both sites. The Lo Rojas and the Matanuska sites have a stratigraphy of alluvial sediments characterized by upper layers of liquefiable sands underlain by clay. For both sites, default values of Cyclid1D were adopted for the contraction, dilation, and liquefaction parameters (c1, c2, d1, d2, liq), and the number of yield surfaces (Elgamal et al. (2015) [26]). The coefficient of lateral pressure was estimated using Jaky's equation (Jaky (1944) [37]). Mass density was estimated from the Standard Penetration Test (SPT) sounding from De la Maza et al. (2017) [22] for the Lo Rojas case, and the SPT sounding from Gillins and Bartlett (2014) [19] for the Matanuska case. Shear wave velocities were obtained from geophysical field tests from Barrueto et al. (2017) [41] for the Lo Rojas case, and using Mayne (2007) [44] correlations for the Matanuska case. For the clays in the Lo Rojas site, peak shear strain, friction angle, and undrained shear strength were obtained from the geotechnical model of De La Maza et al. (2017) [22] and Barrueto et al. (2017) [41] which were based on monotonic and dynamic triaxial tests. For the Matanuska site, we chose pre-defined Cyclic 1D values based on the information from the boreholes documented in Bartlett and Youd (1992) [42] and Gillins and Bartlett (2014) [19].  [22] selected the Rapel (RAP) ground motion. We used the same criterion to select the station to simulate the lateral spreading in the Lo Rojas site in this study. Figure 16 shows the chosen record (both directions) with a significant duration of approximately 34 s and a PGA of 0.2 g.
The 1964 Alaska Mw 9.2 earthquake caused ground failures and collapsing structures from lateral spreading, and the associated tsunami caused about 130 deaths (Bartlett and Youd (1995) [2]). According to Mavroeidis et al. (2008) [40], no strong motion instruments were operative when that destructive seismic event occurred, so no direct measurement of near field ground motions are available. Consequently, we used a simulated ground motion at the Anchorage site shared by Mavroeidis et al. (2008) [40] to reproduce the lateral spreading case. Figure 17 [22] selected the Rapel (RAP) ground motion. We used the same criterion to select the station to simulate the lateral spreading in the Lo Rojas site in this study. Figure 16 shows the chosen record (both directions) with a significant duration of approximately 34 s and a PGA of 0.2 g.
The 1964 Alaska Mw 9.2 earthquake caused ground failures and collapsing structures from lateral spreading, and the associated tsunami caused about 130 deaths (Bartlett and Youd (1995) [2]). According to Mavroeidis et al. (2008) [40], no strong motion instruments were operative when that destructive seismic event occurred, so no direct measurement of near field ground motions are available. Consequently, we used a simulated ground motion at the Anchorage site shared by Mavroeidis et al. (2008) [40] to reproduce the lateral spreading case. Figure 17 shows the simulated record (both directions) with a significant duration of approximately 152 s and a PGA of 0.25 g.  The 1964 Alaska Mw 9.2 earthquake caused ground failures and collapsing structures from lateral spreading, and the associated tsunami caused about 130 deaths (Bartlett and Youd (1995) [2]).
According to Mavroeidis et al. (2008) [40], no strong motion instruments were operative when that destructive seismic event occurred, so no direct measurement of near field ground motions are available. Consequently, we used a simulated ground motion at the Anchorage site shared by Mavroeidis et al. (2008) [40] to reproduce the lateral spreading case. Figure 17 shows the simulated record (both directions) with a significant duration of approximately 152 s and a PGA of 0.25 g.  Tables 6 and 7 show the constitutive models' parameters used for the numerical runs. They were based on the recommendations by Elgamal (2015) [26] and the available geotechnical data of the sites (De la Maza et al. (2017) [22], Barrueto et al. (2017) [41], Bartlett and Youd (1992) [42] and Gillins and Bartlett (2014) [19]). For each field site, horizontal motions in both directions were analyzed and simulated. The average slope of the Lo Rojas site was based on the geotechnical model of the crosssection by De la Maza et al. (2017) [22]. In Matanuska's case, Youd et al. (2002) [5] and Rauch (1997) [48] reported the ground slope of that location in their database. Figures 18 and 19 show the numerical modeling results of the historical cases in terms of acceleration, excess pore pressure ratio, and lateral displacement time history. The results in Table 8 demonstrate the applicability of the proposed model. The simulated displacements were reasonably close to the measured ones, with a maximum difference of about 30%.    Tables 6 and 7 show the constitutive models' parameters used for the numerical runs. They were based on the recommendations by Elgamal (2015) [26] and the available geotechnical data of the sites (De la Maza et al. (2017) [22], Barrueto et al. (2017) [41], Bartlett and Youd (1992) [42] and Gillins and Bartlett (2014) [19]). For each field site, horizontal motions in both directions were analyzed and simulated. The average slope of the Lo Rojas site was based on the geotechnical model of the cross-section by De la Maza et al. (2017) [22]. In Matanuska's case, Youd et al. (2002) [5] and Rauch (1997) [48] reported the ground slope of that location in their database. Figures 18 and 19 show the numerical modeling results of the historical cases in terms of acceleration, excess pore pressure ratio, and lateral displacement time history. The results in Table 8 demonstrate the applicability of the proposed model. The simulated displacements were reasonably close to the measured ones, with a maximum difference of about 30%.

Simulation Results
Tables 6 and 7 show the constitutive models' parameters used for the numerical runs. They were based on the recommendations by Elgamal (2015) [26] and the available geotechnical data of the sites (De la Maza et al. (2017) [22], Barrueto et al. (2017) [41], Bartlett and Youd (1992) [42] and Gillins and Bartlett (2014) [19]). For each field site, horizontal motions in both directions were analyzed and simulated. The average slope of the Lo Rojas site was based on the geotechnical model of the crosssection by De la Maza et al. (2017) [22]. In Matanuska's case, Youd et al. (2002) [5] and Rauch (1997) [48] reported the ground slope of that location in their database. Figures 18 and 19 show the numerical modeling results of the historical cases in terms of acceleration, excess pore pressure ratio, and lateral displacement time history. The results in Table 8 demonstrate the applicability of the proposed model. The simulated displacements were reasonably close to the measured ones, with a maximum difference of about 30%.

Parametric study using Nonlinear Site Response Analysis
We used the validated numerical methodology in a parametric investigation to study the effects of the key parameters that affect lateral spreading in subduction events. The synthetic soil profiles are an idealization of infinite slopes excited by a range of ground motions ( Figure 20). The range of ground motions and soil profiles was selected using various conditions of lateral spreading cases observed in the field. All input ground motions were recorded during large-magnitude subduction earthquakes.

Soil Profiles and Ground Motions
The synthetic soil profile was selected to represent a gently sloping alluvial sand deposit. The soil profiles were analyzed considering the combination of ground slope inclination (2°), nonliquefiable crust thickness (5 m

Results
A total of 120 nonlinear site response analyses were performed. The numerical models were run in Cylic1D, and the primary purpose was to obtain, for each model, the maximum lateral displacement on the surface and correlate it with different intensity measures (IM). For the parametric study, the ground motions were scaled to 0.11 g, 0.22 g, and 0.33 g of PGA to measure the effect of PGA variability with lateral displacement. Figure 21 and Table 9 show the main results of the parametric analysis.

Parametric study using Nonlinear Site Response Analysis
We used the validated numerical methodology in a parametric investigation to study the effects of the key parameters that affect lateral spreading in subduction events. The synthetic soil profiles are an idealization of infinite slopes excited by a range of ground motions ( Figure 20). The range of ground motions and soil profiles was selected using various conditions of lateral spreading cases observed in the field. All input ground motions were recorded during large-magnitude subduction earthquakes.

Parametric study using Nonlinear Site Response Analysis
We used the validated numerical methodology in a parametric investigation to study the effects of the key parameters that affect lateral spreading in subduction events. The synthetic soil profiles are an idealization of infinite slopes excited by a range of ground motions ( Figure 20). The range of ground motions and soil profiles was selected using various conditions of lateral spreading cases observed in the field. All input ground motions were recorded during large-magnitude subduction earthquakes.

Soil Profiles and Ground Motions
The synthetic soil profile was selected to represent a gently sloping alluvial sand deposit. The soil profiles were analyzed considering the combination of ground slope inclination (2°), nonliquefiable crust thickness (5 m), liquefiable layer thickness (10 m), and liquefiable layer SPT resistance (10 and 15 blows/ft). The groundwater level was assumed at the surface in all cases. The base at the bottom of the profile was modeled as rigid bedrock. A

Results
A total of 120 nonlinear site response analyses were performed. The numerical models were run in Cylic1D, and the primary purpose was to obtain, for each model, the maximum lateral displacement on the surface and correlate it with different intensity measures (IM). For the parametric study, the ground motions were scaled to 0.11 g, 0.22 g, and 0.33 g of PGA to measure the effect of

Soil Profiles and Ground Motions
The synthetic soil profile was selected to represent a gently sloping alluvial sand deposit. The soil profiles were analyzed considering the combination of ground slope inclination (2 • ), non-liquefiable crust thickness (5 m), liquefiable layer thickness (10 m), and liquefiable layer SPT resistance (10 and 15 blows/ft). The groundwater level was assumed at the surface in all cases. The base at the bottom of the profile was modeled as rigid bedrock.
A set of 20 ground motions was obtained from different databases (Kit-NET CISMID, SIBERRISK).

Results
A total of 120 nonlinear site response analyses were performed. The numerical models were run in Cylic1D, and the primary purpose was to obtain, for each model, the maximum lateral displacement on the surface and correlate it with different intensity measures (IM). For the parametric study, the ground motions were scaled to 0.11 g, 0.22 g, and 0.33 g of PGA to measure the effect of PGA variability with lateral displacement. Figure 21 and Table 9 show the main results of the parametric analysis. The coefficient of determination (R 2 ), the Pearson correlation coefficient (ρ), and the Spearman rank correlation coefficient (r s ) were used to quantify the correlation between lateral displacements at the surface (LD) and intensity measurements. These correlation coefficients can range from −1 to +1, where −1 means total negative correlation, 0 means no correlation, and +1 means total positive correlation. As Table 9 shows, the most influential parameters were CAV = cumulative absolute velocity (R 2 = 0.87, ρ = 0.93, r s = 0.88), HI = Housner intensity (R 2 = 0.78, ρ = 0.88, r s = 0.79) and SMV = Sustained Maximum Velocity (R 2 = 0.78, ρ = 0.88, r s = 0.79). In this case, the three coefficients are the highest. On the other hand, the least influential parameters were Tp = predominant period (R 2 = 0.01, ρ = −0.05, r s = 0.17) and DRMS = root-mean-square displacement (R 2 = 0.01, ρ = −0.02, r s = 0.06).

Conclusions
Current liquefaction-induced lateral spreading prediction equations exhibit a large margin of error for large-magnitude subduction earthquakes. Based on the results of our parametric study, the main findings are: (1) The numerical methodology used in this study (using Cyclic1D) can properly simulate pore water pressure generation and shear modulus degradation under strong earthquakes. (2) The applicability of the numerical methodology was verified by comparing the simulated responses and the recorded ones in the centrifuge tests of the VELACS and LEAP projects, and those from simulated historical cases. (3) Our parametric study shows that, for the analyzed cases, the residual lateral displacement has a better correlation with CAV = cumulative absolute velocity (R 2 = 0.87, ρ = 0.93, r s = 0.88), HI = Housner intensity (R 2 = 0.78, ρ = 0.88, r s = 0.79) and SMV = sustained maximum velocity (R 2 = 0.78, ρ = 0.88, r s = 0.79). (4) For large-magnitude subduction earthquakes, the use of the distance terms in empirical formulas is still problematic. (5) In future stages of this research, a probabilistic approach that incorporates different sources of uncertainty in seismic loading, soil properties, and soil geometry will be developed. Acknowledgments: The authors would like to thank the anonymous reviewers whose comments greatly enhanced the quality of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.