Numerical Modelling of Seismic Site Response at Large Strains: A Parametric Study

The numerical analysis of seismic site response at large strains should adopt constitutive models able to guarantee not only a correct modelling of stiffness and damping properties but also a compatibility with the shear strength of the materials. The traditional hyperbolic models used in nonlinear analyses are generally calibrated on stiffness and damping curves and therefore does not necessarily match the soil shear strength. An inaccurate modelling of shear strength can lead to unrealistic predictions of the seismic site response with results that are not necessarily conservative: underestimation or overestimation of the computed surface response depends on the difference between the maximum shear stress implied by the adopted hyperbolic nonlinear model and the real soil shear strength. In this paper, over 1900 onedimensional parametric analyses on ideal sand and clay deposits were executed with DEEPSOIL software. A first comparison was undertaken between equivalent linear and nonlinear analyses; then the nonlinear analyses were addressed to study the influence of shear strength as an input parameter on the results of numerical site response analyses. In particular two strategies to take into account the soil shear strength were considered: an adjustment procedure associated to the standard MKZ hyperbolic model and the GQ/H model which allows the shear strength to be explicitly defined as input parameter of the analyses. This parametric study made it possible to define preliminary threshold shear strain values, beyond which it is necessary to execute numerical analyses with more advanced models or procedures, able to capture the real behavior of the soil at large strains. Indicatively above shear strains of 0.1%, traditional nonlinear models neglecting soil strength can provide unrealistic results, with important overestimation of the seismic motion (up to 30% in terms of PGA at the surface).


Introduction
A crucial step in seismic hazard risk reduction is the numerical modelling of local amplification of ground motion due to site conditions (i.e., local geological, geotechnical and morphological characteristics). Site response numerical analyses can be employed: i) in the framework of advanced seismic microzonation studies, generally carried out at urban scale, supporting the urban planning; ii) at building/infrastructure scale in order to define design seismic actions more accurately with respect to simplified approach based on standard spectra generally implemented in the technical codes [1].
A large number of computer programs to perform numerical site response analyses exists; they basically differ in geometry scheme (1D or 2D/3D), domain of analysis (frequency or time domain), numerical scheme, constitutive relationships for cyclic nonlinear and dissipative soil behavior. A list of popular worldwide numerical codes for site response with relevant features can be found in Regnier et al. [2] which reports the results of a recent international benchmark on nonlinear numerical simulations. Two numerical methods of analysis can essentially be chosen to take into account soil nonlinearity: equivalent linear, based on a series of iterative analyses assuming a visco-elastic soil behavior, or a true nonlinear approach, including a broad range of simplified and advanced soil constitutive models [3]. Advanced elasto-plastic models can adequately capture the behavior at large strains including the development of plastic deformations and the shear strength of soils ( [4][5][6][7], among others). However, these advanced models are difficult to calibrate and manage and therefore more simplified models like those belonging to the family of hyperbolic soil models are often used and implemented in well-known nonlinear codes such as DEEPSOIL [8] and DMOD2000 [9]. The use of hyperbolic models involves several issues: the analytical model for the skeleton curve, criteria for loadingunloading-reloading behavior, total stress or effective stress analysis approach, pore water pressure generation/dissipation models. The parameters of hyperbolic models are generally determined by fitting reference curves of normalized shear modulus and damping values as functions of shear strain. In this way the cyclic soil behavior is adequately captured at small-to-medium strains while unrealistic shear stresses (i.e., non-compatible with soil shear strength) can develop at large strains. Performing an analysis without checking that largestrains stresses are compatible with the shear strength of soils can lead to unreasonable soil behavior resulting in overestimation or underestimation of seismic response [10][11][12].
In this paper, over 1900 one-dimensional parametric analyses on ideal sand and clay deposits were executed with DEEPSOIL code, to study the influence of shear strength as an input parameter in numerical site response analyses. The results of numerical analyses carried out including shear strength according to different procedures, a manual correction as proposed by Hashash et al. [3] and the GQ/H model by Groholski et al. [13] explicitly accounting for shear strength, are compared with "standard" hyperbolic nonlinear analyses. The aim of the study is to quantify the effect of neglecting shear strength on site response analyses results and to identify limiting shear strain values, beyond which it is necessary to execute numerical analyses with more advanced models or procedures, able to properly capture the real behavior of the soil at large strains.

