Effect of Nonlinear Soil Model on Seismic Response of Slopes Composed of Granular Soil

A series of two-dimensional finite element analyses are performed to simulate the seismic response of slope composed of granular soil. Four sets of input parameters for the nonlinear soil model are used to fit the reference the shear modulus reduction and damping curves, thereby to evaluate the influence of the nonlinear soil model. ,e first set is fitted to the shear modulus reduction curve. ,e second and third sets are fitted simultaneously to both shear modulus reduction and damping curves. ,e final set applied the shear strength adjustment to adequately capture the nonlinear soil response at large strains.,e accuracy of each set of parameters are evaluated through comparison with centrifuge model test measurements. It is observed that the nonlinear soil model has a marginal influence on the acceleration response. On the contrary, the vertical settlement is highly influenced by the nonlinear soil model. ,e discrepancy is shown to increase with an increase in the intensity of the input ground motion. It is demonstrated that the adjustment for the shear strength is important in performing seismic analyses of slopes, which is most often ignored in practice. Based on the results, practical guidelines on how to select the parameters for the nonlinear soil model are provided.


Introduction
e prediction of the slope stability under severe earthquake loading is one of the primary interests in the field of geotechnical engineering. It is particularly relevant for slopes in the vicinity of nuclear power plants, which can severely damage the facility in case of a slope failure. However, the limit equilibrium procedure based pseudo-static approaches are most often used even for such critical scenarios, despite the fact that the intrinsic limitations of this approximation have been well recognized and documented [1][2][3][4][5][6]. Terzaghi stated that a slope could fail, even though the factor of safety calculated from pseudo-static analysis is greater than unity [7]. e critical limitations include the assumptions of rigid-perfectly plastic behavior of soil and the simultaneous mobilization of the shear strength along the failure surface. Additionally, the failure surface is assumed to be independent of the frequency contents and intensity of the ground motion. e selection procedure of pseudo-static coefficient lacks rigorous theoretical basis and may provide an overly conservative prediction [8]. A Pseudo-static analysis was reported to be a poor for slopes composed of soils that lose their shear strength by more than 15% during earthquake shaking [9]. Although the pseudo-static analysis is a simple procedure for estimating the seismic stability of slopes, it cannot simulate the complex dynamic effects of earthquake. Because of the shortcomings of the pseudo-static method, there is a need to employ more advanced numerical tools, such as two-dimensional (2D) or three-dimensional (3D) dynamic continuum models, to study the seismic stability of slopes. Moreover, because large shear strain is likely to develop along the sliding surface due to the static shear stress present in slopes, the nonlinear soil behavior needs to be accounted for.
Makdisi and Seed [10] performed dynamic finite element analyses using the equivalent linear soil model to develop a simplified procedure to estimate the permanent displacement of dams. Simplified and approximate assumptions were implied, which resulted in conservative estimates. Bouckovalas and Papadimitriou [11] performed a series of numerical analyses with the linear viscoelastic soil to explore the effect of excitation frequency and slope geometry on seismic ground motion. An analytical solution was used to verify the accuracy of their proposed numerical model. However, the simulated results were not calibrated with experimental recordings. Rizzitano et al. [12] studied the topographic amplification of ground motion through 2D finite element analysis using linear and equivalent linear soil properties. e comparisons depicted an underestimation of topographic amplification when the nonlinear behavior of soil is unaccounted for. Du et al. [13] performed sensitivity analyses to evaluate the influence of the slope property variability on slope displacement predictions based on Newmark's sliding block method and fully coupled equivalent linear analysis. e input nonlinear soil properties and shear wave velocity were reported to be the major source of uncertainty that significantly affect the calculated displacement. Song et al. [14] performed dynamic analysis of slopes with different inclinations and potential sliding surfaces to investigate the coupling effect of interaction between the soil topography and the soil layer on slope displacement. To validate the numerical model, the calculated slope surface responses were compared with the results of previously performed numerical studies. Lee et al. [15] performed a series of 2D dynamic nonlinear finite difference (FD) analyses to calculate the permanent seismic displacement of dry mountain slopes. Tsai and Lin [16] performed 2D equivalent linear dynamic analyses to predict the slope displacement assuming a flexible sliding mass. Hailemikeal et al. [17] studied the influence of the subsurface structure and topography on slope response. It was reported that the amplification is dependent on the topography and shear wave velocity structure. Pelekis et al. [18] performed 1D and 2D equivalent linear analyses to study the effect of stratigraphy and surface topography on the seismic response of slopes. Song et al. [19] studied the effect of bedrock depth on the seismic displacement of slopes accounting for the effect of soil properties below the slip surface. Ma et al. [20] examined the different earthquakes to distinguish the effects of the topography and material contrast on the slope amplification. Luo et al. [21] performed 2D and 3D numerical simulations to investigate the seismic response of slopes. It was demonstrated that the local morphology has a primary influence on topographic amplification. e analysis of field monitoring data demonstrated that the topographic amplification is not linearly proportional to the slope height. Cho and Rathje [22] performed a series of nonlinear finite element analysis to develop a predictive model to calculate the slope displacement of clay slopes and used these predictive models for displacement hazard curves. Equivalent linear analysis has been widely used in the seismic response analysis. However, the nonlinear analysis can better capture the behavior of soft soil at large strains [23][24][25]. At shear strains greater than 0.4%, equivalent linear analysis may underpredict the motion. To more accurately capture the soil behavior, a nonlinear analysis needs to be used [26]. e extensive literature review highlights that although numerical models have been widely used to evaluate the seismic response of slopes; none of the studies have thoroughly validated their models through comparison with model tests or field recordings. Considering the sensitivity of the simulated response on a number of factors including the nonlinear soil properties and boundary conditions, the outputs from the computational simulations cannot be fully accepted without proper calibration.
In this study, the centrifuge test measurements of a slope composed of granular soil were utilized to validate the numerical model. e outputs compared herein are the spectral acceleration and the vertical settlement. Moreover, the influence of the parameters selected for the nonlinear model on the calculated response is investigated. Specifically, the influence of the shear strength correction is investigated, the influence of which on the seismic slope stability analysis result has not yet been investigated. e influence of the ground motion intensity is also characterized.

