# Bulletin of the Seismological Society of America

## Abstract

Nonlinear site‐response analyses are becoming an increasingly important component of simulated ground motions for engineering applications. For regional‐scale problems for which geotechnical data are sparse, the challenge lies in computing site response using a very small number of input parameters. We developed a nonlinear soil model that, using only the shear‐wave velocity profile, captures both the low‐strain stiffness and large‐strain strength of soils and yields reliable predictions of soil response to weak and strong shaking. We here present the formulation of the model and an extensive validation study based on downhole array recordings, with peak ground acceleration (PGA) ranging from 0.01*g* to 0.9*g*. We also show that our model, referred to as hybrid hyperbolic (HH), outperforms existing nonlinear formulations and simplified site‐response analyses widely used in practice for ground motions that induce more than 0.04% of soil strain (roughly equivalent to PGA higher than 0.05*g*). In addition to site‐specific response predictions at sites with limited site characterization, the HH model can help improve site amplification factors of ground‐motion prediction equations (GMPEs) by complementing the empirical data with simulated site‐response analyses for very strong ground shaking, as well as physics‐based ground‐motion simulations, particularly for deeper sedimentary sites with low resonant frequencies.

*Electronic Supplement:*Figures showing goodness‐of‐fit score versus peak ground acceleration (PGA) plots for all the nine KiK‐net stations.

## Introduction

Site response has long been known to play an important role in modifying the amplitude, frequency, and duration of earthquake shaking. For nonliquefiable sites, studies showed that site response frequently amplifies the low‐frequency components of weak‐to‐medium intensity motions and introduces complex patterns of amplification and deamplification for higher intensity motions associated with extensive soil yielding (Beresnev and Wen, 1996; Field *et al.*, 1997; Hartzell, 1998; Hartzell *et al.*, 2004).

Near‐surface site response primarily affects the high‐frequency components of ground shaking (>1 Hz). Because these components are frequently captured by stochastic methods in earthquake simulations (Liu *et al.*, 2006; Graves and Pitarka, 2010; Mai *et al.*, 2010), site‐response analyses have been traditionally performed separately from the source and path components of simulated ground motions. For well‐characterized soil sites, the ground surface motion can be computed using elaborate constitutive soil models (e.g., Seidalinov and Taiebat, 2014). For practical applications, however, for which geotechnical site characterization data are scarce and sparse, the use of elaborate geotechnical models introduces unavoidably large uncertainties in the selection of input parameters. Simplified 1D site‐specific models with few input parameters do exist, but these models achieve simplicity by focusing on a specific strain range of soil behavior. For example, elastic–perfectly plastic models such as the Mohr–Coulomb model (de Coulomb, 1776) or the Drucker–Prager model (Drucker and Prager, 1952) have been developed to match the material strength measured in quasi‐static laboratory tests; whereas the Ramberg–Osgood (Ramberg and Osgood, 1943), hyperbolic (Kondner and Zelasko, 1963), and modified hyperbolic models (Matasovic and Vucetic, 1993) have been developed with emphasis on the low‐to‐medium strain range response of soils to cyclic loading (e.g., on the basis of resonant column test results).

In this article, we present a new 1D total stress analysis model that addresses the limitations of existing simplified site‐specific models to simultaneously capture the low‐strain (stiffness) and the large‐strain (strength) response of soils; and furthermore, it does so using the shear‐wave velocity (*V*_{S}) profile as the only input parameter. Given *V*_{S}, the so‐called hybrid hyperbolic (HH) model captures the fundamental physics of strain‐dependent stiffness and strength by means of empirical correlations that have been validated through laboratory experiments and field tests.

Not limited to a specific strain range, our model is thus applicable for site‐response analyses to low‐, medium‐, and strong‐intensity shaking, including near‐field motions, provided that the site conditions and incident motions can be approximately captured by a 1D model. In addition to site‐specific problems, HH can be used as part of hybrid ground‐motion prediction equations (GMPEs) to numerically extend site amplification factors to long‐return‐period events. In the following sections, we present the formulation of the model and an extensive validation study of its performance using 2756 seismic records from nine KiK‐net strong‐motion stations in Japan. Using this statistically significant database of ground motions, we compare HH’s predictions to the widely used equivalent linear and nonlinear constitutive models for site‐specific response analyses.

## A Hybrid Hyperbolic Stress–Strain Model

### Nonlinear Stress–Strain Models for Site‐Response Analyses: State of the Art

Numerous 1D shear stress–strain models have been proposed in the last 50 years. Among others, the hyperbolic model (also known as the KZ model) originally proposed by Kondner and Zelasko (1963), has been extensively used because its formulation is simple and its parameters *A* and *B* reflect physical material properties (1)in which *τ* is the shear stress, and *γ* is the shear strain. Figure 1 shows the example of a KZ stress–strain curve. In this formulation, the tangent slope of *τ*(*γ*) equals to *A* at *γ*=0, and it asymptotically converges to *B* as *γ*→+∞. Thus, when used to approximate shear stress–strain soil behavior, *A* can be set equal to the initial (or maximum) shear modulus of the soil *G*_{max} and *B* equal to the shear strength *τ*_{f}. Equation (1) can then be rewritten as (2)When calibrated to match both stiffness (*G*_{max}) and strength (*τ*_{f}), however, KZ often lacks the necessary geometric flexibility to capture the intermediate strain range. To address this issue, Matasovic and Vucetic (1993) proposed a modified hyperbolic model (hereafter, MKZ), with two additional curvature parameters *β* and *s*. The functional form of MKZ is shown as follows: (3)It is worth noting that the popular study by Darendeli (2001), which provided the functional form of modulus reduction and damping curves for generic soils as a function of soil properties (such as the plasticity index [PI], the *in situ* overburden stress [], and the overconsolidation ratio [OCR], namely, the ratio between the maximum past overburden pressure and ), is based on MKZ’s functional form, with *β*=1 and *s*=0.9190.