Some Remarks on Shear Strength
In 1D site response analyses the soil deposit is subjected to vertical propagating shear waves loading the soil elements in simple shear conditions (Figure 1a, b). Before earthquake, the soil element is subjected to geostatic vertical and horizontal stresses σ' v and σ' h = k 0 σ' v ; in 1D condition (horizontal ground surface and soil layering), vertical and horizontal direction are principal direction of stress state and therefore there are no shear stress on the vertical and horizontal planes (Figure 1a). During S-waves propagation the soil element is subjected to shear stress τ hv on the vertical and horizontal planes under constant normal stresses ( Figure  1b).
The Mohr circle for the initial stress conditions is shown in Figure 1c as dashed circle. During seismic loading, at failure (solid circle in Figure 1c) the maximum shear stress acting on lateral and horizontal planes is (τ hv ) max which can be derived from simple trigonometric considerations [14]: The maximum shear stress reached in the soil element is τ max (Figure 1c): while the shear stress on the failure plane (τ f ) is lower ( Figure 1c): Different values of shear stress can be therefore assumed as limit values, i.e. as shear strength of the soils (see also [15]). In the following, as in 1D site response analyses the maximum computed shear stresses are the stresses acting on horizontal and vertical planes of soil element, we assumed τ lim = (τ hv ) max as shear strength. Please note that assuming as limit value: leads to a severe overestimation of soil strength (Figure 1c). For purely cohesive materials (i.e, clays in total stress analyses) adopting a Tresca failure criterion the limit value of the shear stress is simply the undrained shear strength [15]:

Methodology
Ideal soil deposits characterized by shear wave velocity increasing with the depth have been defined ( Figure 2). The seismic bedrock (characterized by unit weight γ= 22 kN/m 3 and Vs=800 m/s) was located at the following depths H=20, 40, 60 and 80m from ground surface. Two main lithologies were considered: clay with γ=18 kN/m 3 and nonlinear behavior defined by Darendeli [16] curves for PI=30 and OCR=1 and sand characterized by γ=20 kN/m 3 and Seed and Idriss [17] curves. Darendeli curves are pressuredependent being function of the confining stress here computed assuming 0.6 for the at-rest earth pressure coefficient k 0 . Also for sand the dependence of nonlinear curves from confining stress has been roughly considered by assuming the mean or upper bound of the Seed and Idriss curves depending on the thickness of the soil deposit (Table 1). Different Vs profiles were adopted by setting reference values of equivalent shear wave velocity V S,30 or V S,eq representative of subsoil classes B-C-D-E for the Italian technical code. The equivalent shear wave velocity V S,eq is defined by: being H the thickness of soil deposit, h i and V S,i thickness and shear wave velocity of layer i, N the total number of soil layers. For H higher than 30 m, the equivalent velocity is computed by using H=30 in equation (6) thus obtaining the standard V S,30.
In particular V S,30 /V S,eq = 150 m/s and 250 m/s was assumed for clayey soil deposit while 300 m/s and 450 m/s were assigned to sandy soil deposits (Table 1). By combining the different Vs profiles, thickness and lithology of soil deposits, a total of 16 ideal soil deposits were defined and subjected to numerical analyses. Large Strains: A Parametric Study The input motion has been defined in terms of natural accelerograms of increasing level of input energy. The signals were divided into 6 classes, each consisting of 5 recordings, based on the maximum peak acceleration: 1-10 cm/s 2 , 10-50 cm/s 2 , 50-100 cm/s 2 , 100-200 cm/s 2 , 200-400 cm/s 2 , > 400 cm/s 2 . The accelerograms were extracted from the Italian ITACA database (http://itaca.mi.ingv.it/) and, for the higher energy classes, also from the international PEER database (https://ngawest2.berkeley edu /), assuming the following ranges for magnitude and distance: M=4.5-6.9, D=5-50 km. The analyses were performed with the monodimensional code DEEPSOIL 7.0.20 [8]  A total of 16 (deposits) X 30 (input motions) X 4 (modelling approaches) = 1920 site response analyses were therefore carried out in the present work.
As discussed in the previous section, the shear strength of each layer has been assigned according to the Hardin and Drnevich formulation in effective stress for sands (equation 1) and the Tresca criterion for clays (equation 5).
In order to assume a trend of shear strength parameters with depth somewhat compatible with the stiffness gradient, for the sand profiles the friction angle of each layer has been obtained according to the Robertson and Campanella correlation [18]: Where σ' v = effective vertical stress in the midpoint of the layer and q c is computed from stiffness profile assumed for the deposit by the following relation [19]: To obtain static shear strength values coherent with the shear wave velocities for clays, the following relation has been employed [20]: The Figure 3 shows the comparison of the different shear strength formulations discussed in the previous section along the depth for one of the sand profiles here studied (H=80m, V S,30 = 300m/s).
To calibrate the backbone curve, including the shear strength, the procedure proposed by Hashash et al. [3] and the GQ/H model by Groholski et al. [13] have been used in the present work. The two options are quickly described in the following.

MKZ model + Hashash et al. procedure
The procedure allows to adjust the normalized shear modulus reduction curve to consider the shear strength. It is based on the MKZ model having the following hyperbolic backbone curve: where G max = initial shear modulus, τ= shear stress, γ = shear strain. β, s, and γ r are the model parameters to be determined by fitting the normalized shear modulus and damping curves.
The extended unload-reload Masing Rules are them used to model soil behavior under irregular cyclic loading: 1. For initial loading, the stress-strain curve follows the backbone curve: = b cc d (11) in which τ is the shear stress, and Fbb(γ) is the backbone curve function provided by equation (10); 2. If a stress reversal occurs at a point (γ rev , τ rev ), the stressstrain curve follows a path given by: e e af& = b cc H^ ^a f& J 3. If the unloading or reloading curve intersects the backbone curve, it follows the backbone curve until the next stress reversal.
4. If an unloading or reloading curve crosses an unloading or reloading curve from the previous cycle, the stress-strain curve follows that of the previous cycle.
It is well known that Masing rules [21] results in overestimation of damping at large strains. The MRDF-UIUC approach [22] which modify the Masing rules by applying a reduction factor, has been implemented to provide better agreement with damping curves at large strains.
The MKZ-MRDF model is generally calibrated on modulus reduction and damping curve without control of maximum shear stress. It has been observed that while some target modulus reduction curves underestimate the soil strengths, others overestimate the strength. A procedure to adjust the modulus reduction curve to fit with shear stress at large strains (i.e., shear strength) has been proposed by Hashash et al. [3]. The procedure consists of the following steps: 1. Calibrate the backbone to fit the reference stiffness and damping curves (G/G 0 -γ and D-γ) with MKZ-MRDF model; 2. Compute the implied shear strength as the maximum shear stress value calculated using the following equation where G/G 0 is the modulus reduction curve obtained from backbone curve: 3. Compare the implied shear strength obtained from the previous expression, and the real dynamic shear strength of the soil, estimated as 1.1-1.4 of the static shear strength [23].
4. In the case of the implied shear strength is greater/smaller than the real shear strength, modify the backbone curve by manually decreasing/increasing the reference modulus reduction curve at large strain respectively (indicatively for strains higher than 0.1%) 5. Repeat iteratively the procedure from the step 1 until the implied shear strength match the real soil strength.
The example in Figure 4 shows as the shear strength, computed without the illustrated procedure, is significantly overestimated.
GQ/H model In order to capture both the small and large strain soil behavior, a backbone curve should include: i) an initial (or maximum) shear modulus G max at zero strain, ii) a limiting shear stress at large shear strains (i.e., soil shear strength), iii) a rule controlling the nonlinear behavior between the previous two "boundary conditions". In a stress-strain space these two linear boundaries can be joined by a quadratic model producing a continuous curve. Starting from this assumption, Groholski (14) where τ is the shear stress, τ max is the shear stress at failure, γ is the shear strain, γ r is the reference shear strain and θ t is a curve-fitting parameter. Following Hardin and Drnevich [14], the reference stain γ r is defined as: The coefficient θ t is a model-fitting function adjusting the curvature of the stress-strain curve and does not affect the boundary conditions of the model. It is function of the normalized shear strain:  (16) where parameters p through p R are chosen to provide the best fit to normalized shear modulus versus shear strain curves over a defined strain range. The extended unload-reload Masing Rules are used to models soil behavior under irregular cyclic loading. As said before this choice results in overestimation of damping at large strains. The MRDF approach which modify the Masing rules by applying a reduction factor, developed for the MKZ model, can be implemented also in GQ/H to provide better agreement with damping curves at large strains.
The GQ/H-MRDF model is implemented in DEEPSOIL code, where an automatic fitting curve scheme is available; in the code only the shear velocity, unit weight, nonlinear curves, and the shear strength must be assigned by the user. The Figure 5 shows that the calibration of stiffness and damping curves matching the soil shear strength by GQ/H model and Hashash procedure lead to the same results. In the following the results of the nonlinear analyses including shear strength only refer to the GQ/H model. The GQ/H model directly incorporate the shear strength as input parameter of the constitutive model and is therefore preferable to MKZ model coupled with Hashash procedure; this latter model cannot directly represent large-strain shear strength and the iterative curve adjustment is time consuming and somewhat subjective.