Centrifuge Model Test
e measurements from the dynamic beam centrifuge model tests performed at Korean Advanced Institute of Science and Technology (KAIST) were used to calibrate the numerical model. e effective radius of the centrifuge is 5 m and has a maximum capacity of 240 g-tons. Earthquake loading is applied by an in-flight earthquake simulator equipped with an electrohydraulic system. e earthquake simulator can apply a maximum ground acceleration of 0.5 g into a prototype at a centrifugal acceleration of 55 g [27,28]. e equivalent shear beam (ESB) model container, which has been reported to provide a more representative lateral boundary condition of free-field than a rigid walled container [29], was used. e ESB model container was first proposed by Schofield and Zeng [30]. It was built with a stack of light-weight aluminum frames separated by rubber. e model was designed to produce identical deformation and natural frequency as the soil model, the dynamic stiffness of which is controlled with the flexible frictional end walls.
e performance of the ESB model container was further explored in Madabhushi [31] and Zeng and Schofield [32]. Teymur and Madabhushi [33] performed dynamic centrifuge tests on their ESB model container to investigate the boundary effects. To validate the ESB model container constructed at KAIST, Lee et al. [29] performed a series of model tests and compared the measurement with parallel one-dimensional (1D) site response analyses. It was found that the ESB model container recordings compare favorably with the 1D simulations, thus demonstrating that the freefield boundary condition is well represented in the ESB model container. Technical details of the base plate and the container at KAIST experimental facility used in this study can be found in [29]. e schematic of the centrifuge test model is displayed in Figure 1, including the location of accelerometers and laser sensors. e centrifuge model was prepared at 1/55 scale, and the tests were performed at an acceleration of 55 g.

Advances in Civil Engineering
In prototype scale, the slope is 10 m in height and sloped at an angle of 45°. It is underlain by flat soil with a thickness of 24.65 m. e soil used in the model test was the in situ soil extracted from a slope in the vicinity of a nuclear power plant in Korea. e characteristics of the soil are plotted in Figure 2. e soil is classified as a silty sand (SM) with the Unified Soil Classification System. e resonant column tests were performed to determine the dependency of the shear wave velocity (Vs) on confining pressure, as plotted in Figure 3. e results of these measurements were used to develop the Vs profile. e profile for the flat ground section is shown in Figure 4. e shear wave velocity was varied from 90 m/s at the top to 287 m/ s at the bottom of the soil profile. e shear strength was measured from simple shear tests, as shown in Figure 5. Table 1 lists the properties of soil used in the centrifuge model test.

