Prediction of Nonlinear Soil Response
Prediction of Nonlinear Soil Response
Prediction of Nonlinear Soil Response
Abstract Mathematical models of soil nonlinearity in common use and recently developed nonlinear codes are compared to investigate the range of their predictions. We consider equivalent linear formulations with and without frequency-dependent moduli and damping ratios and nonlinear formulations for total and effective stress. Average velocity proles to 150 m depth with midrange National Earthquake Hazards Reduction Program site classications (B, BC, C, D, and E) in the top 30 m are used to compare the response of a wide range of site conditions from rock to soft soil. Nonlinear soil models are compared using the amplication spectrum, calculated as the ratio of surface ground motion to the input motion at the base of the velocity prole. Peak input motions from 0.1g to 0.9g are considered. For site class B, no signicant differences exist between the models considered in this article. For site classes BC and C, differences are small at low input motions (0.1g to 0.2g), but become signicant at higher input levels. For site classes D and E the overdamping of frequencies above about 4 Hz by the equivalent linear solution with frequencyindependent parameters is apparent for the entire range of input motions considered. The equivalent linear formulation with frequency-dependent moduli and damping ratios under damps relative to the nonlinear models considered for site class C with larger input motions and most input levels for site classes D and E. At larger input motions the underdamping for site classes D and E is not as severe as the overdamping with the frequency-independent formulation, but there are still signicant differences in the time domain. A nonlinear formulation is recommended for site classes D and E and for site classes BC and C with input motions greater than a few tenths of the acceleration of gravity. The type of nonlinear formulation to use is driven by considerations of the importance of water content and the availability of laboratory soils data. Our average amplication curves from a nonlinear effective stress formulation compare favorably with observed spectral amplication at class D and E sites in the Seattle area for the 2001 Nisqually earthquake.
Introduction
Since the occurrence of the 1994 Northridge, 1995 Hyogoken-Nanbu (Kobe), and the 2001 Nisqually earthquakes, soil nonlinearity has assumed a prominent role in the analysis of strong ground motion. After the Northridge earthquake, several studies reported widespread nonlinear sediment response (Field et al., 1997, 1998; Trifunac and Todorovska, 1996, 1998; Beresnev et al., 1998; Hartzell, 1998). Previous to this earthquake, clear measurements of soil nonlinearity in strong ground motion were generally limited to the extreme cases of liquefaction. Aguirre and Irikura (1997), Fukushima et al. (2000), and Frankel et al. (2002) have reported varying degrees of soil nonlinearity, dependent on the site conditions, for the Hyogoken-Nanbu and Nisqually earthquakes. Although there has been a general realization among seismologists that nonlinear effects are more common than previously believed, there is uncertainty concerning the appropriate mathematical model to use for predicting these effects. This article addresses the variation that one can expect from one model to another and makes recommendations for their appropriate use. Soil nonlinearity and its correct prediction have major importance for seismic hazard assessment and building design requirements. As an example, National Earthquake Hazards Reduction Program (NEHRP) (2000) recommendations for soil site amplication factors decrease to 1.0, reecting expected nonlinear effects, for site classes C and D as effective peak acceleration increases to 0.5g at 5 Hz. However, these amplication factors are being continually reevaluated as new data become available (Crouse, 1995; Borcherdt, 1996; Abrahamson and Silva, 1997; Joyner and Boore, 2000). These studies generally indicate that the current amplication values may imply too much nonlinearity
1609
1610 for stiff-soil sites, leading to an underprediction of ground motion, and that nearly constant factors with increasing effective peak acceleration are more consistent with the current data set. These amplication factors have been determined by modeling data and making forward predictions; therefore, it is important to use the correct theoretical model. Numerical schemes for predicting soil nonlinearity can be divided into equivalent linear and nonlinear approaches. The equivalent linear method is more simplistic and assumes that a damped linear model can approximate dynamic soil behavior through the appropriate choice of material parameters. Nonlinear formulations track a seismic load through stress-strain space using different assumptions for the details of the stress-strain relationship. In addition, nonlinear methods make a distinction between considering pore-water pressure (effective stress calculation) or assuming dry conditions (total stress calculation). Many studies have compared the performance of subsets of these formulations (Constantopoulos et al., 1973; Streeter et al., 1974; Finn et al., 1978; Yoshida, 1994; Yoshida et al., 1995; Yoshida and Iai, 1998; Hashash and Park, 2002; Kausel and Assimaki, 2002). Other studies have considered various aspects of the effects of nonlinear soil response on ground motion (Yu et al., 1993; Ni et al., 1997, 2000; Hartzell, 1998; Su et al., 1998; Silva et al., 1998; Hartzell et al., 2002). In this work we compare the full range of methodologies from equivalent linear, with and without recently developed frequency-dependent moduli and damping ratios, to nonlinear formulations with total and effective stress conditions. The nonlinear calculations are based on newly developed nite difference codes by Bonilla et al. (1998) and Bonilla (2000).
Table 1
NEHRP Site Classications
NEHRP Site Class Vs30 (m/sec)
A B C D E
Figure 1. Average shear-wave velocity proles to 150-m depth for sites with midrange NEHRP site classications (Silva et al., 1998). B/C boundary is also the reference site condition adopted for the 1996 National Seismic Hazard Maps (Frankel et al., 1996). Proles BC, C, D, and E in Figure 1 range from soft rock to soft soil with nonlinearity allowed to extend to a depth of 152.5 m (500 ft). Below this depth, linear wave propagation is assumed. Class B is assumed to be nonlinear to a depth of 22.25 m (73 ft), where the shear-wave velocity reaches 1650 m/sec. For all the proles except B, a smooth gradient extends from a depth of 152.5 m down to a depth of 305 m, bringing the velocities up to 1950 m/sec, the maximum shear-wave velocity for prole B. The depth to which nonlinearity persists inuences the extent of nonlinear behavior observed at the surface, with these effects becoming more pronounced as the depth increases. We do not consider variations in the thickness of the nonlinear column; however, this topic has been investigated elsewhere (Yu et al., 1993; Ni et al., 1997; Silva et al., 1998; Hashash and Park, 2002). In addition to shear-wave velocity information, we need to quantify the degree of nonlinearity within each prole. For this purpose we adopt the widely used modulus reduction, G/Gmax, and damping ratio, n, curves as a function of shear strain developed by the Electric Power Research Institute (EPRI) (1993) for rock and sand and the curves of Vucetic and Dobry (1991) for clay. These curves are shown in Figure 2. The EPRI (1993) curves are an improvement over those previously available in that they are a function of
Method of Comparison
Nonlinear soil models are often compared in terms of stress-strain plots or the variation in maximum stress or strain with depth. We prefer to make our comparisons using a more accessible quantity, namely, the Fourier spectral ratio of ground motion at the surface with respect to the input motion at depth. The NEHRP (1994) soil classication scheme is used to consider the effects of different soil types. Site categories A, B, C, D, and E are based on the average shear-wave velocity in the top 30 m, Vs30, as given in Table 1. Silva et al. (1998) developed average shear-wave velocity proles for these site classes based on available shear-wave velocity measurements. They selected velocity proles with Vs30 values closest to the midrange of the NEHRP site categories. These proles are averaged and shifted so that their Vs30 values agree with the midrange values for the different NEHRP site categories. We have adopted these proles for this study (Fig. 1). Site class A is not considered here because of its limited occurrence in the western United States and because the high surface velocities are expected to result in approximately linear behavior. Class BC is midway between classes B and C and lls the relatively large jump in material properties between these two classes. The NEHRP
1611
Figure 2. EPRI (1993) modulus reduction and damping ratio curves for rock and sand, and the Vucetic and Dobry (1991) curves for clay. The EPRI curves incorporate reduced nonlinearity with depth, and the Vucetic and Dobry curves are for a range in plasticity indexes.
1612 depth, reecting reduced nonlinearity with depth of burial. The EPRI (1993) rock curves are used for site classes B and BC, the sand curves for site classes C and D, and the Vucetic and Dobry (1991) curves for a plasticity index of 30 for site class E. By this association we consider C and D sites to be cohesionless soils and E sites to be cohesive soils, such as, bay mud. In addition to these classications we also consider the Peninsular Range modulus reduction and damping ratio curves developed by Silva et al. (1997) for site classes C and D, referred to here as Cpen and Dpen. These curves were empirically developed to reect lower levels of nonlinearity for the 1994 Northridge earthquake than predicted by the EPRI (1993) cohesionless soils curves. They are formed by using the 15.5 to 36.5 m (51 to 120 ft) depth range EPRI curve for depths 0 to 15.2 m (0 to 50 ft), and the 152.7 to 304.8 m (5011000 ft) depth range EPRI curve for depths below 15.2 m (50 ft). Modulus reduction and damping ratio curves sufce to quantify the nonlinear response in an equivalent linear formulation and in a nonlinear formulation that is tied to a prescribed modulus reduction curve. However, other more general nonlinear formulations do not utilize modulus reduction curves. We consider both types in this article. For the input motion to our soil proles we use the northsouth component of the rock site record from station LA00 (Stone Canyon Reservoir, 34.106 N, 118.454 W) for the 1994 Northridge earthquake (Fig. 3). This station is located in the Santa Monica Mountains on Mesozoic rock and has a peak acceleration of 0.26g. We consider this record to be typical of rock site ground motion for a magnitude 6 to 7 earthquake, without any unusual behavior such as isolated large spikes or pulses. To obtain peak input motions from 0.1 to 0.9g, this record is scaled up or down. Other studies of nonlinear soil behavior have used an impulse base excitation (Yu et al., 1993; Ni et al., 1997). Because we are investigating nonlinear processes where the solution is computed in the time domain, we consider an actual earthquake record to be a more appropriate input. For the equivalent linear formulations, because they are linear, we can calculate the transfer function of the soil column directly, and thus the amplication values. For the nonlinear methods we take the Fourier spectral ratio of the surface motion with respect to the input motion at depth. The spectral ratios are then smoothed using a 0.25-Hz-wide running mean lter to remove spikes caused by spectral holes in the denominator. In all cases we remove the factor-of-two free-surface effect (Aki and Richards, 1980).
Figure 3.
Acceleration record for the 1994 Northridge earthquake for the Stone Canyon Reservoir rock site, northsouth component. This record is scaled and used as input to the average velocity proles in Figure 1.
Soil Models
Equivalent Linear: SHAKE91 The rst site response formulation considered is the code SHAKE91 (Idriss and Sun, 1992), based on the original code by Schnabel et al. (1972). We modied the SHAKE91 algorithm making parts double precision to improve stability
at high-input ground motions. Because of its ease of application, the equivalent linear approach, and particularly the SHAKE91 code, has been widely used. With the equivalent linear approach the material properties of shear modulus, G, and damping, n, are constant within each layer and for the entire duration of the earthquake, regardless of whether the strains are large or small at a given time. An iterative process is used to determine the values of G and n for each layer that are consistent with the level of strain induced in the layer. Because experimental soils data are based on simple harmonic loading with constant peak shear strain, an effective shear strain is used for earthquake records that is a fraction, , of the maximum shear strain. A value of 0.65 is typically used for and is assumed in this study. The effective shear strain is then used to calculate G and n from empirical modulus reduction and damping ratio curves. Yoshida and Iai (1998) have pointed out two weaknesses of the equivalent linear method. First, damping is constant, independent of frequency. For an unbounded, homogeneous medium, the shear strain is given by the particle velocity divided by the shear wave velocity of the medium. The spectrum of particle velocity, and therefore of shear strain, decreases in amplitude approximately as 1/f at high frequency. Then, the damping obtained from standard laboratory curves should also decrease at high frequencies. Viewed alternatively, the energy dissipated per loading cycle is proportional to the area enclosed by the hysteresis loop. Therefore, low-amplitude, high-frequency components of strain should have much less damping. The overattenuation of high frequencies by the equivalent linear model becomes more pronounced at high levels of strain. Recently, Sugito et al. (1994), Joyner and Boore (1998), and Joyner (1999) have formulated equivalent linear methods that use frequency-dependent damping. Second, equivalent linear analysis can overestimate resonant-frequency amplitudes and shear stresses compared with nonlinear analysis. Because soil stiffness varies as a function of time in nonlinear analysis, resonances are less likely to develop, compared with
1613 ergy at a strain level of ca is ca f (ca)/2 and the energy loss per cycle is obtained by integrating f (c) over a complete cycle yielding (Ishihara, 1996) DW 8(
ca 0
the constant-layer-parameter approach of equivalent linear analysis. Equivalent Linear: TREMORKA Recognizing the overdamping that occurs with the traditional equivalent linear approach, Kausel and Assimaki (2002) have recently advanced a formulation with frequencydependent moduli and damping. Except for the manner in which the frequency-dependent moduli and damping are chosen, their approach is the same as Schnabel et al. (1972). A two-parameter, smooth functional form is used to t the strain spectrum within each layer. This function is then used to read off strain values as a function of frequency, which in turn are used to specify the frequency-dependent moduli and damping from standard laboratory curves. The process is iterated in the same way as SHAKE91 to obtain internally consistent soil parameters. We have implemented this method in the code TREMORKA. There is value in obtaining an improved equivalent linear code. The soil parameters needed by some nonlinear codes can be more numerous and difcult to obtain than the modulus reduction and damping curves used by equivalent linear analysis. The frequency-dependent equivalent linear formulation requires only two parameters, the shear modulus reduction and fraction of damping as a function of shear strain. Both of these quantities are routinely measured in laboratory studies of soils. Also, in some cases the equivalent linear method can be more numerically stable than a nonlinear calculation. Nonlinear: NOAHW The NOAHW code is the rst of three newly developed, second-order, staggered-grid nite difference codes considered in this study. This code, and other nonlinear methods, operates in the time domain by tracking the earthquake load through stress-strain space. In NOAHW, the Iwan (1967) model is implemented to take into account soil nonlinearity. Furthermore, in this model the stress-strain relation is specied by a given modulus reduction curve. The main difference between this formulation and the one implemented by Joyner and Chen (1975) is that real G/Gmax laboratory data can be assigned to each layer. In addition, the implemented Iwan (1967) model is one that uses a series of parallel elastic springs that permits the computation of stress from the strain. For dynamic nonlinear wave propagation with a nite difference formulation, it is easier to obtain the strain from the node displacement gradient and then introduce it into the constitutive equation to compute the stress. This type of formulation is advantageous when the only data available are the G/Gmax curves. In our comparisons, the same modulus reduction curves are used as in the preceding equivalent linear runs. Nonlinear damping comes from the hysteretic loops, however, not from a damping ratio curve. The damping ratio is given by DW/(4pW), where W is the maximum stored energy and DW is the energy dissipated per cycle (Ishihara, 1996; Kramer, 1996). The maximum stored en-
f(c)dc
W),
where f (c) is the functional relationship between the stress and the strain. For the Iwan (1967) model the hysteresis loops follow the so-called Masing criteria (Masing, 1926). If f (c) is the hyperbolic model, as we assume, the damping approaches 2/p 0.63 at large strain (Ishihara, 1996). Nonlinear: NOAHH In the previous section we mentioned that the hysteresis loops follow the Masing (1926) criteria. In the original Masing (1926) rules the initial loading or backbone curve is dened by a functional relationship between the stress, s, and the strain, c, as s f (c). Subsequent unloading and reloading cycles may be expressed as, (1/jH)(s sa) f [(1/jH)(c ca)],
where (ca, sa) is the unload/reload reversal point and jH is a hysteresis scale factor. In Masings original formulation, the hysteresis scale factor is equal to 2. Subsequently, some additional constraints have been added to prevent the computed stress from exceeding the maximum strength of the material, smax. The ensemble of these criteria is known as the extended Masing rules (Kramer, 1996). Several articles have dealt with the implementation of the extended Masing rules in a phenomenological way (Pyke, 1979; Li and Liao, 1993; Finn, 1982; Vucetic, 1990; Iai et al., 1990a). In comparison, the Iwan (1967) model fully complies with these criteria by a mechanical model of a series of parallel or serial elastic springs. Laboratory tests of soils under cyclic shear deformation show that the maximum damping ratio varies between 0.2 and 0.4 depending on the type of material. The Iwan (1967) model with a hyperbolic backbone, on the other hand, may overdamp when large strains are reached. Thus, we would like to nd a compromise between the shear modulus reduction and the damping ratio to better characterize the soil nonlinearity. Bonillas (2000) formulation utilizes a variable hysteresis scale factor that assures the stress-strain path during each loading/reloading is bounded by the maximum shear strength. In addition, subsequent unloading or reloading curves follow the backbone curve if the previous maximum shear stress is exceeded. This phenomenological representation is called the generalized Masing rules because its formulation contains Pykes (1979) and the original Masing jH 2 models as special cases. Furthermore, this formulation allows, by controlling the hysteresis scale factor, the reshaping of the backbone curve as suggested by Ishihara et al. (1985) so that the hysteresis path follows a prescribed damping ratio. Given the laboratorymeasured soil damping, and since we know the damping as a function of the reference strain and hysteresis shape factor,
1614 we equate them to nd a new reference strain compatible with the damping value. This reference strain is used to compute a new backbone curve. Thus, the hysteresis loop will comply with the observed damping ratio. This procedure is followed at each time step. More details of this model are given by Lavallee et al. (2003). The functional form of the backbone curve for NOAHH is the hyperbolic model given by (Kondner and Zelasko, 1963), s Gmaxc Gmax 1 c smax Gmaxc , c 1 cref
Table 2
Water Table Depths
NEHRP Site Class NOAHH Depth (m) NOAHB Depth (m)
B BC C D E
Not considered Not considered Not considered 15.0 and 2.0 5.0 and 2.5
Ko
(1
We assume a value for of 30 and normal consolidation (OCR 1), so that Ko is then 0.5. Gmax is calculated from the density, q, and the shear-wave velocity, Vs, by the expression qVs2. rv is the vertical effective stress, rv rt qwg(z zw),
where Gmax is the low-strain shear modulus, smax is the maximum shear stress that the material can support in the initial state, and cref smax /Gmax is the reference strain. The hyperbolic model has been extensively used because it is relatively easy to implement and because of the low number of parameters necessary for its description (Hardin and Drnevich, 1972; Pyke, 1979; Ishihara, 1996). This code, and the following one, uses the multispring model plane strain formulation for the stress-strain relationship in each soil layer (Towhata and Ishihara, 1985; Iai et al., 1990a). Each spring follows the hyperbolic model and implements hysteresis by the generalized Masing rules described previously. The multispring model was introduced by Towhata and Ishihara (1985) to model cyclic mobility and dilatancy in sands. However, the NOAHH code is for total stress analysis; no pore-water pressure is considered. Under these conditions, the multispring model gives the same results as the single hyperbolic element model customarily applied. NOAHH uses the multispring model for total stress analysis because it is relatively easy to implement, it has a two-dimensional representation (plane strain), and it takes into account the anisotropy of soils (vertical and horizontal initial stresses may be different). In NOAHH we still specify typical water-table depths from the San Fernando Valley, California (Gibbs et al., 1999) (Table 2), because the water table inuences the calculation of the over burden. The maximum shear stress for cohesionless soils is calculated by the expression, smax
where rt is the total vertical stress, qw is the density of water, g is the acceleration of gravity, z is the sample depth, and zw is the depth of the water table measured from the ground surface. Although the nonlinear formulation in NOAHH does not follow a particular modulus reduction curve, as NOAHW does, we produce a soil column with the same strength implied by the modulus reduction curves used previously through the addition of a cohesion term. With the addition of cohesion, we are able to model stiff soils, such as encountered with site class B, and still have reasonable values of Ko and . The MohrCoulomb relation is given by smax c rn tan ( ), where c is the cohesion and rn is the normal stress, or alternatively (Jaeger and Cook, 1979), smax rm sin c(cos ), where rm (rv /3)(1 2Ko) is the effective mean stress. This expression for smax is the same as the one given by Hardin and Drnevich (1972) with an axis transformation of the Mohrs circle. Thus we have an expression for the cohesion in terms of the reference strain, the coefcient of the earth at rest, and the angle of internal friction, c (cref Gmax rm sin u)/cos u.
((1
where Ko is the coefcient of lateral stress at rest, and is the angle of internal friction. Ko may be represented in terms of the angle of internal friction and the overconsolidation ratio as (Kulhawy and Mayne, 1990),
Then, the procedure followed is to rst determine cref from the desired modulus reduction curve as the value of strain when G/Gmax 0.5. The cohesion is determined using the preceding expression. Finally, smax is calculated using the expression of Jaeger and Cook (1979). For site classes B, BC, C, and D (rock and sand) a maximum damping ratio of 0.30 is used. For site class E (clay) the maximum damping ratio is 0.20 (Fig. 2). Realistic intrinsic attenuation with constant Q is incorporated by the rheology of the generalized Maxwell body (Day, 1998). In addition, a small amount of viscous damping may be added at low-strain levels for numerical stability (Li, 1990).
1615
Nonlinear: NOAHB The NOAHB code is the third and nal level of complexity in the nonlinear formulations we consider. NOAHB treats undrained conditions of effective stress and also utilizes the Bonilla (2000) generalized Masing rules with the hyperbolic model discussed previously. The algorithm is based on the multispring mechanism model rst implemented by Towhata and Ishihara (1985) to simulate pore pressure generation in sands under cyclic loading and undrained conditions, and further developed by Iai et al. (1990a,b) to take into account the cyclic mobility and dilatancy of sands. Under cyclic loading, as the pore pressure increases, the effective stress decreases, and as a result, Gmax and smax decrease. As Gmax and smax change, so does the shape of the backbone curve. An empirical approach developed by Towhata and Ishihara (1985) and Iai et al. (1990a) is utilized to model the decrease of effective stress due to the increase in pore pressure. Six parameters ( p, p1, p2, w1, S1, and c1), together with the porosity are used to describe the dilatancy of the material. Table 3 lists the values we have adopted in this study. They are assumed to be the same for all layers. Iai et al. (1990b) give a more detailed description of these model parameters. We consider two values for the angle of internal friction, 30 and 40 (corresponding to a medium-dense sand and a dense sand, respectively), to cover the range of typical values, and Ko 0.5, and a maximum damping ratio of 0.3. With the introduction of pore pressure and effective stress it is no longer possible to make a detailed comparison with the previous equivalent linear formulations or the total stress nonlinear formulations. Signicant buildups in porewater pressure make a distinctly different problem. However, for completeness we include this case. Because the effects of pore pressure generation and liquefaction are signicant only in cohesionless soils with low to medium strength, we limit our tests to site classes D and E. To this point class E has been considered a clay soil with the appropriate modulus reduction curves. For the NOAHB calculations we drop this association because a clay soil is not expected to undergo liquefaction. For NOAHB, site class E is a cohesionless sandy soil similar to class D, only with lower shear-wave velocities. We consider two water-table depths (Table 2) to investigate the effect of the depth to the saturated column. Measurement of resistance to liquefaction of a sand is a standard laboratory test. In this procedure, an undrained sample is subjected to cyclic loading. The pore-water pressure builds up steadily as the cyclic stress is applied, eventually approaching the initial conning pressure, leading to an axial strain of about 5% peak-to-peak (or double amplitude, DA). This state is referred to as the start of liquefaction. Figure 4 shows the analytical liquefaction resistance curves for cyclic simple shear and soil classes D and E and our assumed dilatancy parameters listed in Table 3. These curves plot the number of cycles to produce 5% DA shear deforp
Table 3
Dilatancy Parameters for NOAHB
Porosity p1 p2 w1 S1 c1 28.0 0.45 0.5 0.6 5.0 0.01 1.0
mation versus the cyclic stress ratio, sxy /rm0, where sxy is the cyclic shear stress and rm0 is the initial effective conning stress. Figure 4 shows the numerical response of the multispring model in NOAHB for our hypothetical proles. All simulations are done with a reference effective mean stress, rma, equal to the initial effective conning stress, rm0, which is computed at the middle of each layer for soil models D and E. The shaded area for a given angle of internal friction represents the ranges of the liquefaction resistance curves for all the layers within the prole. We can observe that the stronger the soil is, that is 40, the larger the cyclic stress ratio is needed to deform the material. In addition, the effect of the water table is more evident when 40. Indeed, a deeper water table produces less buoyancy, and thus the material becomes stronger. Conversely, when the soil is relatively weak, 30, and combined with a shallow water table, we can observe that there is little difference between the liquefaction resistance curves for the shallow and deeper layers. We recognize that the results strongly depend on the chosen dilatancy parameters. However, in the absence of specic geotechnical data on the liquefaction resistance curves, we show the effect of pore-water pressure for these two soil classes.
Results
Figures 5 through 9 plot amplication as a function of frequency for the site classes discussed previously (B, BC, C, Cpen, D, Dpen, and E) and peak base input motions of 0.1g to 0.9g in 0.1g increments. The linear response curves, designated by an L, are also shown for each site classication. Consider rst the SHAKE91 results in Figure 5. As expected, the low input acceleration curves most closely follow the linear curve, but they diverge more as the stiffness of the soil decreases, from class B (rock) to E (clay). Most of the equivalent linear curves lie below the linear response because of higher damping at higher strain levels. However, resonant peaks at lower frequencies (below about 3 Hz) can rise above the linear response caused by the familiar downshift in energy from the reduction of the shear modulus. Also, for site class B, higher input motions yield higher amplications than the linear case. This result comes from the competing effects of reduction in the shear modulus and increased damping. For class B the reduction of the shear modulus leads to an increase in amplication that is greater than
1616
Figure 4. Liquefaction resistance curves for soil classes D and E. The strength of the material is related to the angle of internal friction, , shown in each plot. For all models the reference effective mean stress, rma, is equal to the initial effective conning stress, rm0. The dilatancy parameters are listed in Table 3. The shaded area represents the range of liquefaction resistance curves for all the layers within the prole.
the attenuation caused by the increased damping. However, the most striking feature of the SHAKE91 curves is the severe attenuation for softer soils and larger input motions. For site classes D and E, amplication values are in the range from 10 3 to 10 4 for frequencies greater than 6 Hz and input motions above 0.5g. This response has the effect of removing nearly all high frequencies in areas with deep soil proles and runs counter to observations (Durward et al., 1996). As discussed previously, and shown in Figure 6, the TREMORKA amplication curves have signicantly less attenuation for soft soils and large strains than the SHAKE91 results. The moduli and damping ratios are now a function of frequency based on the spectrum of the strain induced in each layer by the earthquake motion. For site class B the difference between SHAKE91 and TREMORKA is negligible, indicating that constant moduli and damping ratios are a good approximation for this site class. For all the other site classes, however, there are signicant differences. Site classes C and Cpen have nearly constant amplication at frequencies above 4 Hz for different input motions. The amplication for class D never goes much below 1.0 and is about 0.1 for class E at high frequencies and large input motions. Figure 7 shows amplication spectra for NOAHW, a nonlinear formulation that utilizes the same modulus reduction curves as the previous equivalent linear models to describe the backbone stress-strain curve. These amplication curves are more complicated because they are calculated from the spectral ratios of the output to the input motions and because of the more complicated response produced by tracking the load through stress-strain space. Smoothing out the peaks in the curves and considering only the general amplication levels, the NOAHW results for class B are very similar to the TREMORKA results. For the softer site classes, C and softer, the NOAHW results give amplication values up to a factor of 2 less than those from TREMORKA at frequencies above about 2 Hz. This difference is due to the different way in which damping is incorporated into the two codes. TREMORKA uses damping values picked from an empirical curve. NOAHW uses the dissipation of energy in each hysteresis loop following the Masing criteria and the hyperbolic model, thus at large strains the maximum damping is 2/p. The NOAHH amplication curves are shown in Figure 8. This code uses a hyperbolic backbone curve and adds cohesion to obtain the same strength as implied by the empirical modulus reduction curves. For the stiffer proles (B, BC, C, and Cpen) there is very little difference between the NOAHW and NOAHH results for both small- and large-input motions. For site classes D and E there is a divergence between the two solutions, with NOAHH yielding lower amplication values above 2 Hz by 50% at low input motions and up to a factor of 2 at large input motions. Recall that NOAHH constrains the maximum damping ratio to 20% for clays and 30% for sands. This difference is discussed further next when we consider time domain records.
1617
Figure 5. Equivalent linear SHAKE91 spectral amplication curves as a function of frequency for the average NEHRP velocity proles in Figure 1 and for peak input motions from 0.1g to 0.9g. An L indicates the linear response curve. Cpen and Dpen use the site class C and D velocity proles, respectively, but modulus reduction and damping curves with reduced nonlinearity (see text and Silva et al., 1997). The results for class E are plotted on both linear and log scales to reveal details of the curves.
1618
Figure 6. Same as Figure 5 for the equivalent linear code TREMORKA with frequency-dependent moduli and damping ratios.
1619
Figure 7. Same as Figure 5 for the nonlinear total stress code NOAHW based on the Iwan (1967) model and hysteresis loops according to the Masing (1926) criteria with prescribed modulus reduction curves.
1620
Figure 8. Same as Figure 5 for the nonlinear total stress code NOAHH with the hyperbolic model backbone curve and hysteresis loops following generalized Masing rules and added cohesion to match given modulus reduction curves.
1621
Figure 9. Same as Figure 5 for the nonlinear effective stress code NOAHB with two assumptions for the depth of the water table and two different values for the angle of internal friction, phi.
of pore-water pressure. We consider only soft soils with class D and E velocity proles. A deeper and a shallower water table are considered for both site classes (Table 2); 15 m and 2 m for class D, and 5 m and 2.5 m for class E. We also consider two angles of internal friction, 30 and 40, which covers the range of typical values. As the water table shallows the trend is for greater attenuation at higher frequencies above 1 to 2 Hz, and somewhat larger resonant peaks below 1 Hz. Decreasing the angle of internal friction has the effect of increasing the attenuation at nearly all frequencies. This effect is expected because decreasing the angle of internal friction means a reduction of the initial strength as well. Figure 10 summarizes the preceding comparisons by plotting a low-strain curve (0.2g) and a high-strain curve (0.7g) from each of the ve different nonlinear soil models. For site class B differences between the different formulations are not signicant. For this comparison we consider mean levels and average over higher-frequency oscillations. All the curves are also similar for site class C for an input motion of 0.2g. However, there are signicant differences between the equivalent linear solutions and the nonlinear solutions for an input motion of 0.7g. The two equivalent linear solutions overestimate the amplitude of the resonant peak at 1.0 Hz. This tendency has been noted by other authors (Yoshida and Iai, 1998). For the softer soil proles (classes D and E), the SHAKE91 solution becomes more divergent from the others, particularly at larger strains. The other solutions group more closely together, with the TREMORKA curves generally lying above the nonlinear solutions. The NOAHB curves for an angle of internal friction of 30 are also included for site classes D and E. An angle of internal friction of 30 was also used for the total stress code NOAHH. NOAHW does not use this parameter. However, as we mentioned previously, when signicant pore-water pressures develop, results from the effective stress calculation of NOAHB will depart from the total stress analysis. We include it here to see the effect of pore-water pressure. The NOAHB curves in Figure 10 are for the deeper water-table values (15 m for class D and 5 m for class E). For these parameters the NOAHB results are generally intermediate to the two total stress calculations. Figure 11a,b shows the same comparisons in the time domain as in Figure 10. As noted previously, there is very little difference between the different solutions for site class B at both low and high input motions. The same can be said for class C at lower input motions (0.2g). However, starting with class C at high input motions (0.7g), differences develop and grow with the reduction in the strength of the soil column. For classes D and E, the SHAKE91 solutions clearly lack high frequencies that are present in the TREMORKA and nonlinear solutions. The TREMORKA and SHAKE91 solutions yield similar peak accelerations for all site classes and input motions, even for classes D and E. This observation is
explained by the fact that the peak acceleration is controlled by the frequency band 1 to 3 Hz, where the two solutions have similar amplitudes. The nonlinear solutions, on the other hand, all have signicantly lower amplitudes in this frequency band, and thus lower peak accelerations. The peak accelerations for the nonlinear formulations are all similar. The increased attenuation of higher frequencies, above about 2 Hz, for the NOAHH solution compared with the NOAHW solution is apparent for site classes D and E. The input record (Fig. 3) shows a progression from lower-amplitude, higher-frequency motion at the beginning to larger-amplitude, lower-frequency motion with the arrival of the main shear-wave energy. This character is typical of earthquake ground motion. The nonlinear soil responses tend to amplify the pattern of higher frequencies early and lower frequencies late, because the larger, late motions cause greater nonlinearity with resulting lower rigidity and higher damping. This pattern of ground motion has also been observed in earthquake records on soft soils (Satoh et al., 1995; Frankel et al., 2002). However, large input motions with soft soils can result in cyclic mobility with late, high-frequency, cusped waveforms (Archuleta, 1998). Both characteristics are seen in Figure 11a,b. The fact that there is damping control in NOAHH does not necessarily mean that the result will be less damped than the Iwan (1967) model. Indeed, what happens, especially at shallow depths or where there is a low-strength layer, is that large strains may develop. This can be understood because, to keep the same energy (hysteresis loop area) in a damped controlled system, the material has to develop larger strains than in one following the traditional Masing criteria. These deformations, typically obtained at shallow depths, result in reduced rigidity that consequently damps the higher frequencies in the ground motion. This effect has been pointed out by Ishihara et al. (1985).
Nisqually
The 28 February 2001 M 6.8 Nisqually earthquake produced a ground motion data set that shows signicant soil nonlinearity (Frankel et al., 2002), manifested as shifts in resonant spectral peaks at lower frequencies, increased attenuation at higher frequencies, and cusped waveforms after the arrival of the S wave (Archuleta, 1998). An ML 3.4 aftershock 14 hr after the mainshock and recorded at the same stations provides a linear response reference point. We use these records for comparison with our model predictions. Figure 12 shows the stations we consider. In terms of the surcial geology, the stations fall into the two general categories of young, soft sediments, consisting of articial ll and Holocene alluvium (BOE, HAR, KDK, and SDS), and older, stiffer sediments of Pleistocene age (BHD, LAP, and THO). The soft sediment sites had peak accelerations from 0.19g to 0.22g and the Pleistocene sediments had peak accelerations from 0.11g to 0.16g. Following Frankel et al. (2002), we estimate the amplication at each site by calcu-
1623
Figure 10. Comparison of amplication spectra for peak input motions of 0.2g and 0.7g for the different equivalent linear and nonlinear soil response formulations considered in this article. Black, SHAKE91; red, TREMORKA; green, NOAHW; blue, NOAHH; magenta, NOAHB. The NOAHB curves assume the deep-water table and an angle of internal friction of 30.
1624
Figure 11. Comparison of time domain surface records for the spectral amplication curves in Figure 10. (a) peak input motion of 0.2g. (b) peak input motion of 0.7g.
1625
Figure 12. Seattle area map showing surcial geology and locations of accelerometer stations that recorded the 2001 Nisqually earthquake and its large aftershock, and for which there is also a measure of the shallow shear-wave velocity.
lating the Fourier spectral ratios with respect to site SEW on Tertiary sedimentary rock (Fig. 12). The spectra are based on the rms of the two horizontal components of acceleration and 40 sec of the S-wave record. This duration ensures that we have the average spectral levels over all signicantamplitude shear-wave energy. The ratios are smoothed with the same 0.25-Hz running mean lter as our nonlinear theoretical amplication spectra. Because both the mainshock and aftershock were deep (52 km and 54 km, respectively) and within a few kilometers of each other, differences in the source-to-station pathlengths are small, and no propagation path corrections are considered necessary. Figure 13 shows mainshock and aftershock spectral ratios for the soft-soil and stiff-soil sites. For the stiff-soil sites on Pleistocene sediments there is little difference between the mainshock and aftershock spectral ratios. For the soft-soil sites, however, frequency downshifts of resonant spectral peaks at lower frequencies are apparent, as well as greater attenuation at higher frequencies. Williams et al. (1999) have determined Vs30 values (average shear-wave velocity in the top 30 m) for the seismic stations (Table 4). These values allow us to specify the NEHRP site class for each station. In Figure 14
we use this information to compare the observed mainshock spectral ratios with the nonlinear predictions. Our objective in this comparison is to see how well the average NEHRP site class velocity proles perform without detailed velocity information at the site. For site classes D and E, we use the effective stress results from NOAHB and the shallow-water table prole with 40. These sites all lie within the Duwamish River Valley, which is known from P-wave seismic refraction data to have a water table in the 2- to 4-m depth range (Williams et al., 1999). The site class C stations are on Pleistocene sediments with a deeper water table. For these stations we use the total stress NOAHH results. Recall, however, that all the formulations considered give similar results for class C with lower input motions. Assuming that the Tertiary sedimentary rock site SEW is site class B, we calculate the theoretical amplication spectra by dividing by the site class B response. Two input ground motion levels are considered in Figure 14, 0.1g and 0.2g, which bracket the range of values in the Seattle area for the Nisqually earthquake. At lower frequencies, below 2 Hz, there are differences between the frequencies and amplitudes of resonant spectral peaks that are due to the details of the velocity struc-
1626
Figure 13. Nisqually mainshock (solid curves) and aftershock (dashed curves) shear-wave spectral ratios relative to the Tertiary rock site SEW. Spectra are based on 40 sec of the rms of the two horizontal components of acceleration.
Table 4
Vs30 for Ground Motion Recording Stations in Seattle, Washington
Station Site Class Vs30 (m/sec) Distance of Vs30 from Station (m)
Figure 14. Comparison of observed mainshock spectral ratios from Figure 13 (dashed and dotted curves) with model predictions for site classes C, D, and E for two peak input motions of 0.1g and 0.2g (solid curves). The model predictions for class C are from the nonlinear total stress code NOAHH. The model predictions for classes D and E are from the nonlinear effective stress code NOAHB with the shallow water table and angle of internal friction of 40.
E E E D C C C
ture at each station, not contained in our smooth, average velocity proles. Using high-resolution S-wave seismic reection data near stations BOE, SDS, and KDK, Williams et al. (2000) and Frankel et al. (2002) found that the resonant spectral peaks below 2 Hz are consistently observed during weak and strong motion and are likely caused by ll/allu-
1627 class C velocity prole to depths of 100 m or more in the Seattle Pleistocene glacial deposits.
vium or young alluvium/bedrock impedance boundaries located in the 16- to 50-m depth range. Because these impedance boundaries are not captured in the average velocity proles used in this study, the nonlinear models underpredict the observed peak amplications by up to a factor of 2 at these resonant peaks for site classes E and C. However, the amplication levels over broader frequency bands are generally well represented, especially considering that no iterations were performed to ne tune the theoretical results to the Nisqually data. There is a trend for lower observed amplication at the class C sites. Two recent boreholes in the Pleistocene sediments on the hills east of downtown Seattle show nearly uniform shear-wave velocities of 500 m/sec to a depth of at least 100 m (Odum et al., 2003). These velocities are signicantly lower than our average site class C prole and would explain the greater attenuation.
Acknowledgments
We thank Walt Silva for supplying his average velocity proles for the NEHRP site categories. S.H. beneted from discussions with Pengcheng Liu and Alena Leeds. This article was improved by reviews from Chris Cramer, Charles Mueller, and two anonymous reviewers.
References
Abrahamson, N. A., and W. J. Silva (1997). Empirical response spectral attenuation relations for shallow crustal earthquakes, Seism. Res. Lett. 68, 94127. Aguirre, J., and K. Irikura (1997). Nonlinearity, liquefaction, and velocity variations of soft soil layers in Port Island, Kobe, during the Hyogoken-Nanbu earthquake, Bull. Seism. Soc. Am. 87, 12441258. Aki, K., and P. G. Richards (1980). Quantitative Seismology Theory and Method, W. H. Freeman, San Francisco. Archuleta, R. J. (1998). Direct observation of nonlinearity in accelerograms, in The Effects of Surface Geology on Seismic Motion, Vol. 2, K. Irikura, K. Kudo, H. Okada, and T. Sasatani (Editors), Balkema, Rotterdam, 787792. Beresnev, I. A., E. H. Field, P. A. Johnson, and K. E.-A. Van Den Abeele (1998). Magnitude of nonlinear sediment response in Los Angeles basin during the 1994 Northridge, California, earthquake, Bull. Seism. Soc. Am. 88, 10791084. Bonilla, L. F. (2000). Computation of linear and nonlinear site response for near eld ground motion, Ph.D. Thesis, University of California, Santa Barbara. Bonilla, L. F., D. Lavallee, and R. J. Archuleta (1998). Nonlinear site response: Laboratory modeling as a constraint for modeling accelerograms, Proc. of Second International Symposium on the Effects of Surface Geology on Seismic Motion, Editors: K. Irikura, K. Kudo, H. Okada, and T. Sasatani, vol. 2, 793800. Borcherdt, R. D. (1996). Preliminary amplication estimates inferred from strong ground-motion recordings of the Northridge earthquake of January 17, 1994, in Proc. of International Workshop on Site Response Subjected to Strong Earthquake Motions, Yokosuka, Japan, 1617 January 1996. Constantopoulos, I. V., J. M. Roesset, and J. T. Christian (1973). A comparison of linear and exact nonlinear analysis of soil amplication, in Proc. of 5th World Conference on Earthquake Engineering, Rome, Italy, 18061815. Crouse, C. B. (1995). Site response studies for purpose of revising NEHRP seismic provisions, Data Utilization Report CSMIP/9503, California Strong Motion Instrumentation Program, 68 pp. Day, S. M. (1998). Efcient simulation of constant Q using coarse-grained memory variables, Bull. Seism. Soc. Am. 88, 10511062. Durward, J. A., D. M. Boore, and W. B. Joyner (1996). The amplitude dependence of high-frequency spectral decay: constraint on soil nonlinearity, in Proc. of International Workshop on Site Response Subjected to Strong Earthquake Motions, Yokosuka, Japan, 1617 January 1996, Vol. 2, 82103. Electric Power Research Institute (EPRI) (1993). Guidelines for determining design basis ground motions, Electric Power Research Institute Technical Report EPRI TR-102293. Field, E. H., P. A. Johnson, I. A. Beresnev, and Y. Zeng (1997). Nonlinear ground-motion amplication by sediments during the 1994 Northridge earthquake, Nature 390, 599602. Field, E. H., Y. Zeng, P. A. Johnson, and I. A. Beresnev (1998). Pervasive nonlinear sediment response during the 1994 Northridge earthquake: observations and nite-source simulations, J. Geophys. Res. 103, 26,86926,883.
Conclusions
Nonlinear soil effects have been examined as measured by Fourier amplication spectra for equivalent linear formulations with and without frequency-dependent moduli and damping ratios and nonlinear formulations with and without pore-water pressure. No signicant differences exist between these solutions for NEHRP site class B and peak input motions from 0.1g to 0.9g. Amplication spectra within the frequency band 1 to 2 Hz and average spectral levels above 2 Hz are essentially the same for site classes BC and C at lower peak input motions of 0.1g to 0.2g. As the input motion increases, differences develop between the solutions. For the softer soils (C, D, and E), the equivalent linear solutions can have up to 50% larger amplitude resonant peaks below 2 Hz than nonlinear formulations because they include no time dependence in the solution. For site classes D and E, the overdamping of the SHAKE91 solution above about 4 Hz is apparent for the entire range of input motions considered. The equivalent linear formulation with frequency-dependent moduli and damping ratios under damps relative to the nonlinear models considered for site class C with larger input motions and most input levels for site classes D and E. At larger input motions the underdamping for site classes D and E is not as severe as the overdamping with the frequency-independent formulation, but there are still signicant differences in the time domain. Based on these comparisons, a nonlinear approach is recommended for site classes D and E, and for site class C for input motions greater than a few tenths of the acceleration of gravity. Our average site class D and E amplication spectra based on a nonlinear effective stress calculation and shallow water table are consistent with observed soil-to-rock spectra for the Nisqually earthquake, excluding low-frequency resonances. Class C sites for the Nisqually earthquake show greater attenuation at frequencies above 2 Hz than the predictions of a total stress nonlinear calculation. This difference is attributed to velocities lower than the average site
1628
Finn, W. D. L. (1982). Dynamic response of saturated sands, in Soil Mechanics: Transient and Cyclic Loads, G. N. Pande and O. C. Zienkiewicz (Editors), J. Wiley and Sons, New York, 105131. Finn, W. D. L., G. R. Martin, and M. K. W. Lee (1978). Comparison of dynamic analysis of saturated sand, in Proc. of ASCE Geotechnical Engineering Division of Specialty Conference on Earthquake Engineering and Soil Dynamics, Pasadena, California, 1921 June, 472 491. Frankel, A. D., D. L. Carver, and R. A. Williams (2002). Nonlinear and linear site response and basin effects in Seattle for the M 6.8 Nisqually, Washington, earthquake, Bull. Seism. Soc. Am. 92, 2090 2109. Frankel, A., C. Muller, T. Barnhard, D. Perkins, E. V. Leyendecker, N. Dickman, S. Hanson, and M. Hopper (1996). National seismic hazard maps: documentation, June 1996, U.S. Geol. Surv. Open-File Rept. 96-532, 110 pp. Fukushima, Y., K. Irikura, T. Uetake, and H. Matsumoto (2000). Characteristics of observed peak amplitude for strong ground motion from the 1995 Hyogoken-Nanbu (Kobe) earthquake, Bull. Seism. Soc. Am. 90, 545565. Gibbs, J. F., J. C. Tinsley, D. M. Boore, and W. B. Joyner (1999). Seismic velocities and geological conditions at twelve sites subjected to strong ground motion in the 1994 Northridge, California, earthquake: a revision of OFR 96-740, U.S. Geol. Surv. Open-File Rept. 99-446, 142 pp. Hardin, B. O., and V. P. Drnevich (1972). Shear modulus and damping in soils: design equations and curves, Proc. Am. Soc. Civil Eng. 98, 667 692. Hartzell, S. (1998). Variability in nonlinear sediment response during the 1994 Northridge, California, earthquake, Bull. Seism. Soc. Am. 88, 14261437. Hartzell, S. A. Leeds, A. Frankel, R. A. Williams, J. Odum, W. Stephenson, and W. Silva (2002). Simulation of broadband ground motion including nonlinear soil effects for a magnitude 6.5 earthquake on the Seattle fault, Seattle, Washington, Bull. Seism. Soc. Am. 92, 831853. Hashash, Y. M. A., and D. Park (2002). Viscous damping formulation and high frequency motion propagation in non-linear site response analysis, Soil Dyn. Earthquake Eng. 22, 611624. Iai, S., Y. Matsunaga, and T. Kameoka (1990a). Strain space plasticity model for cyclic mobility, Rep. Port Harbor Res. Institute 29, 2756. Iai, S., Y. Matsunaga, and T. Kameoka (1990b). Parameter identication for cyclic mobility model, Rep. Port Harbor Res. Institute 29, 5783. Idriss, I. M., and J. I. Sun (1992). Users manual for SHAKE91: a computer program for conducting equivalent linear seismic response analyses of horizontally layered soil deposits, Department of Civil Engineering Technical Report, University of California, Davis. Ishihara, K. (1996). Soil Behavior in Earthquake Geotechnics, Clarendon, London, 350 pp. Ishihara, K., N. Yoshida, and S. Tsujino (1985). Modeling of stress-strain relation of soils in cyclic loading, Fifth International Conference on Numerical Methods in Geomechanics, Nagoya, Japan, 373380. Iwan, W. D. (1967). On a class of models for the yielding behavior of continuous and composite systems, J. Appl. Mech. 34, 612617. Jaeger, J, and N. G. W. Cook (1979). Fundamentals of Rock Mechanics, Third Ed., Chapman and Hall, London, 593 pp. Joyner, W. B. (1999). Equivalent-linear ground-response calculations with frequency-dependent damping in Proc. of 31st Joint Meeting of the US-Japan Panel on Wind and Seismic Effects (UJNR), Tsukuba, Japan, 1114 May, 258264. Joyner, W. B., and D. M. Boore (1998). Equivalent-linear ground response calculations with frequency-dependent damping, in Proc. of 2nd International Symposium on Seismic Hazards and Ground Motion in the Region of Moderate Seismicity, Seoul, Korea, 910 November. Joyner, W. B., and D. M. Boore (2000). Recent developments in earthquake ground motion estimation, in Proc. of Sixth International Conference on Seismic Zonation, Earthquake Engineering Research Institute, Palm Springs, CA, 1215 November 2000.
1629
Yoshida, N., and S. Iai (1998). Nonlinear site response and its evaluation and prediction, in Proc. of Second International Symposium on the Effects of Surface Geology on Seismic Motion, K. Irikura, K. Kudo, H. Okada, and T. Sasatani (Editors), vol. 1, 7190. Yoshida, N., T. Tazoh, and H. Suzuki (1995). Comparison of nonlinear dynamic response analysis methods, in Proc. of 23rd JSCE Earthquake Engineering Symposium, JSCE: 4952 (in Japanese). Yu, G., J. G. Anderson, and R. Siddharthan (1993). On the characteristics of nonlinear soil response, Bull. Seism. Soc. Am. 83, 218244. U.S. Geological Survey Denver Federal Center Box 25046 MS 966 Denver, Colorado 80225 (S.H., R.A.W.) Institute de Radioprotection et de Surete Nucleaire IRSN/DEI/SARG/BERSSIN BP 17-92262 Fontenay-aux-Roses Cedex, France (L.F.B.) Manuscript received 18 June 2004.