Numerical modeling of site response at large strains with simplified nonlinear models: Application to Lotung seismic array

2.1 Deepsoil computer code

Deepsoil [13] is a one-dimensional site-response analysis program mainly used for nonlinear (NL) time-domain analysis with or without pore pressure generation. The equivalent linear frequency domain analysis is also available as implemented in other popular codes like SHAKE [19] or STRATA [20]. For NL analysis, a multi-degree-of-freedom lumped model is used in which each individual layer is represented by a corresponding mass, NL spring, and a dashpot for viscous damping. The soil response is obtained from a constitutive model that describes the cyclic behavior. To represent the initial stress–strain backbone curve, Deepsoil employs an updated version of the MKZ model proposed by Matasovic [21]; moreover, a GQ/H model [22], explicitly considering soil strength, is available for backbone curve. Deepsoil uses both standard and a modified version of the extended Masing rules [23] to simulate the unloading-reloading branches.

In order to model the soil damping at small strains where damping from hyperbolic models approaches zero values, Deepsoil provides in addition to the two-control frequencies full Rayleigh formulation [24], an extended Rayleigh formulation [25], in which the user must specify four control frequencies, thus extending the frequency range in which the damping variation is small, and a frequency independent damping formulation [13].

The main features of MKZ and GQ/H models employed in the present study are described in Sections 2.2 and 2.3, respectively.