Numerical Model
e commercial finite element code LS-Dyna [34] was used to perform dynamic analyses of the centrifuge model slope presented in the previous section. e computational model is depicted in Figure 6. e fixed condition was applied for the lateral and bottom boundaries. e equal-degree-offreedom (EDOF) constraint, which is typically used for the lateral boundaries of soil profiles to simulate free-field conditions, was not applied because of the differences in the height of the boundaries. e width of the soil model was selected based on a sensitivity study such that the waves reflected at the lateral boundaries do not influence the seismic response of the slope. e four-node plain strain elements were used for the soil.
e confining pressure dependency of Vs was accounted for in the model. e height of the soil element was set to 0.5 m. It is smaller than λ/8 recommended by Kuhlemeyer and Lysmer [35], where λ is the wavelength of the maximum frequency of interest. e maximum frequency was set to 40 Hz. e hysteretic model (MAT_079) was used to simulate the nonlinear soil response. It is a nested surface plasticity model capable of defining up to 10 yield surfaces. e model is composed of a series of parallel elastic-perfectly plastic materials to produce a nonlinear shear-stress curve [36][37][38]. e parameters for the model were selected to fit the shear modulus reduction and damping curves of Darendeli [39] at the mid depth of each soil layer, the details of which are presented in the following section. In development of the nonlinear curves for each layer using the formulation of Darendeli [39], the overconsolidation ratio and plasticity index were set to 1.0 and 0, respectively. e number of cycles and frequency were set to 10 and 1 Hz, respectively. e small-strain damping ratios of the soil layers were also selected from the Darendeli [39] curves. e small-strain damping was modeled using the frequency-independent damping formulation [34]. e lower and upper frequency bounds of the assigned damping frequency range were set to 0.1 Hz and 30 Hz, respectively, as proposed in Hashash et al. [38]. e measured motion during the Ofunato earthquake was used as input. e motion was amplitude scaled to three peak ground accelerations (PGAs), which are 0.17 g, 0.3 g, and 0.5 g. e acceleration time history of the motion with PGA � 0.17 g, which is set as the baseline input motion, is shown in Figure 7. It should be noted that in the centrifuge test, only this motion was applied.

Evaluation of Nonlinear Soil Models Used in Numerical Simulations
In a seismic analysis, the earthquake ground motion induces cyclic motion, producing nonlinear hysteretic stress-strain curve. e nonlinear soil response is reported to initiate at very small strains, and therefore capturing this is important in a dynamic simulation. e nonlinear soil behavior is typically represented by the normalized shear modulus reduction and damping ratio curve. e shear modulus reduction curve represents the decrease of the secant shear modulus with an increase in strain. e damping ratio curve plots the increase in the area of the hysteretic curve with an increase in shear strain. e performance of a nonlinear constitutive model is calibrated by comparing the shear modulus reduction and damping curves derived from the nonlinear model with the target curves. e nonlinear model was fitted to the shear modulus reduction and damping curves of Darendeli [39], outlined below using three procedures. e first procedure matches the shear modulus reduction curve, denoted as the modulus reduction fit (MR) model. e second procedure fits both modulus reduction and damping curves; hence, it is termed the modulus reduction and damping fit (MRD) model. e third procedure fits the shear modulus reduction curve in addition to adjusting the curve to match the target shear strength. It is labeled as the shear strength fit (SF) model. is additional adjustment is necessary because the shear modulus reduction curves in their original form are reported to provide unreliable predictions for high levels of shear strains. e generalized quadratic/hyperbolic (GQ/H) constitutive model [40] was used to apply the shear strength correction. Mohr-Coulomb failure criterion was used to calculate the shear strength of sand. e friction angle was     determined through simple shear test, as summarized in Table 1. Figure 8 compares the Darendeli and numerically derived curves calculated with 4 sets of parameters when subjected to an effective vertical stress of 80 kPa. For the MR model, a favorable match is obtained with the target shear modulus reduction curve, whereas the damping is significantly overestimated at strains exceeding 0.01%. For the MRD model, two sets of parameters were selected, denoted as MRD-1 and MRD-2, respectively. Among the two models, the MRD-1 model better fits the shear modulus reduction curve but overestimates the damping ratio at shear strains exceeding 0.1%. In comparison, the MRD-2 model better matches the damping ratio curve but overestimates the shear modulus at shear strains exceeding 0.1%. For the SF model, the shear modulus reduction curve is well fitted at small strains, although it overestimates the modulus at strains exceeding 0.03%. e damping curve is favorably approximated up to a strain of 1.0%. However, at higher strains, the damping is overestimated. It should be noted that the Darendeli curves do not provide reliable predictions for shear strains beyond 1%, because they were developed from resonant column and torsional shear test measurements which can only apply shear strains lower than 1%. Consequently, the target curves should only be matched up to a shear strain of 1%. Figure 9 compares the shear stress plotted against shear strain for four levels of effective vertical stresses. Additionally, the target shear strengths calculated with the Mohr-Coulomb model are provided. e MRF model produces significantly lower shear stresses compared with the target shear strength. Moreover, it reaches ultimate resistance at low shear strains, producing significant overestimation of the damping (Figure 8(b)) whereas the shear strength is underestimated. e comparisons illustrate that MRF should not be used in a slope stability analysis. Although the MRD-1 and MRD-2 curves exhibit similarities when compared in the framework of the shear modulus reduction and damping curves, they are demonstrated to be quite different when compared in the form of the shear stress and shear strain curves. e MRD-1 model produces significantly lower shear strength compared with the MRD-2 model. In addition, the residual between these two curves is also dependent on the effective vertical stress. e SF model obviously provides excellent fit with the target shear strengths at all vertical stresses. is model produces high levels of damping at shear strain exceeding 1%, because the shear stress-strain curve levels off as the shear strength is reached. e MRD-2 curve is revealed to favorably fit with the SF model up to a shear strain of 2%. It is therefore concluded from this demonstration that it is possible to approximate the shear strength adjusted model with a conventional nonlinear model, provided that its input parameters are selected appropriately.