Processing of Results: Synthetic Parameters
The results have been processed by executing two comparisons to evaluate the reliability at the large strains of the approaches typically used in local response analyses, the equivalent linear and nonlinear analyses without shear strength inclusion: 1) equivalent linear vs. nonlinear (hyperbolic MKZ + Masing criteria) to assess the performance of the equivalent linear approach 2) nonlinear MRDF-UIUC vs nonlinear MRDF-UIUC with shear strength inclusion to assess the performance of the nonlinear analysis without including the shear strength.
The second comparison also allowed to quantitatively assess the influence of shear strength on nonlinear site response analyses.
In the following the results of the 960 analyses for sandy profiles and 960 for clay profiles (for a total of 1920 analyses) are reported in terms of Peak Ground Acceleration (PGA) profiles with depths and the maximum acceleration computed at surface PGA z=0 . The comparisons were made using the following synthetic parameters, d and ∆, respectively for PGA profile and PGA at the surface. (18) in which A1 and A2 are the output of the two analyses to be compared and N the number of the PGA values along with the depth (a value has been computed each meter therefore we have H values for a deposit characterized by thickness H).
In particular, in the equations, the "reference" analysis A2 is the approach typically used in local response analyses; therefore, it corresponds to the Linear Equivalent in the first comparison and the Nonlinear without shear strength inclusion in the second one. This assumption is sketched in Figure 6.
To investigate the performance of the different approaches at increasing nonlinear levels, synthetic parameters are represented as a function of the highest strain level reached in the maximum shear strains profile. This peak value of shear strain has been computed from the "reference" A2 analysis.