In Deepsoil, semi-empirical PWP models are available to compute the PWP variation during seismic loading. As stated in the introduction, these models are coupled with the NL hyperbolic models working in total stress (“loosely coupled approach” according to Chiaradonna [26] as an alternative to more complex approach based on advanced constitutive effective-stress models (“fully couped approach”). A summary of available models and codes for both loosely and fully coupled methods is reported in the study by Chiaradonna [26] while main features of strain-based loosely coupled models used in the present study are summarized here.

2.2 MKZ model

The Modified Kondner Zelasko model, known as MKZ, is a constitutive model based on a modification of the hyperbolic model first suggested by Kondner [27] that uses a hyperbolic shear stress–strain relationship according to the following equation:

where G ₀ = initial shear modulus, τ = shear stress, γ = shear strain, while β, s, and γ _r are the model parameters to be determined by fitting the normalized shear modulus and damping curves. Moreover, Deepsoil extends the MKZ through the “pressure-dependent hyperbolic model” proposed by Hashash and Park [28] allowing the reference strain confining to be pressure-dependent (i.e., depending on the effective vertical stress).

The extended unload-reload Masing rules [23] are then used to model soil behavior under irregular cyclic loading:

1. For initial loading, the stress–strain curve follows the backbone curve:

   τ = F bb ( γ )

   where F _bb(γ) is the backbone curve function.

2. If a stress reversal occurs at a point (γ _rev, τ _rev), the stress–strain curve follows a path given by:

   (2) τ − τ rev 2 = F bb γ − γ rev 2 .

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 results in overestimation of damping at large strains [29]. To overcome this problem, the modulus reduction and damping curve fitting (MRDF) Pressure-Dependent Hyperbolic model introducing a reduction factor which effectively alters the Masing rules has been implemented in the code. The reduction factor F(γ _m) is calculated as a function of γ _m, the maximum shear strain experienced by the soil at any given time, according to two alternative formulations modulus reduction and damping curve fitting-University of Illinois at Urbana-Champaign and MRDF-Darendeli described in the study by Park and Hashash [13]. Non Masing MRDF-UIUC approach has been employed in this study. The hyperbolic unload-reload equation is modified with the reduction factor, F(γ _m), as follows:

The MKZ model is generally calibrated on modulus reduction and damping curve without control of maximum shear stress, i.e., the soil shear strength is not explicitly considered in the model. A four-step 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. [6] and followed in the present study:

1. Calibrate the backbone to fit the reference stiffness and damping curves (G/G ₀-γ and D-γ) with MKZ model;

2. Compute the implied shear strength as the maximum shear stress value calculated using the following equation where G(γ)/G ₀ is the modulus reduction curve obtained from backbone curve:

   (4) τ = ρ · V s 2 · G ( γ ) G 0 · γ .

3. Compare the implied shear strength obtained from the previous expression and the real shear strength of the soil, estimated in dynamic conditions as 1.1–1.4 of the static shear strength [30].

4. In the case of the implied shear strength greater/smaller than the real shear strength, manually change backbone curve parameters to decrease/increase, respectively, the G(γ)/G ₀ curve at large strain (indicatively for strains higher than 0.1%).

The procedure is iteratively repeated until the implied shear strength match the real soil strength.

The limiting values of shear stresses in soils can be estimated assuming simple shear conditions in the soil element subjected to vertical propagating shear waves. Before earthquake, the soil element is subjected to geostatic vertical and horizontal stresses σ v ′ and σ h ′ = k 0 · σ v ′ being k ₀ coefficient of earth pressure at rest; in 1D condition (horizontal ground surface and soil layering), vertical and horizontal directions are principal directions of stress state and therefore, there are no shear stress on the vertical and horizontal planes. During S-waves propagation, the soil element is subjected to shear stress τ _hv on the vertical and horizonal planes under constant normal stresses. Assuming a Mohr-Coulomb failure criterion, when the soil comes to failure, the maximum shear stress acting on vertical and horizontal planes (τ _hv)_max that cannot be exceeded during 1D site response analyses, is [31,32]:

with c′ and φ′ cohesion and friction angle of soil in terms of effective stresses.

2.3 GQ/H model

Groholski et al. [22] developed a new general/hyperbolic model characterized by the following backbone curve:

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. This model allows to fit both an initial (or maximum) shear modulus G ₀ at zero strain and a limiting shear stress at large shear strains (i.e., soil shear strength), assuming a quadratic model controlling the NL behavior between these two “boundary conditions.” Following Hardin and Drnevich [31], 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 a function of the shear strain according to the following expression:

where parameters θ ₁−θ ₅ are chosen to provide the best fit to normalized shear modulus vs shear strain curve over the whole strain range.

As for traditional hyperbolic models like MKZ, the extended unload-reload Masing rules as well as MRDF approach, which modify the Masing rules by applying a reduction factor, can be used in conjunction with GQ/H model. An automatic fitting curve scheme is available in Deepsoil and only the shear velocity, unit weight, and the shear strength must be assigned by the user. The shear strength can be estimated as (τ _hv)_max according to the formula reported in the previous section as a function of soil resistance parameters.

2.4 PWP models

As stated in the introduction, the PWP variation induced by seismic loading can be predicted by semi-empirical PWP models coupled with hyperbolic models working in total stress. Pore pressure generation models can generally be categorized into stress-based, strain-based, and energy-based models [6]; the first models are based on the results of cyclic stress-controlled tests while further studies showed a better correlation between PWP build-up and the shear strain level or the energy dissipated within the soil deposit [6,26]. In Deepsoil code, the stress-based model developed by Park and Ahn [33] for sands has been implemented in the code while for strain-based relationships, the Dobry/Matasovic [21] and Matasovic/Vucetic [34] models, for sands and clays, respectively, are available. Moreover, the user can apply the energy-based model GMP [35] for sandy soils and energy-based Generalized model [35,36] for any soil.

In the present study, strain-based relationships have been employed for PWP build-up; these models are briefly described here while more details, including the formulation of stress- and energy-based models, can be found in the study by Park and Hashash [13].

The Dobry/Matasovic model has been originally developed by Dobry et al. [37] and later modified by other authors [21,38,39]. The normalized excess pore pressure (r _u = Δu/ σ vo ′ ) is generated using the following equation:

where N _eq is the equivalent number of cycles; for uniform strain cycles, the equivalent number of cycles is the same as the number of loading cycles while for irregular strain cycles, N _eq is calculated incremented by 0.5 at each strain reversals; γ _c is the current reversal shear strain; γ _tvp is the volumetric threshold: the excess PWP is different from zero only if the shear strain level is greater than the threshold shear strain value; this threshold value is generally determined by means of cyclic or dynamic laboratory tests;

f is a parameter introduced by Vucetic [40,41] in order to simulate the cyclic loading in multiple dimensions: f is equal to 1 in 1D motion and 2 for 2D motion;

p, s, F, and curve are fitting parameters and can be obtained by calibration procedures based on undrained cyclic strain-controlled laboratory tests [21]. In the absence of laboratory data, for practical purposes, p can be assumed equal to 1 while Carlton [42] proposed empirical correlations for the curve fitting parameters F and s for sands:

where V _s is the shear wave velocity in m/s, and FC is the percentage of fines content. Mei and Olson [43] developed further correlation for parameter F using 123 cyclic shear test results compiled from literature referring to sub-angular to sub-rounded clean sands. Two soil index properties, relative density (Dr), and uniformity coefficient (Cu) are used in the correlation.

For clayey soils, Matasović [39] proposed the following equation for the normalized excess PWP:

where N _eq and γ _tvp are the same parameters described in the Dobry/Matasovic model for sands, while the s, r, A, B, C, and D are fitting parameters that can be obtained from the empirical correlations proposed by Carlton [42] as a function of plasticity index (PI) and overconsolidation ratio (OCR) of the clay [13].

The degradation of shear strength and shear stiffness due to the generation of excess PWP during dynamic loading can be considered in the hyperbolic model (MKZ or GQ/H) by introducing the following two degradation indices [21]:

where δ _G is the shear modulus degradation factor, δ _τ is the shear stress degradation factor, r _u is the normalized excess PWP, and ν is a curve-fitting parameter to better model the degradation of shear strength with excess pore pressure generation. These degradation parameter formulations are implemented for all soil models except the Matasovic and Vucetic models for clays where the degradation parameters are defined [21] as follows:

Based on Terzaghi 1D consolidation theory, the code allows to also include PWP dissipation simultaneously to generate during ground shaking.

4.1 Earthquake data

In the framework of a LSST program at a site near Lotung, a seismically active region in northeast Taiwan, the construction of ¼ and ½ scale models of a nuclear plant containment structure was carried out for soil-structure interaction research. The installed instrumentation (Figure 3) involves three linear surface arrays (arms 1, 2, and 3) and two downhole arrays (DHA and DHB). Down-hole array recordings are extremely useful to investigate soil amplification phenomena and to validate constitutive models and numerical modelling techniques for site response analyses [50]. Since its establishment in 1985, the Lotung site has been extensively studied by a large number of authors ([51,52,53,54,55,56,57] among others).

Figure 3

Lotung array: (a) Plain view and (b) schematic cross-section of the site with the deployment of instruments at the ground surface and at the depths of 6, 11, 17, and 47 m.

The location of the downhole arrays, with respect to the edge of the ¼ scale model, is illustrated in Figure 3b. Extensive instrumentations was deployed to record seismic structural and ground response and to monitor soil PWP build-up. More specifically, the arrays were equipped with three-component accelerometers located at the ground surface (FA1-5) and approximately at depths of 6, 11, 17, and 47 m; in addition, a total of 27 pore pressure transducers have been installed. In this study, the free-field site response was studied using the downhole array DHB considering the recordings at borehole sensor located at 47 m depth as input motion.

During a 6-year operation from 1985 to 1990, 30 earthquakes of magnitude comprised between 4.0 and 6.5 triggered this array. In this study, in order to explore different levels of shaking, two events were considered: a weak event (LSST8) and one strong event (LSST7) to investigate small and high strain soil behavior, respectively. Events 12 and 16 (LSST12-16) were also used to validate the semi-empirical PWP models comparing numerical results with the pore water generation recorded at 6 m depth. As matter of fact, to the Authors knowledge, PWP recordings are available in the literature only for these two events.

The records of the transducers installed at 6 m are illustrated in Figure 4, showing a corresponding maximum normalized excess PWP r _u of about 0.17 for event #12 and 0.25 for event #16.

Figure 4

PWP build-up recorded during event LSST12 and LSST16 (after Zeghal et al., [65]).

However, for these events (LSST12-16), the registration of the acceleration at 47 m depth is not available in the downhole array DHB; therefore, the 47 m depth recordings at downhole array DHA were employed as input in these analyses. Main seismological and ground motion parameters of the events selected for the analyses are listed in Table 2.

4.2 Geophysical and geotechnical data

Lotung geology consists of recent alluvium and Pleistocene materials approximately 40–50 and 350 m thick, respectively, overlying Miocene basement. Upper layers comprise about 30 m of predominantly silty sands and sandy silts with some gravels followed by 20 m of interlayered sandy silts, clayey silts, and silty clays [58]. The water table is at about 1 m depth from the ground surface.

As said before, several authors have studied this site, the shear wave velocity profiles measured and adopted by [52,55,59,60,61] are summarized in Figure 5a.

Figure 5

(a) V _s profiles from various studies and (b) numerical viscoelastic linear amplification function.

In this research, the soil profile and parameters assumed in the geotechnical model are referred from the study by Lee et al. [55]. The numerical viscoelastic linear amplification function between surface and the depth of 47 m computed assuming this V _s profile and a uniform damping of 3% is reported in Figure 5b; the obtained fundamental frequency f ₀ = 1.4 Hz is in good agreement with the experimental one in the linear range [53,62], thus confirming the reliability of the V _s profile.

In Figure 6, normalized shear stiffness modulus and damping curves obtained from laboratory testing reported in the study by Tsai [58] and EPRI [63] are illustrated.

Figure 6

Normalized modulus and damping curves for the Lotung site from various literature studies.

In particular, EPRI [63] reports resonant column and torsional shear tests conducted on eight “undisturbed” samples of silty sands and silts at the University of Texas at Austin. The range of experimental curves represented as a grey region in Figure 6 suggests dynamic soil properties close to sands and low plasticity clays as expected being the plastic index of silt of about 7–8% (see the Seed curves [47] for reference). The researchers working on Lotung site have executed numerical analyses assuming different sets of curves. Chang et al. [59] performed equivalent linear analyses by two sets of curves assumed, respectively, for 0–6 and 6–17 m depth ranges, more recently used also by Lee et al. [55] in NL simulations of the Lotung site (Figure 6). Elgamal et al. [64] also computed three material curves from the in situ recorded seismic response for 0–6, 6–11 and 11–17 m layers. Since data were only available to a depth of 17 m, from 17 to 47 m were assumed the same curves as that at 11–17 m depth. This set of curves was also adopted by Borja et al. [52] and Stewart [61] in their seismic response verification studies. As may be observed, considering the stiffness modulus curves, the set proposed by Elgamal et al. [64] fall within the range proposed for sandy soils by Seed and Idriss [47], whereas the Chang et al. [59] curves, also employed by Lee et al. [55], fall along with the lower bound. Analogously, damping curves by Zeghal et al. [65] fall along the upper bound of the Seed and Idriss range, while curves by Chang et al. [59] exhibit much lower damping. As shown in Figure 6, a significant variation in the stiffness and damping curves is evident, confirming a large uncertainty in setting NL cyclic properties of soils.

4.3 Soil parameters and geotechnical model

The soil parameters assumed in the analyses are summarized in Table 3. As shown in Section 4.2, given the significant variation in terms of stiffness and damping curves, numerical analyses have been performed by using different sets of curves with reference to a weak (LSST8) and a strong event (LSST7) to calibrate the NL models in Deepsoil.

The curves employed for the NL calibration involved both a unique set for the whole deposit consisting of Seed & Idriss Lower or Seed & Idriss Mean (set curves 1 and 2 in Table 4) and two sets of curves, respectively, for the 0–17 and 17–47 m layers (set curves 3–7). A summary of all curves is reported in Table 4 and Figure 7 where the employed curves are compared with the range defined by the experimental results reported in Figure 6.

Figure 7

Normalized modulus and damping ratio curves for the Lotung site used in the calibration analyses.

The Figures 8 and 9 show the results in terms of computed vs recorded acceleration response spectra at 17, 11, 6 m and at the surface (0 m).

Figure 8

Recorded and numerical acceleration response spectra (ξ = 5%) for the Event LSST8; numerical analyses carried out using motion at 47 m as input (DHB47); East-West component on the left side and North-South component on the right side. Label 1–7 refer to the different sets of curves reported in Table 4.

Figure 9

Recorded and numerical acceleration response spectra (ξ = 5%) for the Event LSST7; numerical analyses carried out using motion at 47 m as input (DHB47); EW comp. on the left side and NS comp. on the right side. Label 1–7 refer to the different sets of curves reported in Table 4.

For the weak event, LSST8, there is no substantial differences in the numerical results of the analyses carried out using the different sets of curves; each set provides a good agreement with the recordings at all depths. On the contrary, the analyses of the strong event, LSST7, highlight two groups of curves providing a substantial different performances: sets 1–2–3 show a poor prediction of ground motion with a significant underestimation of recorded spectral accelerations in the whole period range and at all depths while a satisfactory agreement with recordings is obtained with curves 4–5–6–7. The selection of the best set is not an easy task (Huang et al., [66]) and is somewhat based on subjective expert judgment. A simple objective statistical analysis has been performed to reach a decision. In the following analyses, set 6 that correspond to Suwal et al. – 4b curves for 0–17 m and Suwal et al. – 4c curves for 17–47 m, has been employed as it provides the lower mean residual (Figure 10):

where S _a,num and S _a,rec are the numerical and recorded spectral acceleration, respectively, and N is the number of points at which the response spectrum is defined between 0 and 2 s.

Figure 10

Performance of the different NL curves: residuals between recorded and numerical spectra for events LSST8 and LSST7.

4.4 Influence of shear strength as input parameter in numerical site response analyses

The parameters employed in the analysis are reported in Table 3. The shear strength was computed by the Hardin & Drnevich formulation equation (6) in which the friction angle and cohesion in terms of effective stress (c′ and φ′) have been assigned according to the modelling assumed by Amorosi et al. [57] based on available experimental data of the Lotung site. In Figure 11 is reported, for Event 7, a comparison of the MKZ (no shear strength included in the analyses) and the GQ/H model (accounting for τ _max) to assess the influence of shear strength as an input parameter in numerical analyses. The results indicate a slight improvement of the numerical prediction at all sensors by introducing the shear strength in the analyses only for the EW component; no appreciable differences can be noted for NS component. Results can be interpreted by observing the stress–strain (τ−γ) curves associated with MKZ and GQ/H models for representative soil elements at 7.5, 20.5, 26.5, and 31.5 m (Figure 12). The GQ/H stress–strain curve obviously perfectly matches the shear strength while MKZ overestimates the maximum shear stress especially for strains higher than 1%. In this analysis, the maximum shear strain experienced by soils is lower than 0.2% for EW component and even lower (<0.1%) for NS component. In this strain range, the discrepancies between MKZ and GQ/H stress curves are quite limited; therefore, the explicit inclusion of shear strength in the modelling lead only to a slight improvement of the prediction. It should be noted that the calibration of NL reference curves reported above allowed us to select curves (set #6) providing a satisfactory numerical prediction of event 7 and therefore modelling a response already compatible with the soil behavior at medium-to-high strains. The explicit inclusion of shear strength in the modelling with the GQ/H analyses allowed a further refinement of the prediction with respect to the standard MKZ NL modelling.

Figure 11

Recorded and numerical acceleration response spectra (ξ = 5%) for the Event LSST7; numerical analyses carried out with MKZ (no shear strength included in the analyses) and the GQ/H (accounting for τ _max) models using motion at 47 m as input (DHB47); EW comp. on the left side and NS comp. on the right side.

Figure 12

Normalized shear modulus curves from MKZ and GQ/H models (on the left) and stress–strain relationship (on the right) for the soil elements at 7.5, 20.5, 26.5, and 31.5 m depth. The reference curves are Suwal et al. (2015) 4b and c for layers 0–17 and 17–47 m, respectively.

It is important to note that the accuracy of the numerical prediction obtained in this study is consistent with previous research utilizing more advanced elasto-plastic NL models [57,68,69]. These studies have also shown that the numerical prediction for event 7 is less accurate for the NS component where an underestimation of recorded response is generally observed like in the present study.

4.5 NL analyses with PWP models

NL analyses have finally been performed by coupling to the GQ/H model with two different semi-empirical PWP models (described in Section 2.4): the Dobry/Matasovic model for the silty sand and sand with gravel layers (0–17, 17–23, 29–36 m), and the Matasovic and Vucetic model for the silty clays (23–29 and 36–47 m layers). A summary of the parameters of PWP models is reported in Table 5. The analyses were carried out for LSST12 and LSST16 events, for which pore pressure recordings are available.

In particular, for the Dobry/Matasovic formulation, the f parameter is used to account for loading in multiple dimensions; in this analysis, it has been set equal to 2 to take into correctly the dimensionality of this problem. The parameter p is set to be 1, as suggested for practical purposes. F and s parameters have been computed according to equations (11) and (12) based on the shear wave velocity and estimated fine content of the materials. Regarding the Matasovic and Vucetic models, the fitting parameters have been obtained by the suggested empirical correlations proposed by Carlton [42] as a function of OCR and PI of soils. No information is available for the volumetric threshold γ tvp , it was set equal to 0.01% for both models, which is a typical value for sandy and low plasticity soils [44,45].

Figure 13 shows the comparison between the recorded and numerical response spectra at the different sensors for LSST12 and LSST16 events using GQ/H total stress analysis and GQ/H coupled with PWP models (GQ/H PWP models). Reference is made to NS component; similar considerations can be made for EW component. The addition of the PWP generation models does not significantly affect the seismic response prediction in terms of response spectra with respect to total stress analyses.

Figure 13

Recorded and numerical acceleration response spectra (ξ = 5%) for Events LSST12 (left side) and LSST16 (right side): comparison between GQ/H and GQ/H PWP analyses (NS component).

The pore pressure comparison in terms of normalized excess PWP time history at 6 m depth is finally shown in Figure 14 for both LSST12 and LSST16 events, with reference to the NS component. For both events, the excess PWP is underestimated even if the order of magnitude is captured by the simplified model. A satisfactory prediction can be observed for LSST16: the recorded pore-pressures reach a peak value of 0.25 with a 0.2 residual value after a dissipation; the numerical model shows a residual r _u of about 0.15 with a 25% error with respect to the experimental value. For comparison, for the same event, a residual value of 0.45 at 6 m (almost double with respect to the recorded one) is predicted by Elgamal et al. [70] with a multi-surface plasticity model and u-p solid-fluid formulation. Even if more research for other events is necessary, the simplified NL strategy proposed here show an accuracy comparable with more advanced models also in terms of PWP build-up.

Figure 14

Build-up pore pressure normalized by the effective vertical stress for events LSST12 and LSST16: computed vs recorded time history at 6 m depth.

The low values of excess PWP reached in the analyses for LSST12 and LSSTR16 justify the negligible differences shown by GQ/H and GQ/H PWP analyses in terms of response spectra observed in Figure 13.