However, although the tangent slope of equation (3) at *γ*=0 is *G*_{max} (the low‐strain soil stiffness), *τ*(*γ*) is unbounded for *γ*→+∞. Because equation (3) does not converge to a finite strength at infinite strain, Matasovic and Vucetic (1993) introduced the soil shear strength through an auxiliary reference strain defined as *γ*_{ref}=*τ*_{f}/*G*_{max}, and equation (3) becomes (4)Recognizing that *τ*_{MKZ} is unbounded, Matasovic and Vucetic (1993) and Darendeli (2001) validated MKZ only in the low‐to‐medium strain range, namely, below 0.5%. (0.5% is approximately the largest strain that can be mobilized in resonant column [RC] tests, which were the kinds of tests used by Darendeli, 2001, to derive empirical correlations for the MKZ parameters; Matasovic and Vucetic, 1993, on the other hand, used data from cyclic direct simple shear [cDSS] tests that can mobilize up to 1% strain to validate their model.)

Medium‐to‐strong ground shaking, however, can induce strains beyond 0.5%. For example, as we will show in more detail later in this article, the 26 September 2003 *M*_{w} 8.3 Hokkaido earthquake induced more than 3% strain at seismic station KSRH10 (in Hamanaka, Hokkaido, Japan). To simulate motions of such intensity, the stress–strain model needs to be able to predict not only the small‐strain soil behavior but also (a) the mobilized shear stress at larger strains (beyond the range of applicability of MKZ) and (b) the soil shear strength; thus imposing a physical upper bound to the maximum ground response to seismic shaking.

Recognizing these issues, several recent studies proposed either entire new stress–strain models or modifications of previous formulations based on MKZ. For example, Hashash *et al.* (2010) proposed to manually adjust the functional form of MKZ beyond 0.1% strain so that it matches the correct soil shear strength. Although this approach does capture the correct soil strength, the stress–strain behavior at intermediate strains is artificially (geometrically) constructed to merge the low and large strain ends, and thus may not reflect the true soil behavior. Stewart and Kwok (2008) and Yee *et al.* (2013) proposed and employed a similar geometry‐based modification to MKZ in the 0.3%–0.5% transition strain zone between low‐strain stiffness and material strength. This method also captures the soil strength, but its stress prediction for larger strains has not been validated against laboratory data. Additionally, by geometric construction, this method is prone to discontinuities (e.g., stress singularity, negative stress, or stress decreasing with strain), especially if the MKZ‐predicted stress exceeds the target shear strength at very small strains (which could be, for example, the case for the more confined, deeper layers of soil profiles).

More recently, Gingery and Elgamal (2013) and Groholski *et al.* (2015) proposed two additional shear‐strength‐incorporated models, but neither study has shown that the models are flexible enough to fit laboratory data over the entire strain range. On the other hand, the strength‐incorporated hybrid model proposed by Hayashi *et al.* (1994) was both geometrically flexible and validated against laboratory stress–strain data; however, they did not provide correlations between the model parameters and physical soil properties. This limitation constrains the usability of their model to cases where laboratory stress–strain data are available and reliable. Then, when such data are not available (as is most frequently the case in site‐response analyses), we need a stress–strain model that is not only geometrically flexible enough but also is previously validated against laboratory data, and formulated on the basis of parameters with physical meanings. This was the motivation for the development, calibration, and validation of the HH model presented next.

### The Formulation of the Hybrid Hyperbolic Model

Motivated by the need to capture the dynamic soil behavior over the entire strain range (from stiffness to strength), and in light of the above‐mentioned limitations of existing models, we developed a nonlinear stress–strain model formulated as a composite of MKZ in the low‐to‐medium strain range and a higher‐order, flexible KZ model (hereafter, FKZ) in the large strain range. As such, the new model, coined HH, takes advantage of the flexibility of MKZ in the low‐to‐medium strain range, gradually transitions into FKZ, and asymptotically converges to the shear strength of soil. The novelty of HH lies not only in the flexibility of its functional form to fit a wide range of laboratory stress–strain data but also in its ability to be calibrated with only shear‐wave velocity (*V*_{S}) information, which can yield realistic stress–strain curves when laboratory data are unavailable.

The functional form of the HH model is (5)in which *τ*_{MKZ}(*γ*) is the MKZ stress (defined in equation 4), *w*(*γ*) is a transition function to be defined later, and *τ*_{FKZ}(*γ*) is the aforementioned new flexible hyperbolic (FKZ) model, intended for modeling large‐strain behaviors: (6)Comparing equation (6) with equation (2) (KZ model), we see that FKZ has two additional parameters *μ* and *d*. FKZ is thus a generalized hyperbolic model that collapses to the original two‐parameter KZ model when *μ*=*d*=1. The physical meaning and evaluation of *μ* and *d* will be presented in the HH Prediction of Large‐Strain Shear Stress and Shear Strength section.

The transition function from the stress–strain response of MKZ to FKZ is defined as (7)in which *γ*_{t} and *a* are two parameters, the transition strain and rate of transition, respectively. The formulation of *w*(*γ*) comes from modifying an S‐shape function *s*(*x*)=1/(1+10^{−x}), and fixing the point where *s*(*x*) starts to deviate from the asymptote (the two coefficients 4.039 and 1.036 in equation 7 come from this process).

Referring to Figure 2, *w*(*γ*)=1 when , and it transitions to *w*(*γ*)=0 at a rate controlled by *a*. This means that *τ*_{HH}=*τ*_{MKZ} when and transitions into *τ*_{FKZ} when *γ*>*γ*_{t}. On the basis of the functional form of HH, *γ*_{t} could be physically interpreted as the strain beyond which MKZ is not considered a reliable representation of soil behavior (∼0.01%–3%). (The hybrid model by Hayashi *et al.*, 1994, also has a transition function e^{−αγ}, but contrary to HH, it deviates from 1 at *γ*=0.) Given the extensive literature available on the development and calibration of MKZ (Matasovic and Vucetic, 1993; Darendeli, 2001; Vardanega and Bolton, 2013), we will next focus on the derivation and calibration of FKZ and the transition function. We will base these derivations on laboratory tests on large‐strain soil behavior, namely, direct simple shear (DSS) and triaxial (TX) tests, which usually measure shear strain from ∼0.05% to 20%.

Combining the functional forms of MKZ and FKZ, each formulated to capture a different range of shear strain, we circumvent the dilemma of flexibility versus simplicity without the parameters losing their clear physical meanings. Table 1 summarizes the nine parameters of the HH model, and Figure 3 shows an example of the hybrid stress and hybrid modulus reduction curve predicted by HH.

### HH Prediction of Large‐Strain Shear Stress and Shear Strength

As mentioned above, laboratory tests that measure low strain behavior (e.g., RC test) and medium‐to‐large strain behavior (e.g., DSS, TX tests) only overlap over a very narrow strain range, in the vicinity of 0.5%. Furthermore, in RC tests, soil behavior is affected by large strain rate (i.e., dynamic) effects, whereas DSS tests are performed quasi‐statically. Although RC and DSS data can be merged by careful consideration of rate effects (Vardanega and Bolton, 2013), it is rare to find RC and DSS test data performed on the same soil, subjected to the same state of stress and stress history. Given the constraints above, we initially validated HH by separately evaluating MKZ and FKZ in their corresponding strain ranges against laboratory experiments.

Figure 4 shows the calibration of FKZ using four sets of high‐quality stress–strain data measured in DSS tests (from Ladd and Edgers, 1972; McCarron *et al.*, 1995). Two of the FKZ parameters can be readily obtained from data: *τ*_{f} is the peak shear stress (FKZ and hence HH assumes no post‐peak stress reduction), and *G*_{max} is calculated as follows: (8)in which *G*^{*} is the initial secant modulus of the DSS test, which is different from the initial tangent modulus *G*_{max} measured in RC or TS tests. We estimate *G*^{*}/*G*_{max} by the following MKZ formula: (9)in which *γ*^{*} is the smallest shear strain in the DSS test dataset, and *γ*_{ref} can be calculated using PI, OCR, and (all provided in the DSS test dataset) from Darendeli’s correlation (equation A5).

The third parameter of FKZ, *μ*, can be calculated as follows: Vardanega and Bolton (2011) related the steepness of the DSS stress–strain curve to OCR and (defined under equation A5), (10)in which *γ*_{M=2} is the strain at which 50% of the shear strength is mobilized. Combining equations (10) and (6), temporarily setting *d* to 1.0, we have(11)in which is in kPa, and *G*_{max}/*τ*_{f} and OCR are dimensionless. This indicates that *μ* controls how fast shear stress accumulates when shear strain increases. Then *μ* and *τ*_{f} represent, respectively, two important soil properties in higher strain range: (a) stress accumulation and (b) shear strength.

With *G*_{max}, *τ*_{f}, and *μ* calculated as above, *d* remains the only free parameter in FKZ. Then we estimate *d* by least‐squares curve fitting. As shown in Figure 4, the fit is very satisfactory. The values of *d* fall in a narrow range (mean = 1.03, standard deviation = 0.12), and within this range, *d* did not show clear correlations with any of the soil properties we considered. Because *d* is the power of *μ* in the FKZ formula, we associated the variation of *d* with either uncertainties not captured in equation (11) or errors in the DSS test data.

To further test the validity of equation (11), we perform a blind stress–strain prediction (fixing *d* at the mean value, 1.03) and compare the results with an independent set of DSS data (from Koutsoftas, 1978), as shown in Figure 5. Although these stress–strain curves have been constructed using only *G*_{max}, *τ*_{f}, PI, OCR, and (all provided in Koutsoftas, 1978), the FKZ prediction compares very well with laboratory data.

When stress–strain laboratory data are available in the large strain range, one can directly obtain two of the FKZ parameters *τ*_{f} and *G*_{max}, whereas *μ* and *d* can be estimated from curve fitting of the same data in lieu of equation (11). On the other hand, in absence of laboratory data, we propose an HH calibration (HHC) procedure that constructs empirical *τ*_{FKZ} and thus *τ*_{HH}. We should highlight here that HHC requires only the shear‐wave velocity (*V*_{S}) profile of the site as input—a situation often encountered by engineers and seismologists performing site‐response analysis, just as in the KiK‐net validation study in the A Comparative Study of 1D Site Response Methodologies section. The detailed formulas of the HHC procedure are presented in the Appendix, and they are also coded into our site‐response analysis code SeismoSoil (see Data and Resources).

We should also note that the data in Figure 4 are from static/quasi‐static tests, with strain rate on the order of ∼1% per hour. This strain rate is four orders of magnitude slower than the typical strain rate of dynamic site‐response analyses. To account for the effects of strain rate in site response (Richardson and Whitman 1963; Ladd and Edgers 1972), Vardanega and Bolton (2013) suggested a simple rate‐effect correction factor *Z* (see the Appendix for formula) as part of the following correction scheme, which we adopted in the simulations presented in the following sections: (12)Applying this correction to FKZ, equation (6) becomes (13)As we can see, only static *G*_{max} and static *τ*_{f} need to be scaled up by *Z*, whereas *μ* and *d* are not affected. This means that we can directly use dynamic *G*_{max} and dynamic *τ*_{f} as HH parameters, without the need to recalibrate *μ* and *d*. The procedures for determining dynamic *G*_{max} and dynamic *τ*_{f} are also listed in the Appendix. In short, in the typical case in which only *V*_{S} profiles are available, *G*_{max} need not be corrected by *Z*, whereas *τ*_{f} should be corrected by *Z*.

### HH Prediction of Hysteretic Soil Damping

Another aspect of the validity of a dynamic stress–strain soil model is its capability to properly capture the hysteretic damping ratio as a function of shear strain. Hysteretic damping ratio (or intrinsic attenuation) is defined as the ratio between the energy absorbed in one loading cycle to the maximum elastic energy stored over the same cycle. Some recent studies such as Phillips and Hashash (2009) and Li and Assimaki (2010) proposed hysteretic rules that can fit stress–strain models (such as MKZ) to actual damping data or empirically derived damping curves (e.g., Seed and Idriss, 1970; Vucetic and Dobry, 1991; Electric Power Research Institute [EPRI], 1993; Darendeli, 2001). Such hysteresis rules are capable of simultaneously matching shear modulus and damping data, thus resolving the issue of damping overestimation at higher strains that has been the major drawback of traditional Masing rules (Masing, 1926).

Still, even when such non‐Masing rules are implemented, the numerical scheme cannot accurately represent the actual damping behavior if the stress–strain model is not geometrically flexible enough to fit the damping data. An example is shown in Figure 6, in which we observe the poor fit of MKZ to damping of sands (measured data from Matasovic and Vucetic, 1993) and clays (design curves from EPRI, 1993) at large strains (>1%). Shown in the same figure is the fit of HH to the same data/curve, which is significantly better than MKZ, particularly in the large strain range (>1%) for which some soils exhibit damping reduction under certain combinations of soil type and overburden stress. This damping reduction has been documented in EPRI (1993) and is attributed by Matasovic and Vucetic (1993) to the dilative behavior of soils at higher strains. With a monotonically increasing damping prediction with strain, MKZ overestimates soil damping for certain soil types at larger strains, which in turn could lead to an underestimation of ground response to very strong input motions (namely, in the context of this study, motions that mobilize shear strains higher than 1%). In contrast, HH has more parametric flexibility in capturing damping reduction, which ensures a satisfactory damping representation in the numerical scheme.

We should note here that, in absence of measured data, we recommend the empirical formulas proposed by Darendeli (2001; which are documented in the Appendix) to estimate the damping ratios and calibrate the HH model parameters.

## A Comparative Study of 1D Site‐Response Methodologies

In this section, we implement the MKZ and the HH models in nonlinear site‐specific response analyses and evaluate their performance by comparing simulation results to surface strong‐motion recordings at nine downhole arrays of the Japanese strong‐motion seismograph network KiK‐net (Aoi *et al.*, 2004; see also Data and Resources). At the same time, we also perform a series of simplified site‐response analyses (linear and equivalent linear), and we quantify the strain beyond which their performance diverges from the nonlinear predictions and from the recorded ground motions. In each case, we use the subsurface shear‐wave velocity profile at each station as input, and we employ the same empirical correlations to evaluate the modulus reduction and damping curves. Before we proceed with the quantitative comparative results, we first provide an overview of the characteristics of the three families of site‐response analysis methods used in this study.

In the linear method, the material properties remain constant during shaking. Our linear method is in the frequency domain with linear viscoelastic material behavior. It has been repeatedly shown that the linear method is not suitable for site‐response analyses to strong ground motions, with the possible exception of hard‐rock sites (National Earthquake Hazards Reduction Program [NEHRP] classes A and B) or very weak motions—peak ground acceleration (PGA) at rock outcrop smaller than ∼0.1*g* (Hartzell *et al.*, 2004; Assimaki and Li, 2012; Kaklamanos *et al.*, 2013).

The equivalent linear method, originally proposed by Seed and Idriss (1970), accounts for material yielding (modulus reduction) and hysteretic attenuation (damping) by iteratively matching the soil modulus and damping to a characteristic strain level. Nonetheless, this method is essentially still a linear method because material properties remain constant throughout an iteration, although the stiffness is reduced and damping is increased compared to the linear method. This method yields satisfactory results for relative stiff sites subjected to intermediate levels of strain (<0.1%) but severely underestimates the high frequencies (>5 Hz) of the ground motion, as will be shown in the following sections.

The nonlinear method is performed in the time domain, in which the material properties are adjusted instantaneously to the strain level and loading path. The nonlinear scheme used in this study was developed by Li and Assimaki (2010) and Assimaki and Li (2012) and has the following features:

It uses a memory‐variable technique proposed by Liu and Archuleta (2006) to simulate the frequency‐independent small‐strain soil damping.

It can incorporate any stress–strain model with a closed‐form expression, including (but not limited to) the MKZ and HH models implemented in this article.

It uses the Li and Assimaki (2010) hysteresis rule that is based on the hysteresis model proposed by Muravskii (2005) and can simultaneously match the

*G*/*G*_{max}and damping curves. This goal is achieved using two sets of stress–strain parameters as input for the nonlinear scheme: the first set describes loading/reloading (obtained from fitting a certain stress–strain model to*G*/*G*_{max}measurements, or directly from empirical correlations that generate the stress–strain model parameters), and the second set describes unloading (obtained from fitting the stress–strain model to damping data). Through the geometric representation of narrower and more realistic hysteresis loops, the Li and Assimaki (2010) rule yields a better fit than the extended Masing rules (Pyke, 1979; Kramer, 1996), which seek a compromised matching of*G*/*G*_{max}and damping curves. It should be pointed out that the Li and Assimaki (2010) hysteresis rule can incorporate any closed‐form monotonic stress–strain model, including (but not limited to) MKZ and HH.

### KiK‐Net Strong‐Motion Stations

The KiK‐net strong‐motion seismograph network in Japan consists of ∼670 stations with a ground surface and a downhole array instrument pair (Aoi *et al.*, 2004). To evaluate the performance of the HH model relative to MKZ and the simplified linear and equivalent linear analyses methods, we identified nine stations for which, according to the taxonomy proposed by Thompson *et al.* (2012), 1D wave propagation is a valid assumption for site‐response analyses. In Thompson’s taxonomy, such an assumption is valid when the theoretical transfer function (TTF; using the Thomson–Haskell method; Thomson, 1950; Haskell, 1953) of a horizontally layered profile is in good agreement with the average empirical transfer function (ETF) at the same site, evaluated as the spectral ratio between the surface and borehole measurements (from weaker motions with PGA_{surface}<0.1*g* only). The site characteristics of the nine stations are shown in Table 2, and their *V*_{S} profiles are shown in Figure 7.

### Soil Profiles by Waveform Inversion

Although shear‐wave velocity profiles are available at the KiK‐net stations from suspension logging and downhole tests performed at the time of installation, previous studies (e.g., Assimaki *et al.*, 2006; Assimaki and Steidl, 2007; Thompson *et al.*, 2009) showed that these profiles are sometimes too coarse to capture the effects of site amplification on the higher frequency components of ground shaking, even for very weak ground motions. To obtain a more detailed description of the soil profile at the selected stations, we employed the waveform inversion technique proposed by Assimaki *et al.* (2006), using the *V*_{S} profiles provided by KiK‐net as initial trial profile and five weak ground motions (PGA_{surface}<0.01*g*) across which we averaged the *V*_{S} results.

Figure 8a shows an example at station FKSH14 of the TTF before and after waveform inversion and compares them with the corresponding ETF. We should note here that the ETF in Figure 8a is computed from all 1697 weak ground motions (PGA_{surface}<0.1*g*) recorded at FKSH14, which is a statistically significant amount to demonstrate the effectiveness of the inversion technique. Figure 8b then compares the *V*_{S} profile before and after waveform inversion. Finally, Figure 9 shows the optimized TTF (after waveform inversion) versus the averaged ETF calculated from all weak ground motions at all nine stations.

### Ground‐Motion Data and Goodness‐of‐Fit Criteria

For the simulations presented in the following sections, we used ground motions with PGA higher than 0.01*g* recorded at the nine KiK‐net stations between the year 2000 and 2013. The detailed number of records from each station is shown in Table 2. For each record, east–west and north–south components were first rotated into *SH* and *SV* components according to the azimuth between the station and the epicenter. Successively, the *SH* borehole and ground surface components were used as input motion and benchmark response, respectively.

The simulated ground surface motions were quantitatively compared with the benchmark, using a goodness‐of‐fit (GoF) scheme synthesized from Anderson (2004). As in the work of Kristeková *et al.* (2009), Olsen and Mayhew (2010), Taborda and Bielak (2013), among others, our GoF gauntlet differs from Anderson’s original scheme in the following three aspects:

We omitted the cross correlation between simulation and recording, due to the sensitivity of this GoF measure to small misalignments between waveforms.

Instead of comparing PGA (peak ground acceleration), PGV (peak ground velocity), and PGD (peak ground displacement), we compare the root mean square (rms) of the acceleration, velocity, and displacement, respectively, namely the dominant rather than peak amplitude of the time series.

The range of the score is changed from [0, 10] in the original scheme (with 10 representing perfect match) to [−10, 10]. In our scheme, 0 represents perfect fit between simulations and recordings, positive values indicate overprediction, and negative values indicate underprediction. The formulas for mapping differences to scores are also changed from an exponential‐decay‐type function to the error function scaled up by 10 (as shown in Table 3), because the error function is directional (i.e., under‐ and overprediction distinguishable) and symmetric around zero.

The physical meanings of each score (*S*_{1},*S*_{2},…,*S*_{9}) and the formulas are summarized in Tables 3 and 4. The average of *S*_{1} through *S*_{9}, and can be calculated for one pair of simulated and measured waveforms band‐pass filtered through frequency band [*B*]. In other words, provides GoF at frequency range [*B*]. The final score is the average over five different frequency bands:(14)

## Results and Discussion

We performed site‐response predictions for all the motions with measured PGA_{surface}>0.01*g* at the nine KiK‐net stations, using the following analyses: (a) linear (LN), (b) equivalent linear (EQ), and (c) nonlinear (NL). For the EQ and NL analyses, we employed two stress–strain relations: MKZ and HH. Because *V*_{S} profiles are the only available soil properties useful for site‐response analyses at the KiK‐net stations, we employ the HHC procedure to obtain MKZ and HH parameters (MKZ parameters are a subset of HH parameters). For the NL analyses, we provide MKZ and HH parameters to the numerical scheme, whereas for the EQ analyses, we construct MKZ and HH modulus reduction curves from their parameters and use them as reference points for 1D interpolation to iteratively determine the strain‐compatible soil properties. Damping values come from Darendeli’s empirical correlations (documented in the Appendix), and both the MKZ and the HH models are fitted to the same values, yielding different sets of input parameters. (As is shown earlier in Fig. 6, HH’s fit to the same damping data is better than MKZ’s.)

Evidently, there are five different model‐analyses pairs (hereafter, methods):

LN: linear method;

EQ

_{MKZ}: equivalent linear method using the MKZ modulus curves and Darendeli’s damping curves;EQ

_{HH}: equivalent linear method using the HH modulus curves and Darendeli’s damping curves;NL

_{MKZ}: nonlinear method using the MKZ modulus and damping parameters; andNL

_{HH}: nonlinear method using the HH modulus and damping parameters.

For each of the five methods above, we then calculated the GoF scores, and , as defined in equation (14). In what follows, we synthesize and present the scores for each method at each station against two metrics: (1) measured surface PGA and (2) maximum shear strain within the soil column (as estimated by NL_{HH}) *γ*_{max}.

### Goodness‐of‐Fit Scores versus PGA

An example of the GoF score for each event plotted versus PGA_{surface} (PGA recorded by seismometers on ground surface) at FKSH14 is shown in Figure 10 for the five methods mentioned above. Each point in the figure corresponds to the response of the profile to one event, each calculated with a different method (indicated by different markers).

The LN scores are shown in Figure 10a,b. Clearly, the LN method significantly overestimates the surface ground motion for all frequency bands, a misfit that increases with increasing PGA. This comes hardly as a surprise because the LN method does not consider material degradation (strain‐compatible or instantaneous), and thus it cannot reproduce the reduced site amplification or deamplification observed in the records for stronger ground motions.

Figure 10a next shows the EQ_{MKZ} and EQ_{HH} scores. The values of all three methods are similar for small PGAs but start to deviate for : by EQ_{MKZ} and EQ_{HH} decreases with PGA, and the EQ_{MKZ} GoF is for the most part lower than the EQ_{HH} fit to the ground surface observations. Also, we observe that for the higher frequencies, both and of EQ_{MKZ} decrease with PGA to almost −5; and and of EQ_{HH} at high PGAs also show a slightly decreasing trend with PGA. This observation confirms that the equivalent linear method can severely underestimate high‐frequency motions, especially when employed with non‐strength‐corrected modulus reduction and damping (MKZ in this case) curves.

Figure 10b then shows the scores of NL_{MKZ} and NL_{HH}. For PGA<0.05*g*, the two NL methods yield slightly better predictions than LN, EQ_{MKZ}, and EQ_{HH}. For PGA>0.05*g*, the advantage of NL_{HH} becomes more apparent: its score remains close to 0 (perfect fit) and appears to be independent of PGA. The score by NL_{MKZ}, however, is worse than that of NL_{HH}, yet better than EQ_{MKZ} and similar to EQ_{HH}. Also, higher frequency (above 5 Hz) scores of NL_{MKZ}‐estimated ground motions still indicate an underprediction trend of the observed high‐frequency content.

Assimaki *et al.* (2008) proposed a threshold rock‐outcrop PGA value of 0.2*g*, beyond which nonlinear analyses should be performed in lieu of linear and equivalent linear methods. Similarly, Kaklamanos *et al.* (2013) reported the threshold observed PGA of 0.1*g*. In the present study, the threshold PGA for FKSH14 is shown to be 0.05*g*, significantly lower than in the two previous studies: in this case, nonlinear site response clearly manifests at medium intensity shaking, namely at levels of PGA lower than previously expected.

We should emphasize here that we do not define the threshold PGA where the GoF score crosses over zero, but rather at the PGA where the score starts to deviate from its linear baseline (i.e., the baseline score for very small PGAs). We base this rationale on the fact that this linear baseline reflects all the factors contributing to the misfit in the linear range (e.g., errors in the small‐strain soil properties, deviation from 1D wave propagation conditions, etc.). As the impact of these factors on the GoF usually remains constant with ground‐motion intensity, the contribution of a specific site‐response method to the GoF starts to manifest when a new trend in the score starts to show. The same rationale applies to the threshold maximum strain (*γ*_{max}), which will be presented later.

Figure 11 plots the GoF score versus PGA_{surface} for all the events at all stations. Each marker is the averaged value of all the values within a PGA bin, representing the overall for that specific PGA level. From this figure, we can observe that the by LN is generally the highest (namely, LN overestimates ground motions), except at station TKCH08. For the other four methods, at stations FKSH11, FKSH14, IBRH13, IWTH27, and KSRH10, values are sensitive to increasing PGA and are quite different from method to method: LN overestimates ground motions; EQ_{MKZ} and NL_{MKZ} increasingly underestimate ground motions; EQ_{HH} shows improvement over EQ_{MKZ} and NL_{MKZ}; and NL_{HH} shows even better performance than EQ_{HH}—qualitatively determined because it gives slightly higher than EQ_{HH} and would be considered a conservative (safer) prediction in engineering design. We also observe that NL_{HH} is rather insensitive to the increase of PGA, which can be a significant advantage, because we can know *a priori* that the quality of high‐PGA simulations by NL_{HH} is similar to the quality of low‐PGA simulations, whereas the performance of EQL_{MKZ} or NL_{MKZ} for low‐PGA motions is hardly indicative of that of stronger shaking (higher PGAs). For the other four stations (IBRH10, IBRH17, IWTH08, and TKCH08), by all the methods (except LN) are rather insensitive to PGA and are closer together to each other. Hence, threshold PGA for these nine stations are different, which range from 0.02*g* to more than 0.5*g*. This is because the level of soil nonlinearity at these four stations is not as significant as at other stations, and PGA alone is not always a good indicator of soil nonlinearity. This will be discussed in detail in the Goodness‐of‐Fit Scores versus Maximum Strain section.

### Goodness‐of‐Fit Scores versus Maximum Strain

If we consider the modulus reduction curves used in both EQ and NL methods, the maximum strain level *γ*_{max} induced by earthquake shaking is a strong indicator of the extent of nonlinearity experienced by the soil column during a given event. For example, a hard‐rock site may experience only low strains, even when subjected to strong ground motions, indicating little (if at all) nonlinearity; whereas a soft site may exhibit rather large strains under moderate shaking, indicating that the site response was strongly nonlinear. To examine the relationship between the intensity of ground motion (quantified by PGA) and the extent of nonlinearity (quantified by *γ*_{max}), Figure 12 plots *γ*_{max} versus PGA_{surface} for all stations and events presented above. For a given station, *γ*_{max}‐PGA is a linear correlation with narrow scatter. One could thus conclude that plotting GoF versus *γ*_{max} instead of PGA would result in some scaled version of the GoF–PGA trend presented in the Goodness‐of‐Fit Scores versus PGA section. However, with a closer look at Figure 12, we observe that for a given PGA, *γ*_{max} can vary over a factor of 100 from station to station (e.g., compare KSRH10 to TKCH08). Therefore, for interstation analysis, it is important to factor in different site conditions using *γ*_{max} instead of PGA.

Figure 13 shows the GoF score () plotted against *γ*_{max} for all nine stations and five site‐response methods. The diversity of site conditions provides a wide range of simulated *γ*_{max}, ranging from 0.0004% to 3%, which in turn enables us to better assess the performance of each method. Each subfigure corresponds to one method, and each marker corresponds to one of the 2756 events. The black solid line is the average score within each *γ*_{max} bin, and the dashed lines are standard deviation bounds.

From Figure 13, we clearly observe that the by LN increases from ∼0 to ∼5 with increasing strain. EQ_{MKZ} and NL_{MKZ} remain close to 0 until *γ*_{max} reaches 0.04%, beyond which the score rapidly decreases, indicating that the methods underpredict the ground surface motions for larger strains. This threshold strain of ∼0.04% is lower than the previously reported 0.1%–0.3% threshold by Kaklamanos *et al.* (2013). EQ_{HH} and NL_{HH} show better performance than EQ_{MKZ} and NL_{MKZ}: not only is close to 0, but also ‐*γ*_{max} is approximately constant over the entire strain range. Because the only difference between the HH‐ and MKZ‐based methods is the stress–strain model (HH vs. MKZ), one can conclude that the difference in the GoF is the result of the better soil behavior representation over the entire strain range—from stiffness to strength—of the HH model.

A more detailed examination of Figure 13 reveals that by NL_{HH} is better than the corresponding by EQ_{HH} at high strains (>0.5%). Specifically, although NL_{HH} slightly overpredicts the ground surface motions, EQ_{HH} underpredicts the observed ground shaking, which could lead to unsafe design in engineering practice. Figure 14a depicts only the salient features of Figure 13, by averaging all the values within each *γ*_{max} bin. One can clearly see the advantages of NL_{HH} over the other four methods in the large strain range, in which their predictions become increasingly poor with increased ground‐motion intensity. The strain level at which this prediction divergence between methods is defined as threshold strain, *γ*_{thres}. Figure 14b,c shows the score for higher frequencies and , for which we see that, even though EQ_{HH} has relative satisfactory performance for (averaged across all frequency bands), it severely underestimates higher frequencies when *γ*_{max}>*γ*_{thres}, whereas NL_{HH} still has the best relative performance among five methods. One should note here that *γ*_{thres} for *f*≥5 Hz is even lower, on the order of 0.015%. These results reinforce as well as quantify the well known limitation of the equivalent linear method (regardless of whether the stress–strain curves used in the iterations are corrected for shear strength or not): the use of the equivalent linear method, especially for frequencies above 5 Hz, might lead to severe underestimation of the ground‐motion intensity even for events too weak to be associated with nonlinear site response in practice.

### Case Analysis of Two Strong Events

We next present a detailed analysis of two strong events to further shed light on the difference between the NL and EQ methods and between HH and MKZ constitutive soil models. We specifically discuss the 11 March 2011 *M*_{w} 9.0 Tohoku‐Oki earthquake mainshock recorded at FKSH11 and the 26 September 2003 *M*_{w} 8.3 Hokkaido earthquake mainshock recorded at KSRH10. The 2003 event corresponds to the largest calculated *γ*_{max} of all stations (3.67%); the 2011 event corresponds to the second largest *γ*_{max} (1.13%) of all stations, and the largest *γ*_{max} at FKSH11. Figure 15 shows the waveforms and Fourier spectra (smoothed using a Konno–Ohmachi smoothing window by Konno and Ohmachi, 1998) of the ground surface recording and simulation results. For both events and in line with our previous analyses, LN overestimates ground‐motion intensities. For the 2003 event, EQ_{MKZ} and NL_{MKZ} severely underestimate ground motions, whereas EQ_{HH} and NL_{HH} predict adequate intensity. Still, EQ_{HH} underestimates frequencies higher than 5 Hz (as shown in the Fourier spectra plot). For the 2011 event, all methods except LN underpredict the ground motion. The Fourier spectra plot, however, shows that EQ_{MKZ} and EQ_{HH} severely overdamp higher frequencies (), whereas NL_{HH} provides relatively satisfactory predictions over the full frequency range. (It is worth noting that a good PGA prediction does not necessarily indicate a good match in all frequencies.)

The reason that NL_{HH} does not underestimate surface ground motions as much as NL_{MKZ} is that HH incorporates shear strength, which yields a stiffer response (more elastic) in the shallow layers and a softer response in the deeper layers than MKZ. Figure 16 shows the stress–strain loops calculated by NL_{MKZ} and NL_{HH} at a shallow and a deep layer, for both 2003 and 2011 events. From Figure 16c,f, we can see that at larger depths, NL_{MKZ} and NL_{HH} follow almost the same stress path, with very similar peak stress and peak strain (much lower than *γ*_{t}). However, at shallow depths (Fig. 16a,d), NL_{MKZ} erroneously produces excessively large strain level (over 50% in Fig. 16a) and severely underpredicts stress (nearly no stress in Fig. 16d), whereas NL_{HH} produces stiffer stress–strain loops with higher peak stress and lower peak strain levels. Thus, it is mainly the shallower layers that contribute to the difference between NL_{MKZ} and NL_{HH}.

## Conclusions

We presented a new stress–strain soil model for 1D total stress site‐response analysis, which we coined the HH model. The stress–strain curve of HH compared very well with data from resonance column tests in the small‐strain range and direct simple shear tests in the medium‐to‐large strain range, indicative of the satisfactory geometric flexibility of the functional form of the model. We also showed that, with as little input as the shear‐wave velocity (*V*_{S}) profile of a site, HH can capture the shear strength of soils and also realistically describe the accumulation of stress as a function of strain (before shear strength is reached), based on empirical correlations between *V*_{S} and HH’s parameters.

Subsequently, we presented an extensive validation and comparative study of HH against a widely employed nonlinear model MKZ, both as part of equivalent linear and nonlinear analyses. For this purpose, we used 2756 downhole array recordings at nine KiK‐net strong‐motion stations, for which we used downhole recordings as input and surface recordings as benchmark. Across all stations and events, our findings are briefly summarized below:

The linear method (LN) overestimated the ground‐motion amplitude in the time and frequency domains increasingly with ground‐motion intensity.

The equivalent linear method with the MKZ model (EQ

_{MKZ}) gave satisfactory predictions below*γ*_{thres}=0.04% (or ), whereas it yields increasingly underestimated ground motions with ground‐motion intensity.The equivalent linear method with the HH model (EQ

_{HH}) performed better than EQ_{MKZ}due to the use of the HH model; it underpredicted stronger ground motions (which induced larger strains) to some degree and significantly underestimated high frequencies (above 5 Hz).The nonlinear method with the MKZ model (NL

_{MKZ}) had similar performance and similar trend as EQ_{MKZ}, only slightly better for*γ*>*γ*_{thres}.The nonlinear method with the HH model (NL

_{HH}) had the most satisfactory performance across all ground‐motion intensities (PGA up to 0.9*g*, strain up to 3.67%) and for both broadband and high frequencies (unlike EQ_{HH}); perhaps most importantly, its GoF appeared insensitive to the increase of ground‐motion intensity.

Through this article, we attempt to further advance our understanding of the predictive capabilities (as well as the limitations) of 1D site‐specific response models: a topic that has become the focus of a series of recently published papers. For example, two recent studies, Kaklamanos *et al.* (2015) and Zalachoris and Rathje (2015), employed validation methodologies similar to this study, although both were based on significantly smaller datasets of events. Their GoF metrics were also different from those used in this article: they calculated the residue of response spectra or amplification function between simulations and recordings, whereas we included additional intensity measures, such as rms acceleration, rms velocity, Arias intensity, and so on, each correlating with different aspects of the destructive potential of seismic shaking. Zalachoris and Rathje (2015) concluded that EQ and NL both underpredicted ground motions for strains larger than 0.4%, even when using stress–strain curves corrected for shear strength, which was not consistent with our results. The source of such a discrepancy could be from (1) their use of the strength‐incorporation scheme by Hashash *et al.* (2010), which, as pointed out earlier, employed engineering judgment rather than a rigorous calibration process, or from (2) their choice of validation sites: among the 11 sites used in that study, only four have been shown by Thompson *et al.* (2012) to conform to 1D wave propagation conditions. Kaklamanos *et al.* (2015), on the other hand, used KiK‐net stations that satisfy the 1D wave propagation conditions and concluded that the nonlinear method (equivalent to NL_{MKZ} in the present study) offered a slight improvement over equivalent linear method (equivalent to EQ_{MKZ} in the present study), which is in general consistent with our findings (see Fig. 14a). They did not, however, use strength‐incorporated models in their work.

Comparing the median GoF score of all methods, we determined that 0.04% is the threshold strain beyond which the performance of all the other four methods deteriorate except NL_{HH}. This value is lower than that reported in Kaklamanos *et al.* (2013) (0.1%) and only slightly lower than in Kaklamanos *et al.* (2015) (0.05%). The validation exercise in this study suggests that the use of a strength‐incorporated stress–strain soil model in a time‐domain nonlinear scheme is important not only for site‐responses analyses of rare (and extreme) events but also for medium‐intensity events (with PGA as low as 0.05*g*).

Notwithstanding the rigor of our analyses, the statistical significance of the dataset that we used could be challenged by the scarcity of very strong ground motion recordings. However, the factors that cause imperfect GoF for weaker ground motions, such as errors in initial velocity and attenuation (or damping) profiles or the angle of wave incidence, would become insignificant for stronger events. Thus, we believe that the GoF scores at higher strains—as few as they may be—are still statistically representative. Because global databases of strong‐motion downhole array recordings are limited, centrifuge model tests could complement the available strong‐motion field recordings with carefully calibrated nonlinear site‐response tests to be used as benchmarks for similar studies in the future.

Last, we would also like to note that the KiK‐net ground‐motion dataset represents a group of sites with profiles generally stiffer than the conditions prevailing in highly seismically active regions in the United States, such as California and central and eastern United States. Consequently, the maximum strain level presented here is most likely a lower bound of the strain that sediments in the Los Angeles basin or the Mississippi embayment would experience for the same shaking. Because our analyses showed that the threshold strain beyond which nonlinear analyses become increasingly important is only 0.04%, we believe that the use of nonlinear analyses with strength‐incorporated models such as HH should be considered instead, or at least in addition to equivalent linear analyses, even for ground motions that have been traditionally considered too weak to cause substantial soil yielding.

## Data and Resources

The KiK‐net ground‐motion data used in this study are available online at http://www.kik.bosai.go.jp/ (last accessed May 2016). The model calibration and site‐response analyses were performed using SeismoSoil, a 1D site‐response analysis code with graphical user interface (GUI), available for download from http://asimaki.caltech.edu/resources/index.html (last accessed May 2016).

## Appendix

This provides a guideline for empirically constructing modulus reduction and damping curves using only the *V*_{S} profile as input. We should remind the reader here that all rules of thumb for the selection of parameters should be reserved for the cases in which there is no information on the site conditions, and that measured material properties, if available, are always preferable.

Stress–strain curves can be constructed in the following steps. This is referred to as the hybrid hyperbolic calibration (HHC) procedure:

Evaluation of , the vertical confining pressure, at a specific depth (A1)in which

*g*is the gravitational acceleration,*ρ*_{j}and*h*_{j}are the mass density and thickness of each soil layer, respectively. The summation should be carried out from the ground surface to the layer of interest.*u*is the water pressure at that specific depth. If no information about the water table is available, then the decision lies upon the engineer or scientist to assume dry or saturated soil conditions. (In the case of the present study, we assume the soils are dry.) The mass density can be evaluated, following Mayne*et al.*(1999), as (A2)in which*ρ*is in g/cm^{3},*V*_{S}is in m/s, and*z*is the depth of the soil layer in meters.Evaluation of vertical preconsolidation stress of soil following Mayne

*et al.*(1998) (A3)in which is in kPa, and*V*_{S}is in m/s.Evaluation of OCR (overconsolidation ratio): (A4)

Construct

*τ*_{MKZ}using equation (4), with*β*=1,*s*=0.9190, and*γ*_{ref}evaluated as follows: (A5)(recommended by Darendeli, 2001), in which is the mean effective confining pressure (unit: atm), to be defined in the next step; , , , and are calibrated by Darendeli (2001); and PI is the plasticity index. With only*V*_{S}information available, we use a rule of thumb to determine PI, which is similar in principle to that of Zalachoris (2014): (A6)The mean effective confining pressure, is the average of three stress components (two horizontal , and one vertical ) (A7)in which

*K*_{0}is the ratio of horizontal and vertical confining pressure and is evaluated as (A8)(following Mayne and Kulhawy, 1982), in which is the effective friction angle of soils and is taken as 30° (without better available information).Evaluation of dynamic shear strength of soil

*τ*_{f}. We choose the formula of undrained shear strength (*s*_{u}), following Ladd (1991), for all soil types with*V*_{S}≤760 m/s, because earthquakes’ loading is imposed so quickly that “even coarse‐grained soils do not have sufficient time to dissipate excess pore‐water pressure, and thus undrained condition applies” (Budhu, 2011, p. 267). For materials with*V*_{S}>760 m/s (rocks or very stiff soils), we use Mohr–Coulomb criterion to determine the shear strength. Then we will apply the rate‐effect correction factor*Z*: (A9)in which is the normal effective stress on the failure plane, calculated as (A10)in which and are the larger and smaller one, respectively, between (vertical confining stress) and (horizontal confining stress), and is 30° as before.*Z*, the rate‐effect correction factor, has the form (following Vardanega and Bolton, 2013) (A11)in which is the strain rate of the cyclic loading (unit: s^{−1}), and 10^{−6}s^{−1}is the strain rate of the static tests, such as direct simple shear (DSS). A typical earthquake frequency is taken as 1 Hz, then a simplistic shear strain rate is 10^{−2}s^{−1}(Vardanega and Bolton, 2013), yielding*Z*=1.20. We use 1.20 for all ground motions.Evaluation of

*G*_{max}, the initial shear modulus, from wave propagation theory, (A12)in which*ρ*is in kg/m^{3}, and*V*_{S}is in m/s. This*G*_{max}from equation (A12) is the dynamic shear modulus, thus rate‐effect correction is not necessary. However, when a*G*_{max}is calculated as the initial slope of a static stress–strain test, it needs to be multiplied by*Z*to be used as a HH model parameter.Evaluation of

*μ*, the following formula (identical to equation 11) is used: (A13)Because the original study by Vardanega and Bolton (2011) that provided equation (10) was for clays and silts, we apply equation (A13) only to soils with . For rocks or very stiff soils (*V*_{S}>760 m/s), we simply use*μ*=1 for lack of more information. Inaccurate as this may seem, the shear strain in such stiff layers seldom reaches*γ*_{t}, thus the FKZ part does not play a role in the simulation for these layers.Finally, use an optimization algorithm to find

*d*,*γ*_{t}, and*a*, so that*τ*_{HH}is continuous and monotonically increasing and the difference (or the area) between*τ*_{MKZ}and*τ*_{FKZ}is minimized (in order for MKZ and FKZ to correspond to similar soil types at smaller strains). The search range for*d*is set to 1.03±(3×0.12) (using the 3*σ*rule), in which 1.03 and 0.12 are the mean and standard deviation of the curve fitting in Figure 4. And the search range for*γ*_{t}is [0.01%, 3%], which is slightly wider than [0.1%, 1%] to ensure MKZ and FKZ can merge. In practice,*γ*_{t}usually takes the strain where*τ*_{MKZ}=*τ*_{FKZ}, and*a*usually takes 100 (quick transition). Admittedly, uncertainties would be introduced in this step, but this is unavoidable because we do not have any actual stress–strain measurements, only a shear‐wave velocity value.Damping–strain curve is calculated following Darendeli (2001), (A14)in which

*G*/*G*_{max}=*τ*_{MKZ}/(*γG*_{max}), and(A15)in which (A16)and*ξ*_{min}(or*D*_{min}, in some other literature) is calculated as (A17)in which , , , , and are calibrated by Darendeli, and*f*is the earthquake frequency, taking a nominal value of 1.0 Hz.

## Acknowledgments

This research was supported by the Southern California Earthquake Center (SCEC; Contribution Number 6019). SCEC is funded by National Science Foundation (NSF) Cooperative Agreement EAR‐1033462 and U.S. Geological Survey (USGS) Cooperative Agreement G12AC20038. We are grateful for the financial support that has made this research possible. We would also like to acknowledge the help of Paul J. Vardanega and Malcolm D. Bolton, who provided us with direct shear test data from the literature and from their experimental work; and Rui Wang, who helped us conceptually refine our soil constitutive model. We finally wish to thank two anonymous reviewers for their constructive comments, which led to significant improvements over our initial submission.

- Manuscript received 21 October 2015.