Results
The comparison between the equivalent linear and nonlinear models (A1 = nonlinear analysis, A2 = equivalent linear analysis) is reported for sandy and clayey deposits in Figure 7a and 8a, respectively. For peak strains lower than 0.1% the PGA profiles are virtually coincident and PGA at surface are characterized by small differences (i.e., below 20%). Note that EL approach provides higher PGA values with respect to nonlinear one, being ∆ values smaller than 1; this overestimation of PGA is well-known and extensively reported in the literature ( [24,25] among others).
The overestimation of PGA by EL approach increases with increasing nonlinearity, a shear strain threshold of 0.02% can be identified corresponding to differences in terms of PGA at surface higher than 20% with respect to more accurate nonlinear analyses. Above this threshold also the synthetic parameter d shows an appreciable increasing trend thus indicating that PGA profiles start to become quite different. At large strains the difference in terms of PGA at surface seems to stabilize around 30%-40% while d parameter shows an exponential growth.
For peak shear strains higher than 0.02% a nonlinear approach is therefore preferable; this indication is more restrictive with respect to other literature recommendations: a peak shear strain around 0.1% is generally indicated as upper limit to employ equivalent linear approach [26].
The comparison between the standard nonlinear analysis and the nonlinear analysis with the shear strength inclusion as input parameter (A1 = nonlinear analyses with τ max , A2 = standard nonlinear analyses) is reported for sandy and clayey deposits in Figures 7b and 8b, respectively.
Below a peak shear strain of 0.1% no significant differences do exist in terms of both PGA profiles and PGA at surface. At higher strains, the discrepancy in terms of PGA at surface is more than 20% and reaches about 40% while the d parameter shows an increasing trend with nonlinear level. Note that ∆ is generally lower than 1 for both sandy and clayey deposits thus indicating that the more accurate analyses including shear strength predict lower PGA values. In other words, neglecting the shear strength lead to an overestimation of the PGA at surface. This behavior (i.e., the underestimation or overestimation) is obviously related to nonlinear G/G 0 curves adopted in the nonlinear analysis. More analyses using different set of nonlinear curves are currently in progress.

Conclusions
The accuracy of seismic site response numerical prediction depends primarily on the capability of employed constitutive models to correctly represent nonlinear soil behavior in the whole strain range of interest. The traditional hyperbolic models used in nonlinear site response analyses are generally calibrated on stiffness and damping curves at small-tomedium strains and not necessarily represent the large-strain behavior of soil, i.e. the shear strength. An inaccurate modelling of shear strength can lead to an underestimation or overestimation of the computed surface response.
In this paper over 1900 one-dimensional parametric analyses on ideal sandy and clayey deposits were executed with 1D DEEPSOIL software. Different stiffness profiles were defined for the soil deposits and a large set of natural input motions with increasing PGA were applied to explore a huge strain range (0.0001-2%).
The analyses were performed considering numerical modelling approaches of increasing complexity: linear equivalent approach, nonlinear analyses with the hyperbolic Modified Kondner Zelasko (MKZ) model and Masing rules, nonlinear analyses with MKZ model and Masing rules modified with a reduction factor of damping at large strains (MKZ-MRDF model), nonlinear analyses with modified Masing rules and inclusion of shear strength in the analyses. Regarding the last step (i.e., the nonlinear analyses including shear strength) the simulations have been executed with two models: the GQ/H model, directly incorporate the shear strength as input parameter of the constitutive model, and the MKZ model (not considering explicitly shear strength) coupled with an procedure consisting in a iterative shear modulus decay curve adjustment. The results were processed in terms of PGA profiles with depth and PGA at surface reported as a function of the highest strain level reached in the maximum shear strains profile. A first comparison was undertaken between equivalent linear and nonlinear analyses; a shear strain threshold of 0.02% can be identified corresponding to differences in terms of PGA at surface higher than 20% between the two methods. Above this threshold true nonlinear analyses are recommended to avoid the wellknown PGA overestimation associated to the equivalent linear method.
Then the nonlinear analyses were addressed to study the influence of shear strength as an input parameter on the results of numerical site response analyses. The parametric study highlighted that above shear strains of 0.1%, traditional nonlinear models neglecting soil strength can provide unrealistic results, with important overestimation of the seismic motion (up to 30% in terms of PGA at the surface).
Literature indications [26] provide that the traditional equivalent approach begins to break down above a shear strain of 0.1% and true nonlinear analyses are therefore necessary in this range. The present study show that above 0.1%, i.e. in the case of highly deformable soils and strong seismic motions, it appropriate not only to discard the EL approach but calibrating the nonlinear models taking into account, in addition to the stiffness and damping properties, the shear strength of the soil (Figure 9).
These threshold values can be regarded as preliminary and should be confirmed by additional processing of the results in terms of other ground motion parameters such as spectral accelerations at surface and profiles of shear stresses and strains with depth.