Influence of the Nonlinear Model on the
Seismic Slope Response e calculated responses are compared with the recordings at A5, A12, A14, and A19 of the centrifuge model test, as depicted in Figure 10. It is shown that the MRD and SF models successfully predict the recorded acceleration responses of the slope. e output using the MRF model is not displayed because it failed to converge. Apparently, the low shear strength of the model induced unacceptable levels of shear strains, causing it to diverge. It is revealed that the nonlinear soil model has a marginal influence on the calculated acceleration. Figure 11 compares the vertical settlement calculated at the center of slope (S2), normalized to the slope height. e measured response is shown in a grey line. Due to the large fluctuations observed in the recorded settlement, it is difficult to compare the peak settlements. erefore, the measured settlement was smoothened to capture the median response, indicated by a green line. It is  shown that the soil model has pronounced influence on the calculated settlement. e MRD-1 model produces the lowest estimate of the vertical settlement, most probably because of the underestimation of the shear strength. Both the MRD-2 and SF models provide favorable fits with the recording.
is goodness of fits is achieved due to the similarities in the soil stress-strain response up to a shear strain of 2%, as observed in Figure 9. ese comparisons highlight that it is indeed crucial to capture the wide range of the nonlinear response correctly in a seismic slope analysis. e MRD-1 model produces the lowest acceleration at the crest. However, the calculated accelerations are similar for both the MRD-2 and SF models. Overall, the difference in the PGA profile is not significant for different sets of nonlinear soil parameters. Figure 13 displays the calculated vertical settlement time histories subjected to the baseline motion scaled to 0.3 g and 0.5 g. It is demonstrated that the difference increases with an increase in the input motion intensity. e MRD-1 model produces significantly higher settlement, because of the softer stress-strain response. An underprediction of the shear strength in the nonlinear model would result in an overestimation of the vertical settlement, especially for high-intensity motions producing large shear strains.

Conclusions
e influence of the nonlinear soil model on the calculated slope response is investigated from numerical simulations. A 2D nonlinear FE model was used to perform seismic analyses of slopes. Four sets of input parameters were selected for the nonlinear model to fit the reference of the shear modulus reduction and damping curves. e first set, denoted as the MR model, was selected to only fit the shear modulus reduction curve. e second set, labeled as the MRD model, was selected to fit simultaneously to both shear modulus reduction and damping curves. Distinctively, two variations of the MRD models were used. e MRD-1 model better fits the shear modulus reduction curve, whereas the MRD-2 model produces more favorable match with the damping curve. e final set, referred to as the SF model, applied the shear strength adjustment to fit the measured stress-strain behavior at large strains.

Advances in Civil Engineering
When comparing the shear stress versus strain curves derived with the 4 nonlinear models, the MR model produces a significantly lower estimate of the shear stress at strains exceeding 0.1%. Whereas this strain range probably has a secondary influence on the horizontally layered soil profiles, it is important in a seismic slope analysis where large levels of strains are produced. All other models yield higher shear stresses compared to those obtained from the MR model. e MRD-2 model results in the largest shear stress at strains higher than 2%. e calculated responses using four nonlinear models are compared with the centrifuge test measurements for validation.
e MR model fails to converge, because the significant underestimation of the shear strength produces unacceptably high levels of shear strain. It is therefore recommended not to fit the nonlinear model solely to the shear modulus reduction curve. e MRD-1, MRD-2, and SF models produce favorable predictions of the acceleration response. Whereas the calculated acceleration response spectra using the MRD and SF models are similar, the nonlinear soil model is revealed to have pronounced influence on the calculated settlement.
e MRD-1 model produces the lowest estimate of the vertical settlement, apparently due to the underestimation of the shear strength.
e MRD-2 and SF models provide agreeable fits with the recording. e stiffer MRD-2 model yields similar response with the more rigorous strength adjusted model. e comparisons highlight that it is important to capture the large strain nonlinear soil response. Additionally, this study demonstrates that the shear modulus reduction and damping curves, where the shear strains are represented in log scales, may not provide sufficient information on the fit of the numerical model. e shear stress-shear strain plot is revealed to provide additional data on the performance of the nonlinear model at large strain, which is crucial for a seismic slope stability analysis.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.