Strong Ground Motion Estimation

This progressive sequence of ground motion surprises suggests that the current state of 
knowledge in strong motion seismology is probably not adequate to make unequivocal strong ground motion predictions. However, with these caveats in mind, strong ground motion estimation provides substantial value by reducing risks associated with earthquakes and engineered structures. We present the current state of earthquake ground motion estimation. We start with seismic source characterization, because this is the most important and challenging part of the problem. To better understand the challenges of developing ground motion prediction equations (GMPE) using strong motion data, we present the physical factors that influence strong ground shaking. New calculations are presented to illustrate potential pitfalls and identify key issues relevant to ground motion estimation and future ground motion research and applications. Particular attention is devoted to probabilistic implications of all aspects of ground motion estimation.


Seismic source characterization
The strongest ground shaking generally occurs close to an earthquake fault rupture because geometric spreading reduces ground shaking amplitudes as distance from the fault increases. Robust ground motion estimation at a specific site or over a broad region is predicated on the availability of detailed geological and geophysical information about locations, geometries, and rupture characteristics of earthquake faults. These characteristics are not random, but are dictated by the physical properties of the upper crust including rock types, pre-existing faults and fractures, and strain rates and orientations. Because such information is often not readily available or complete, the resultant uncertainties of source characterization can be the dominant contributions to uncertainty in ground motion estimation. Lettis et al. (1997) showed that intraplate blind thrust earthquakes with moment magnitudes up to 7 have occurred in intraplate regions where often there was no previously known direct surface evidence to suggest the existence of the buried faults. This observation has been repeatedly confirmed, even in plate boundary settings, by numerous large earthquakes of the past 30 years including several which have provided rich sets of ground motion data from faults for which neither the locations, geometries, or other seismic source characterization properties were known prior to the earthquake. Regional seismicity and geodetic measurements may provide some indication of the likely rate of earthquake occurrence in a region, but generally do not demonstrate where that deformation localizes fault displacement. Thus, an integral and necessary step in reducing ground motion estimation uncertainties in most regions remains the identification and characterization of earthquake source faults at a sufficiently detailed scale to fully exploit the full range of ground motion modelling capabilities. In the absence of detailed source characterizations, ground motion uncertainties remain large, with the likely consequence of overestimation of hazard at most locations, and potentially severe underestimation of hazard in those few locations where a future earthquake ultimately reveals the source characteristics of a nearby, currently unknown fault. The latter case is amply demonstrated by the effects of the 1983 M 6.5 Coalinga, 1986M 6.0 Whittier Narrows, 1989M 6.6 Sierra Madre, 1989M 7.0 Loma Prieta, 1992M 7.4 Landers, 1994M 6.7 Northridge, 1999M 7.6 Chi-Chi Taiwan, 2001M 7.7 Bhuj, India, 2010 Canterbury, New Zealand, and 2011 M 6.1 Christchurch, New Zealand, earthquakes. The devastating 2011 M 9.1 Tohoku, Japan, earthquake and tsunami were the result of unusually large fault displacement over a relatively small fault area (Shao et al., 2011), a source characteristic that was not forseen, but profoundly influenced strong ground shaking (NIED, 2011) and tsunami responses (SIAM News, 2011). All these earthquakes occurred in regions where the source faults were either unknown or major source characteristics were not recognized prior to the occurrence of these earthquakes.

Physical basis for ground motion prediction
In this section we present the physical factors that influence ground shaking in response to earthquakes. A discrete representation is used to emphasize the discrete building blocks or factors that interact to produce strong ground motions. For simplicity, we start with linear stress-strain. Nonlinear stress-strain is most commonly observed in soils and evaluated in terms of site response. This is the approach we use here; nonlinear site response is discussed in Section 4. The ground motions produced at any site by an earthquake are the result of seismic radiation associated with the dynamic faulting process and the manner in which seismic energy propagates from positions on the fault to a site of interest. We assume that fault rupture initiates at some point on the fault (the hypocenter) and proceeds outward along the fault surface. Using the representation theorem (Spudich and Archuleta, 1987), ground velocity,   where k is the component of ground motion, ij are the indices of the discrete fault elements, n is the number of fault elements in the strike direction and m is the number of elements in dip direction (Figure 3.1). We use the notation F() to indicate the modulus of the Fourier transform of f(t). It is instructive to take the Fourier transform of (1) and pursue a discussion similar to Hutchings and Wu (1990) and Hutchings (1994) where at each element ij,    S ij  is the source slip-velocity amplitude spectrum,     ij is the source phase spectrum,   G kij  is the Green's function amplitude spectrum, and     kij is the Green's function phase spectrum. The maximum peak ground motions are produced by a combination of factors that produce constant or linear phase variations with frequency over a large frequency band. While the relations in (1) and (2) are useful for synthesizing ground motions, they don't provide particularly intuitive physical insights into the factors that contribute to produce specific ground motion characteristics, particularly large peak accelerations, velocities, and displacements. We introduce isochrones as a fundamental forensic tool for understanding the genesis of ground motions. Isochrones are then used to provide simple geometric illustrations of how directivity varies between dipping dip-slip and vertical strike-slip faults. Bernard and Madariaga (1984) and Spudich and Frazer (1984; developed the isochrone integration method to compute near-source ground motions for finite-fault rupture models. Isochrones are all the positions on a fault that contribute seismic energy that arrives at a specific receiver at the same time. By plotting isochrones projected on a fault, times of large amplitudes in a ground motion time history can be associated with specific regions and characteristics of fault rupture and healing. A simple and reasonable way to employ the isochrone method for sites located near faults is to assume that all significant seismic radiation from the fault consists of first shear-wave arrivals. A further simplification is to use a simple trapezoidal slip-velocity pulse. Let f(t) be the slip function, For simplicity we assume where t r is rupture time, and t h is healing time. Then, all seismic radiation from a fault can be described with rupture and healing isochrones. Ground velocities (v) and accelerations (a) produced by rupture or healing of each point on a fault can be calculated from (Spudich and Frazer, 1984;Zeng et al., 1991;Smedes and Archuleta, 2008)

Isochrones analysis of rupture directivity
x factors can combine to affect ground motions. The arrival times of rupture at a specific receiver are where x is the receiver position, are all fault positions, t  are shear-wave propagation times between the receiver and all fault positions, and t r are rupture times at all fault positions. The arrival times of healing at a specific receiver are where R are the rise times (the durations of slip) at all fault positions. Archuleta (1984) showed that variations in rupture velocity had pronounced effects on calculated ground motions, whereas variations in rise times and slip-rate amplitudes cause small or predictable changes on calculated ground motions. The effect of changing slipvelocity amplitudes on ground motions is strongly governed by the geometrical attenuation (1/r for far-field terms). Any change in the slip-velocity amplitudes affects most the ground motions for sites closest to the region on the fault where large slip-velocities occurred (Spudich and Archuleta, 1987). This is not the case with rupture velocity or rise time; these quantities influence ground motions at all sites. However, as Anderson and Richards (1975) showed, it takes a 300% change in rise time to compensate for a 17% change in rupture time. Spudich and Oppenheimer (1986) show why this is so. Spatial variability of rupture velocity causes the integrand in (3) to become quite rough, thereby adding considerable highfrequency energy to ground motions. The roughness of the integrand in (3)  where T r are the isochrones from (5) and  s is the surface gradient operator. Variations of T r on the fault surface associated with supershear rupture velocities, or regions on the fault where rupture jumps discontinuously can cause large or singular values of c, called critical points by Farra et al. (1986). Spudich and Frazer (1984) showed that the reciprocal of c, isochrone slowness is equivalent to the seismic directivity function in the twodimensional case. Thus, by definition, critical points produce rupture directivity, and as is shown with simulations later, need not be associated strictly with forward rupture directivity, but can occur for any site located normal to a portion of a fault plane where rupture velocities are supershear. It is useful to interpret (3) and (4) in the context of the discrete point-source summations in (1) and (2). When isochrone velocities become large on a substantial area of a fault it simply means that all the seismic energy from that portion of the fault arrives at nearly the same time at the receiver; the summation of a finite, but large number of band-limited Green's functions means that peak velocities remain finite, but potentially large. Large isochrone velocities or small isochrone slownesses over significant region of a fault are diagnostic of ground motion amplification associated with rupture directivity; the focusing of a significant fraction of the seismic energy radiated from a fault at a particular site in a short time interval. In this way isochrones are a powerful tool to dissect ground motions in relation to various characteristics of fault rupture. Times of large ground motion amplitudes can be directly associated with the regions of the fault that have corresponding large isochrone velocities or unusually large slip velocities. From (5) and (6) it is clear that both fault rupture variations, and shear-wave propagation time variations, combine to determine isochrones and isochrone velocities. Boore and Joyner (1989) and Joyner (1991) discussed directivity using a simple line source model. A similar approach is used here to illustrate how directivity differs between vertical strike-slip faults and dipping dip-slip faults. To focus on source effects, we consider unilateral, one-dimensional ruptures in a homogenous half-space (Figure 3.2). The influence of the free surface on amplitudes is ignored. The rupture velocity is set equal to the shearwave velocity to minimize time delays and to maximize rupture directivity. To eliminate geometric spreading, stress drops increase linearly with distance from the site in a manner that produces uniform ground motion velocity contribution to the surface site for all points on the faults. Healing is ignored; only the rupture pulse is considered. Thrust dip-slip faulting is used to produce coincident rake and rupture directions. Seismic radiation is simplified to triangular slip-velocity pulses with widths of one second. For the strike-slip fault, the fault orientation and rupture directional are coincident. But, as fault rupture approaches the site, takeoff angles increase, so the radiation pattern reduces amplitudes, and total propagation distances (rupture length plus propagation distance) increase to disperse shear-wave arrivals in time (Figures 3.2a and 3.2b). The surface site located along the projection of the thrust fault to the surface receives all seismic energy from the fault at the same time, and c is infinity because the fault orientation, rupture, and shearwave propagation directions are all coincident for the entire length of the fault (Figures 3.2c and 2d). Consequently, although the strike-slip fault is 50% longer than the thrust fault, the thrust fault produces a peak amplitude 58% larger than the strike-slip fault. The thrust fault site receives maximum amplitudes over the entire radiated frequency band. High-frequency amplitudes are reduced for the strike-slip site relative to the thrust fault site because shearwaves along the strike-slip fault become increasingly delayed as rupture approaches the site, producing a broadened ground motion velocity pulse. The geometric interaction between dip-slip faults and propagation paths to surface sites located above those faults produces a kinematic recipe for maximizing both isochrone velocities and radiation patterns for surface sites that is unique to dip-slip faults. In contrast, Schmedes and Archuleta (2008) use kinematic rupture simulations and isochrone analyses to show why directivity becomes bounded during strike-slip fault along long faults. Schmedes and Archuleta (2008) consider the case of subshear rupture velocities and use critical point analyses with (3) and (4) to show that for long strike-slip ruptures there is a saturation effect for peak velocities and accelerations at sites close to the fault located at increasing distances along strike relative to the epicenter, consistent with empirical observations (Cua, 2004;Boore and Atkinson, 2008;Campbell and Bozorgnia, 2008;Chiou and Youngs, 2008). Dynamic fault rupture processes during dip-slip rupture complicate dip-slip directivity by switching the region of maximum fault-normal horizontal motion from the hangingwall to the footwall as fault dips increase from 50 to 60 (O'Connell et al., 2007).  (c) and (e) relative to ground motion sites (triangles). Black arrows show the orientation of the faults, red arrows show fault rupture directions, and blue arrows show shear-wave propagation directions (dashed lines) to the sites. Discrete velocity contributions for seven evenly-spaced positions along the fault are shown to the right of each rupture model (b, d, f) as triangles with amplitudes (heights) scaled by the radiation pattern. The output ground motions for each fault rupture are shown in (g). Isochrone velocity, c, is infinity in (d), is large, but finite, in (f), and decreases as the fault nears the ground motion site in (b).

The fundamental difference between strike-slip and dip-slip directivity
Typically, seismic velocities increase with depth, which changes positions of maximum rupture directivity compared to Figure 3.2. For dip-slip faults, the region of maximum directivity is moved away from the projection of the fault to the surface, toward the hanging wall. This bias is dependent on velocity gradients, and the dip and depth of the fault. For strike-slip faults, a refracting velocity geometry can increase directivity by reducing takeoff angle deviations relative to the rupture direction for depth intervals that depend on the velocity structure and position of the surface site (Smedes and Archuleta, 2008). When the two-dimensional nature of finite-fault rupture is considered, rupture directivity is not as strong as suggested by this one-dimensional analysis (Bernard et al., 1996), but the distinct amplitude and frequency differences between ground motions produced by strike-slip and dip-slip faulting directivity remain. Full two-dimensional analyses are presented in a subsequent section. A more complete discussion of source and propagation factors influencing ground motions is presented next to provide a foundation for discussion of amplification associated with rupture directivity. The approach here is to discuss ground motions separately in terms of source and propagation factors and then to discuss how source and propagation factors can jointly interact to strongly influence ground motion behavior.

Seismic source amplitude and phase factors
The flat portion of an amplitude spectrum is composed of the frequencies less than a corner frequency,  c , which is defined as the intersection of low-and high-frequency asymptotes following Brune (1970). The stress drop, , defined as the difference between an initial stress,  0 , minus the dynamic frictional stress,  f , is the stress available to drive fault slip (Aki, 1983). Rise time, R, is the duration of slip at any particular point on the fault. Rise times are heterogeneous over a fault rupture surface. Because the radiation pattern for seismic phases such as body waves and surface waves are imposed by specification of rake (slip direction) at the source and are a function of focal mechanism, radiation pattern is included in the source discussion. Regressions between moment and fault area (Wells and Coppersmith, 1994;Somerville et al., 1999;Leonard, 2010) show that uncertainties in moment magnitude and fault area are sufficient to produce moment uncertainties of 50% or more for any particular fault area. Consequently, the absolute scaling of synthesized ground motions for any faulting scenario have about factor of two uncertainties related to seismic moment (equivalently, average stress drop) uncertainties. Thus, moment-fault area uncertainties introduce a significant source of uncertainty in ground motion estimation. Andrews (1981) and Frankel (1991) showed that correlated-random variations of stress drop over fault surfaces that produce self-similar spatial distributions of fault slip are required to explain observed ground motion frequency amplitude responses. Somerville et al. (1999) showed that a self-similar slip model can explain inferred slip distributions for many large earthquakes and they derive relations between many fault rupture parameters and seismic moment. Their results provide support for specifying fault rupture models using a stochastic spatially varying stress drop where stress drop amplitude decays as the inverse of wavenumber to produce self-similar slip distributions. They assume that mean stress drop is independent of seismic moment. Based on their analysis and assumptions, Somerville et al. (1999) provide recipes for specifying fault rupture parameters such as slip, rise times, and asperity dimensions as a function of moment. Mai and Beroza (2000) showed that 5.3 < M < 8.1 magnitude range dip-slip earthquakes follow self-similar scaling as suggest by Somerville et al. (1999). However, for strike-slip earthquakes, as moment increases in this magnitude range, they showed that seismic moments scale as the cube of fault length, but fault width saturates. Thus, for large strike slip earthquakes average slip increases with fault rupture length, stress drop increases with magnitude, and self-similar slip scaling does not hold. The large stress drops observed for the M 7.7 1999 Chi-Chi, Taiwan thrust-faulting earthquake (Oglesby et al., 2000) suggest that self-similar slip scaling relations may also breakdown at larger moments for dip-slip events.

Factor
Influence Moment rate, Rake and spatial and temporal rake variations scale amplitudes as a function of azimuth and take-off angle. Rake spatial and temporal variations over a fault increase the spatial complexity of radiation pattern amplitude variations and produce frequency-dependent amplitude variability. Rise time, Diffraction at the crack tip introduces a frequency dependent amplitude to the radiation pattern (Madariaga, 1977;Boatwright, 1982;Fukuyama and Madariaga, 1995).
Fault rupture in heterogeneous velocity structure can produce anisotropic slip velocities relative to rupture direction (Harris and Day, 1997) and slip velocities and directivity are a function of rake and dip for dip-slip faults (Oglesby et al., 1998;O'Connell et al., 2007). Frictional heating, fault zone fluids, and melting may also influence radiated energy (Kanamori and Brodsky, 2001;Andrews, 200X High healing velocities increase amplification associated with directivity. Healing velocities interact with stress drop and rise time variations to modify the amplitude spectrum, although to a smaller degree than rupture velocities, since rupture slip veloc ities are typically several times larger than healing slip velocities. Rake, Rake and spatial and temporal rake variations scale amplitudes as a function of azimuth and take-off angle. Rake spatial and temporal variations over a fault increase the spatial complexity of radiation pattern amplitude variations and produce frequency-dependent amplitude variability. Rise time, Diffraction at the crack tip introduces a frequency dependent amplitude to the radiation pattern (Madariaga, 1977;Boatwright, 1982;Fukuyama and Madariaga, 1995).
The same dynamic processes identified in Table 1 produce corresponding phase variability.  Oglesby et al. (1998; showed that stress drop behaviors are fundamentally different between dipping reverse and normal faults. These results suggest that stress drop may be focal mechanism and magnitude dependent. There are still significant uncertainties as to the appropriate specifications of fault rupture parameters to simulate strong ground motions, particularly for larger magnitude earthquakes. O'Connell et al. (2007) used dynamic rupture simulations to show that homogeneous and weakly heterogeneous half-spaces with faults dipping ≲50°, maximum fault-normal peak velocities occurred on the hanging wall. However, for fault dips ≳50°, maximum fault-normal peak velocities occurred on the footwall. Their results indicate that simple amplitude parameterizations based on the hanging wall and/or footwall and the fault normal and/or fault parallel currently used in ground motion prediction relations may not be appropriate for some faults with dips > 50°. Thus, the details of appropriate spatial specification of stress drops and/or slip velocities as a function of focal mechanism, magnitude, and fault dip are yet to be fully resolved. Day (1982) showed that intersonic rupture velocities ( < V r < ) can occur during earthquakes, particularly in regions of high prestress (asperities), and that peak slip velocity is strongly coupled to rupture velocity for non-uniform prestresses. While average rupture velocities typically remain subshear, high-stress asperities can produce local regions of supershear rupture combined with high slip velocities. Supershear rupture velocities have been observed or inferred to have occurred during several earthquakes, including the M 6.9 1979 Imperial Valley strike-slip earthquake (Olson and Apsel, 1982;Spudich and Cranswick, 1984;Archuleta, 1984), the M 6.9 1980 Irpinia normal-faulting earthquake (Belardinelli et al., 1999), the M 7.0 1992 Petrolia thrust-faulting earthquake (Oglesby and Archuleta, 1997), the M 7.3 Landers strike-slip earthquake (Olsen et al., 1997;Bouchon et al., 1998;Hernandez et al., 1999) the M 6.7 1994 Northridge thrust-faulting earthquake (O'Connell, 1999b), and the 1999 M 7.5 Izmit and M 7.3 Duzce Turkey strike-slip earthquakes (Bouchon et al., 2001). Bouchon et al. (2010) find that surface trace of the portions of strike-slip faults with inferred supershear rupture velocities are remarkably linear, continuous and narrow, that segmentation features along these segments are small or absent, and the deformation is highly localized. O'Connell (1999b) postulates that subshear rupture on the faster footwall in the deeper portion of the Northridge fault relative to the hangingwall produced supershear rupture in relation to hangingwall velocities and contributed to the large peak velocities observed on the hangingwall. Harris and Day (1997) showed that rupture velocities and slip-velocity functions are significantly modified when a fault is bounded on one side by a low-velocity zone. The lowvelocity zone can produce asymmetry of rupture velocity and slip velocity. This type of velocity heterogeneity produces an asymmetry in seismic radiation pattern and abrupt and/or systematic spatial variations in rupture velocity. These differences are most significant in regions subject to rupture directivity, and may lead to substantially different peak ground motions occurring at either end of a strike slip fault (Bouchon et al., 2001). Thus, the position of a site relative to the fast and slow sides of a fault and rupture direction may be significant in terms of the dynamic stress drops and rupture velocities that are attainable in the direction of the site. Observations and numerical modeling show that the details of stress distribution on the fault can produce complex rupture velocity distributions and even discontinuous rupture, factors not typically accounted for in kinematic rupture models used to predict ground motions (e.g. Somerville et al., 1991;Schneider et al., 1993;Hutchings, 1994;Tumarkin et al., 1994;Zeng et al., 1994;Beresnev and Atkinson, 1997;O'Connell, 1999c). Even if only smooth variations of subshear rupture velocities are considered (0.6* < Vr < 1.0*), rupture velocity variability introduces ground motion estimation uncertainties of at least a factor of two (Beresnev and Atkinson, 1997), and larger uncertainties for sites subject to directivity. Rupture direction may change due to strength or stress heterogeneities on a fault. Beroza and Spudich (1988) inferred that rupture was delayed and then progressed back toward the hypocenter during the M 6.2 1984 Morgan Hill earthquake. Oglesby and Archuleta (1997) inferred that arcuate rupture of an asperity may have produced accelerations > 1.40 g at Cape Mendocino during the M 7.0 1992 Petrolia earthquake. These results are compatible with numerical simulations of fault rupture on a heterogeneous fault plane. Das and Aki (1977) modeled rupture for a fault plane with high-strength barriers and found that rupture could occur discontinuously beyond strong regions, which may subsequently rupture or remain unbroken. Day (1982) found that rupture was very complex for the case of nonuniform prestress and that rupture jumped beyond some points on the fault, leaving unbroken areas behind the rupture which subsequently ruptured. In the case of slip resistant asperity, Das and Kostrov (1983) found that when rupture began at the edge of the asperity, it proceeded first around the perimeter and then failed inward in a "double pincer movement". Thus, even the details of rupture propagation direction are not truly specified once a hypocenter position is selected. Guatteri and Spudich (1998) showed that time-dependent dynamic rake rotations on a fault become more likely as stress states approach low stresses on a fault when combined with heterogeneous distributions of stress and nearly complete stress drops. Pitarka et al. (2000) found that eliminating radiation pattern coherence between 1 Hz and 3 Hz reproduced observed ground motions for the 1995 M 6.9 Hyogo-ken Nanbu (Kobe) earthquake. Spudich et al. (1998) used fault striations to infer that the Nojima fault slipped at low stress levels with substantial rake rotations occurring during the 1995 Hyogo-ken Nanbu earthquake. This dynamic rake rotation can reduce radiation-pattern coherence at increasing frequencies by increasingly randomizing rake directions for decreasing time intervals near the initiation of slip at each point on a fault, for increasingly complex initial stress distributions on faults. Vidale (1989) showed that the standard double-couple radiation pattern is observable to 6 Hz based on analysis of the mainshock and an aftershock from the Whittier Narrows, California, thrust-faulting earthquake sequence. In contrast, Liu and Helmberger (1985) found that a double-couple radiation pattern was only discernible for frequencies extending to 1 Hz based on analysis the 1979 Imperial Valley earthquake and an aftershock. Bent and Helmberger (1989) estimate a  of 75 MPa for the 1987 Whittier Narrows M 6.1 thrust faulting earthquake, but allow for a  as low as 15.5 MPa. The case of high initial, nearly homogeneous stresses that minimize rake rotations may produce high-frequency radiation pattern coherence as observed by Vidale (1989). These results suggest that there may be a correlation between the maximum frequency of radiation pattern coherence, initial stress state on a fault, focal mechanism, and stress drop. Table 3 lists factors influencing propagation amplitudes, G kij (). Table 4 lists factors influencing propagation phase,  ij . Large-scale basin structure can substantially amplify and extend durations of strong ground motions (Frankel and Vidale, 1992;Frankel, 1993;Olsen and Archuleta, 1996;Wald and Graves;1998;Frankel and Stephenson, 2000;Koketsu and Kikuchi, 2000;Frankel et al., 2001). Basin-edge waves can substantially amplify strong ground motions in basins (Liu and Heaton, 1984;Frankel et al., 1991;Phillips et al., 1993;Spudich and Iida, 1993;Kawase, 1996;Pitarka et al., 1998, Frankel et al., 2001. This is a particular concern for fault-bounded basins where rupture directivity can constructively interact with basin-edge waves to produce extended zones of extreme ground motions (Kawase, 1996;Pitarka et al., 1998), a topic revisited later in the paper. Even smaller scale basin or lens structures on the order of several kilometers in diameter can produce substantial amplification of strong ground motions (Alex and Olsen, 1998;Graves et al., 1998;Davis et al., 2000). Basin-edge waves can be composed of both body and surface waves (Spudich and Iida, 1993;Meremonte et al., 1996;Frankel et al., 2001) which provides a rich wavefield for constructive interference phenomena over a broad frequency range. Critical reflections off the Moho can produce amplification at distances > ~75-100 km (Somerville and Yoshimura, 1990;Catchings and Kohler, 1996). The depth to the Moho, hypocentral depth, direction of rupture (updip versus downdip), and focal mechanism determine the amplification and distance range that Moho reflections may be important. For instance, Catchings and Mooney (1992) showed that Moho reflections amplify ground motions in the > 100 km distance range in the vicinity of the New Madrid seismic zone in the central United States.

Seismic wave propagation amplitude and phase factors
Numerous studies have demonstrated that the seismic velocities in the upper 30 to 60 m can greatly influence the amplitudes of earthquake grounds motions at the surface (e.g. Borcherdt et al., 1979;Joyner et al., 1981;Seed et al., 1988). Williams et al. (1999) showed that significant resonances can occur for impedance boundaries as shallow as 7-m depth.  compared the amplification of generic rock sites with very hard rock sites for 30 m depth averaged velocities. They defined very hard rocks sites as sites that have shear-wave velocities at the surface > 2.7 km/s and generic rock sites as sites where shear-wave velocities at the surface are ~0.6 km/s and increase to > 1 km/s at 30 m depth.  found that amplifications on generic rock sites can be in excess of 3.5 at high frequencies, in contrast to the amplifications of less than 1.2 on very hard rock sites. Considering the combined effect of attenuation and amplification, amplification for generic rocks sites peaks between 2 and 5 Hz at a maximum value less than 1.8 Depending on the dynamic soil properties and pore pressure responses, nonlinear responses can increase or reduce phase dispersion. In the case of coupled pore-pressure with dilatant materials can collapse phase producing intermittent amplification caus tics.
Frequency indepen dent attenuation, Linear hysteretic behavior produces frequency-dependent velocity dispersion that produces frequency dependent phase variations. Scattering, The scattering strength and scattering characteristics determine propagation distances required to randomize the phase of shear waves as a function of frequency. Anisotropy, Complicates shear-wave polarizations and modifies radiation pattern polarizations. Topography, Complicates phase as a function of topographic length scale and near-surface velocities.
A common site-response estimation method is to use the horizontal-to-vertical (H/V) spectral ratio method with shear waves (Lermo and Chavez-Garcia, 1993) to test for site resonances. The H/V method is similar to the receiver-function method of Langston (1979). Several investigations have shown the H/V approach provides robust estimates of resonant frequencies (e.g., Field and Jacob, 1995;Castro et al., 1997;Tsubio et al., 2001) although absolute amplification factors are less well resolved (Castro et al., 1997;. One-dimensional site-response approaches may fail to quantify site amplification in cases when upper-crustal three-dimensional velocity structure is complex. In southern California, Field (2000) found that the basin effect had a stronger influence on peak acceleration than detailed geology used to classify site responses. Hartzell et al. (2000) found that site amplification characteristics at some sites in the Seattle region cannot be explained using 1D or 2D velocity models, but that 3D velocity structure must be considered to fully explain local site responses. Chavez-Garcia et al. (1999) showed that laterally propagating basingenerated surface waves can not be differentiated from 1D site effects using frequency domain techniques such as H/V ratios or reference site ratios. The ability to conduct sitespecific ground motion investigations is predicated on the existence of geological, geophysical, and geotechnical engineering data to realistically characterize earthquake sources, crustal velocity structure, local site structure and conditions, and to estimate the resultant seismic responses at a site. Lack of information about 3D variations in local and crustal velocity structure are serious impediments to ground motion estimation.
It is now recognized that correlated-random 3D velocity heterogeneity is an intrinsic property of Earth's crust (see Sato and Fehler, 1998 for a discussion). Correlated-random means that random velocity fluctuations are dependent on surrounding velocities with the dependence being inversely proportional to distance. Weak (standard deviation, , of ~5%), random fractal crustal velocity variations are required to explain observed short-period (T < 1 s) body-wave travel time variations, coda amplitudes, and coda durations for ground motions recorded over length scales of tens of kilometers to tens of meters (Frankel and Clayton, 1986), most well-log data (Sato and Fehler, 1998), the frequency dependence of shear-wave attenuation (Sato and Fehler, 1998), and envelope broadening of shear waves with distance (Sato and Fehler, 1998). As a natural consequence of energy conservation, the excitation of coda waves in the crust means that direct waves (particularly direct shear waves that dominate peak ground motions) that propagate along the minimum travel-time path from the source to the receiver lose energy with increasing propagation distance as a result of the dispersion of energy in time and space. Following Frankel and Clayton (1986) fractal, self-similar velocity fluctuations are described with an autocorrelation function, P, of the form, where a is the correlation distance, k r is radial wavenumber, n=2 in 2D, and n=3 in 3D. When n=4 an exponential power law results (Sato and Fehler, 1998). Smoothness increasing with distance as a increases in (8) and overall smoothness is proportional to n in (8). This is a more realistic model of spatial geologic material variations than completely uncorrelated, spatially independent, random velocity variations. "Correlated-random" is shortened here to "random" for brevity. Let  denote wavelength. Forward scattering dominates when  << a (Sato and Fehler, 1998). The situation is complicated in self-similar fractal media when considering a broad frequency range relevant to strong motion seismology (0.1 to 10 Hz) because  spans the range  >> a to  << a and both forward and backscattering become important, particularly as n decreases in (8). Thus, it is difficult to develop simple rigorous expressions to quantify amplitude and phase terms associated with wave propagation through the heterogeneous crust (see Sato and Fehler, 1998). O'Connell (1999a) showed that direct shear-wave scattering produced by P-SV-wave coupling associated with vertical velocity gradients typical of southern California, combined with 3D velocity variations with n=2 and a standard deviation of velocity variations of 5% in (8), reduce high-frequency peak ground motions for sediment sites close to earthquake faults. O'Connell (1999a) showed that crustal scattering could substantially influence the amplification of near-fault ground motions in areas subjected to significant directivity. Scattering also determines the propagation distances required to randomize phase as discussed later in this paper. Dynamic reduction of soil moduli and increases in damping with increasing shear strain can substantially modify ground motion amplitudes as a function of frequency (Ishihara, 1996). While there has been evidence of nonlinear soil response in surface strong motion recordings (Field et al., 1997;Cultera et al., 1999), interpretation of these surface records solely in terms of soil nonlinearity is intrinsically non-unique (O'Connell, 1999a). In contrast, downhole strong motion arrays have provided definitive evidence of soil nonlinearity consistent with laboratory testing of soils (Chang et al., 1991;Wen et al., 1995, Ghayamghamain andKawakami, 1996;Satoh et al, 1995Satoh et al, , 1997Satoh et al, , 2001. Idriss and Seed (1968a, b) introduced the "equivalent linear method" to calculate nonlinear soil response, which is an iterative method based on the assumption that the response of soil can be approximated by the response of a linear model whose properties are selected in relation to the average strain that occurs at each depth interval in the model during excitation. Joyner and Chen (1975) used a direct nonlinear stress-strain relationship method to demonstrate that the equivalent linear method may significantly underestimate shortperiod motions for thick soil columns and large input motions. Archuleta et al. (2000) and Bonilla (2000) demonstrated that dynamic pore-pressure responses can substantially modify nonlinear soil response and actually amplify and extend the durations of strong ground motions for some soil conditions. When a site is situated on soil it is critical to determine whether soil response will decrease or increase ground amplitudes and durations, and to compare the expected frequency dependence of the seismic soil responses with the resonant frequencies of the engineered structure(s). When soils are not saturated, the equivalent linear method is usually adequate with consideration of the caveats of Joyner and Chen (1975). When soils are saturated and interbedding sands and/or gravels between clay layers is prevalent, a fully nonlinear evaluation of the site that accounts for dynamic pore pressure responses may be necessary . Lomnitz et al. (1999) showed that for the condition 0.91 1 <  0 , where  1 is the shear-wave velocity of low-velocity material beneath saturated soils, and  0 is the acoustic (compressional-wave) velocity in the near-surface material, a coupled mode between Rayleigh waves propagating along the interface and compressional waves in the near surface material propagates with phase velocity  0 . This mode can propagate over large distances with little attenuation. Lomnitz et al. (1999) note that this set of velocity conditions provides a "recipe" for severe earthquake damage on soft ground when combined with a large contrast in Poisson's ratio between the two layers, and when the resonant frequencies of the mode and engineering structures coincide. Linear 2D viscoelastic finite-difference calculations demonstrate the existence of this wave mode at small strains, but nonlinear 2D finite-difference calculations indicate that long-distance propagation of this mode is strongly attenuated (O'Connell et al., 2010). Anisotropy complicates polarizations of shear waves. Coutant (1996) showed that shallow (< 200 m) shear-wave anisotropy strongly influences surface polarizations of shear waves for frequencies < 30 Hz. Chapman and Shearer (1989) show that quasi-shear (qS) wave polarizations typically twist along ray paths through gradient regions in anisotropic media, causing frequency-dependent coupling between the qS waves. They show that this coupling is much stronger than the analogous coupling between P and SV waves in isotropic gradients because of the small difference between the qS-wave velocities. Chapman and Shearer (1989) show that in some cases, far-field excitation of both quasi-shear wave and shear-wave splitting will result from an incident wave composed of only one of the quasishear waves. The potential for stronger coupling of quasi-shear waves suggests that the influence of anisotropy on shear-wave polarizations and peak ground motion may be significant in some cases. While the influence of anisotropy on strong ground motions is unknown, it is prudent to avoid suggesting that only a limited class of shear-wave polarizations are likely for a particular site based on isotropic ground motion simulations of ground motion observations at other sites. Velocity anisotropy in the crust can substantially distort the radiation pattern of body waves with shear-wave polarization angles diverging from those in an isotropic medium by as much as 90 degrees or more near directions where group velocities of quasi-SH and SV waves deviate from corresponding phase velocities (Kawasaki and Tanimoto, 1981). Thus, anisotropy has the potential to influence radiation pattern coherence as well as ground motion polarization. A common approach is to assume the double-couple radiation pattern disappears over a transition frequency band extending from 1 Hz to 3 Hz (Pitarka et al., 2000) or up to 10 Hz (Zeng and Anderson, 2000). The choice of frequency cutoff for the radiation pattern significantly influences estimates of peak response in regions prone to directivity for frequencies close to and greater than the cutoff frequency. This is a very important parameter for stiff (high-frequency) structures such as buildings that tend to have natural frequencies in the 0.5 to 5 Hz frequency band (see discussion in Frankel, 1999). Topography can substantially influence peak ground motions (Boore, 1972;1973). Schultz (1994) showed that an amplification factor of 2 can be easily achieved near the flanks of hills relative to the flatter portions of a basin and that substantial amplification and deamplification of shear-wave energy in the 1 to 10 Hz frequency range can occur over short distances. Bouchon et al. (1996) showed that shear-wave amplifications of 50% to 100% can occur in the 1.5 Hz to 20 Hz frequency band near the tops of hills, consistent with observations from the 1994 Northridge earthquake (Spudich et al., 1996). Topography may also contribute to amplification in adjacent basins as well as the contributing to differential ground motions with dilatational strains on the order of 0.003 (Hutchings and Jarpe, 1996). Topography has a significant influence on longer-period amplification and groundshaking durations. Ma et al. (2007) showed that topography of the San Gabriel Mountains scatters the surface waves generated by the rupture on the San Andreas fault, leading to lessefficient excitation of basin-edge generated waves and natural resonances within the Los Angeles Basin and reducing peak ground velocity in portions of the basin by up to 50% for frequencies 0.5 Hz or less. These discussions of source and propagation influences on amplitudes and phase are necessarily abbreviated and are not complete, but do provide an indication of the challenges of ground motion estimation, and developing relatively simple, but sufficient ground motion prediction equations based on empirical strong ground motion data. Systematically evaluating all the source and wave propagation factors influencing site-specific ground motions is a daunting task, particularly since it's unlikely that one can know all the relevant source and propagation factors. Often, insufficient information exists to quantitatively evaluate many ground motion factors. Thus, it is useful to develop a susceptibility checklist for ground motion estimation at a particular site. The list would indicate available information for each factor on a scale ranging from ignorance to strong quantitative information and indicate how this state of information could influence ground motions at the site. The result of such a checklist would be a susceptibility rating for potential biases and errors for peak motion and duration estimates of site-specific ground motions.

Nonlinear site response 4.1 Introduction
The near surface geological site conditions in the upper tens of meters are one of the dominant factors in controlling the amplitude and variation of strong ground motion, and the damage patterns that result from large earthquakes. It has long been known that soft sediments amplify the earthquake ground motion. Superficial deposits, especially alluvium type, are responsible for a remarkable modification of the seismic waves. The amplification of the seismic ground motion basically originates from the strong contrast between the rock and soil physical properties (e.g. Kramer, 1996). At small deformations, the soil response is linear: strain and stress are related linearly by the rigidity modulus independently of the strain level (Hooke's law). Mainly because most of the first strong motion observations seemed to be consistent with linear elasticity, seismologists generally accept a linear model of ground motion response to seismic excitation even at the strong motion level. However, according to laboratory studies (e.g. Seed and Idriss, 1969), Hooke's law breaks down at larger strains and the nonlinear relation between strain and stress may significantly affect the strong ground motion at soil sites near the source of large earthquakes. Since laboratory conditions are not the same as those in the field, several authors have tried to find field data to understand nonlinear soil behavior. In order to isolate the local site effects, the transfer function of seismic waves in soil layers has to be estimated by calculating the spectral ratio between the motion at the surface and the underlying soil layers. Variation of these spectral ratios between strong and weak motion has actively been searched in order to detect nonlinearity. For example, Darragh and Shakal (1991) observed an amplification reduction at the Treasure Island soft soil site in San Francisco. Beresnev and Wen (1996) also reported a decrease of amplification factors for the array data in the Lotung valley (Taiwan). Such a decrease has also been observed at different Japanese sites including the Port Island site (e.g. Satoh et al., 1997, Aguirre andIrikura, 1997). On the other hand, Darragh and Shakal (1991) also reported a quasi-linear behavior for a stiff soil site in the whole range from 0.006 g to 0.43g. According to these results there is a need to precise the thresholds corresponding to the onset of nonlinearity and the maximum strong motions amplification factors according to the nature and thickness of soil deposits (Field et al., 1998). Nevertheless, the use of surface ground motion alone does not help to directly calculate the transfer function and these variations. Rock outcrop motion is then usually used to estimate the motion at the bedrock and to calculate sediments amplification for both weak and strong motion (e.g. Celebi et al., 1987;Singh et al., 1988;Darragh et al., 1991;Field et al., 1997;Beresnev, 2002). The accuracy of this approximation strongly depends on near surface rock weathering or topography complexity (Steidl et al., 1996). Moreover, the estimate of site response can be biased by any systematic difference for the path effects between stations located on soil and rock. One additional complication is also due to finite source effects such as directivity. In case of large earthquakes, waves arriving from different locations may interfere causing source effects to vary with site location (Oglesby and Archuleta, 1997). Since these finite source effects strongly depend on the source size, they could mimic the observations cited as evidence for soil nonlinearity. Finally, O'Connell (1999) and Hartzell et al. (2005) show that in the near-fault region of M > 6 earthquakes linear wave propagation in weakly heterogeneous, random three dimensional crustal velocity can mimic observed, apparently, nonlinear sediment response in regions with large vertical velocity gradients that persist from near the surface to several km depth, making it difficult to separate soil nonlinear responses from other larger-scale linear wave propagation effects solely using surface ground motion recordings. Because of these difficulties, the most effective means for quantifying the modification in ground motion induced by soil sediments is to record the motion directly in boreholes that penetrate these layers. Using records from vertical arrays it is possible to separate the site from source and path effects and therefore clearly identify the nonlinear behavior and changes of the soil physical properties during the shaking (e.g. Zeghal and Elgamal, 1994;Aguirre and Irikura, 1997;Satoh et al., 2001;Assimaki et al., 2007;Assimaki et al., 2010;Bonilla et al. 2011).

Nonlinear soil behavior
For years, it has been established in geotechnical engineering that soils behave nonlinearly. This fact comes from numerous experiments with cyclic loading of soil samples. The stressstrain curve has a hysteretic behavior, which produces a reduction of shear modulus as well as an increasing in damping factor.  .1 shows a typical stress-strain curve with a loading phase and consequent hysteretic behavior for the later loading process. There have been several attempts to describe mathematically the shape of this curve, and among those models the hyperbolic is one of the easiest to use because of its mathematical formulation as well as for the number of parameters necessary to describe it (Ishihara, 1996;Kramer, 1996;Beresnev and Wen, 1996) where is the undisturbed shear modulus, and τ is the maximum stress that the material can support in the initial state.
is also known as G max because it has the highest value of shear modulus at low strains. In order to have the hysteretic behavior, the model follows the so-called Masing's rule, which in its basic form translates the origin and expands the horizontal and vertical axis by a factor of 2. Thus, where (γ , τ ) is the reversal point for unloading and reloading curves. This behavior produces two changes in the elastic parameters of the soil. First, the larger the maximum strain, the lower the secant shear modulus obtained as the slope of the line between the origin and the reversal point of the hysteresis loop. Second, hysteresis shows a loss of energy in each cycle, and as it was mentioned above, the energy is proportional to the area of the loop. Thus, the larger the maximum strain, the larger the damping factor. How can the changes in the elastic parameters be detected when looking at transfer functions? We know that the resonance frequencies are proportional to (2 + 1) 4 ⁄ (the fundamental frequency corresponds to = 0). Where is the shear velocity and is the soil thickness. Thus, if the shear modulus is reduced then the resonance frequencies are also reduced because = , where is the material density. In other words, in the presence of nonlinearity the transfer function shifts the resonance frequencies toward lower frequencies. In addition, increased dissipation reduces soil amplification. Figure 4.2 shows an example of nonlinear soil behavior at station TTRH02 (Vs 30 = 340 m/s), KiK-net station that recorded the M JMA 7.3 October 2000 Tottori in Japan. The orange shaded region represents the 95% borehole transfer function computed using events having a PGA less than 10 cm/s 2 . Conversely, the solid line is the borehole transfer function obtained using the data from the Tottori mainshock. One can clearly see the difference between these two estimates of the transfer function, namely a broadband deamplification and a shift of resonance frequencies to lower values. The fact that the linear estimate is computed at the 95% confidence limits means that we are confident that this site underwent nonlinear site responses at a 95% probability level. However, nonlinear effects can also directly be seen on acceleration time histories. Figure 4.3 shows acceleration records, surface and downhole, of the 1995 Kobe earthquake at Port Island (left) and the 1993 Kushiro-Oki earthquake at Kushiro Port (right). Both sites have shear wave velocity profiles relatively close each other, except in the first 30 meters depth. Yet, their response is completely different. Port Island is a man-made site composed of loose sands that liquefied during the Kobe event (Aguirre and Irikura, 1997). Practically there is no energy after the S-wave train in the record at the surface. Conversely, Kushiro Port is composed of dense sands and shows, in the accelerometer located at ground level, large acceleration spikes that are even higher than their counterpart at depth. Iai et al., (1995), Archuleta (1998), and Bonilla et al., (2005) showed that the appearance of large acceleration peak values riding a low frequency carrier are an indicator of soil nonlinearity known as cyclic mobility. Laboratory studies show that the physical mechanism that produces such  phenomenon is the dilatant nature of cohesionless soils, which introduces the partial recovery of the shear strength under cyclic loads. This recovery translates into the ability to produce large deformations followed by large and spiky shear stresses. The spikes observed in the acceleration records are directly related to these periods of dilatancy and generation of pore pressure. These examples indicate that nonlinear soil phenomena are complex. We cannot see the effects of nonlinear soil behavior on the transfer function only, but also on the acceleration time histories. This involves solving the wave equation by integrating nonlinear soil rheologies in the time domain, the subject treated in the next section.

The strain space multishear mechanism model
The multishear mechanism model (Towhata and Ishihara, 1985) is a plane strain formulation to simulate pore pressure generation in sands under cyclic loading and undrained conditions. Iai et al. (1990aIai et al. ( , 1990b modified the model to account for the cyclic mobility and dilatancy of sands. This method has the following strong points:  It is relatively easy to implement. It has few parameters that can be obtained from simple laboratory tests that include pore pressure generation.  This model represents the effect of rotation of principal stresses during cyclic behavior of anisotropically consolidated sands.  Since the theory is a plane strain condition, it can be used to study problems in two dimensions, e.g. embankments, quay walls, among others. In two dimensional cartesian coordinates and using vectorial notation, the effective stress σ′ and strain ϵ tensors can be written as where the superscript T represents the vector transpose operation; σ′ , σ′ , ϵ , and ϵ represent the effective normal stresses and strains in the horizontal and vertical directions; τ and γ are the shear stress and shear strain, respectively. The multiple mechanism model relates the stress and strain through the following incremental equation (Iai et al., 1990a(Iai et al., , 1990b, where the curly brackets represent the vector notation; ϵ is the volumetric strain produced by the pore pressure, and is the tangential stiffness matrix given by The first term is the volumetric mechanism represented by the bulk modulus . The second part is the shear mechanism represented by the tangential shear modulus ( ) idealized as a collection of springs (Figure 4.4). Each spring follows the hyperbolic stress-strain model (Konder and Zelasko, 1963) during the loading and unloading hysteresis process. The shear mechanism may also be considered as a combination of pure shear and shear by differential compression.
In addition, where ∆θ = π I ⁄ is the angle between each spring as shown in Figure 4.4. Towhata and Ishihara (1985) found, using laboratory data, that the pore pressure excess is correlated with the cumulative shear work produced during cyclic loading. Iai et al. (1990aIai et al. ( , 1990b) developed a mathematical model that needs five parameters, called hereafter dilatancy parameters, to take into account this correlation. These parameters represent the initial and final phases of dilatancy, p and p ; overall dilatancy w ; threshold limit and ultimate limit of dilatancy, c and S . These parameters are obtained by fitting laboratory data, from either undrained stress controlled cyclic shear tests or from cyclic stress ratio curves. Details of this constitutive model can be found in Iai et al. (1990aIai et al. ( , 1990b. Fig. 4.4. Schematic figure for the multishear mechanism. The plane strain is the combination of pure shear (vertical axis) and shear by compression (horizontal axis) (after Towhata and Ishihara, 1985).
At this point, this formulation provides only the backbone curve. It is here that the hysteresis is now taken into account by using the generalized Masing rules. In fact, they are not simple rules but a state equation that describes hysteresis given a backbone curve (Bonilla, 2000). They are called generalized Masing rules because its formulation contains Pyke's (1979) and the original Masing 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.

The generalized Masing rules
In previous sections we use the hyperbolic model to describe the stress-strain space of soil materials subjected to cyclic loads. In the hyperbolic model, the nonlinear relation can be written as ⁄ is the reference strain. Introducing the equation above into = , where is the shear stress and is the shear strain; and adding the hysteresis operator, we have where is the backbone curve, and (. ) is the hysteresis operator (application of the generalized Masing rules). Hysteresis behavior can be implemented in a phenomenological manner with the help of the Masing and extended Masing rules (Vucetic, 1990;Kramer, 1996). However, these rules are not enough to constrain the shear stress τ to values not exceeding the strength τ . This happens when the time behavior of the shear strain departs from the simple cyclic behavior, and of course, noncyclic time behavior is common in seismic signals. Inadequacy of the Masing rules to describe the hysteretic behavior of complicated signals has been already pointed out and some remedies have been proposed (e.g. Pyke, 1979;Li and Liao, 1993). The Masing rules consist of a translation and dilatation of the original law governing the strain-stress relationship. While the initial loading of the material is given by the backbone curve ( ), the subsequent loading and unloading, the strain-stress relationship is given by: where the coordinate ( , ) corresponds to the reversal points in the strain-stress space, and is the so-called hysteresis scale factor . In Masing's original formulation, the hysteresis scale factor is equal to 2. A first extension to the Masing rules can be obtained by releasing the constraint = 2. This parameter controls the shape of the loop in the stress-strain space (Bonilla et al., 1998). However, numerical simulations suggest spurious behavior of for irregular loading and unloading processes even when extended Masing rules are used. A further generalization of Masing rules is obtained choosing the value of in such way to assure that the path , at a given unloading or reloading, in the strain-stress space will cross the backbone curve, and becomes bounded by the maximum strength of the material . This can be achieved by having the following condition, The equation above represents a general constraint on the hysteresis scale factor, so that the computed stress does not exceed depending on the chosen maximum deformation that the material is thought to resist. The limit → ∞ corresponds to the Cundall-Pyke hypothesis (Pyke, 1979), while → is similar to some extent to a method discussed in (Li and Liao, 1993). In the following section, we will see an example of application of this soil constitutive model (Towhata and Ishihara, 1985;Iai et al., 1990aIai et al., , 1990b together with the Generalized Masing hysteresis operator (Bonilla, 2000).

Analysis of the 1987 Superstition Hills Earthquake
On 24 November 1987, the M L 6.6 Superstition Hills earthquake was recorded at the Wildlife Refuge station. This site is located in southern California in the seismically active Imperial Valley. In 1982 it was instrumented by the U.S. Geological Survey with downhole and surface accelerometers and piezometers to record ground motions and pore water pressures during earthquakes (Holzer et al., 1989). The Wildlife site is located in the flood plain of the Alamo River, about 20 m from the river's western bank. In situ investigations have shown that the site stratigraphy consists of a shallow silt layer approximately 2.5 m thick underlain by a 4.3 m thick layer of loose silty sand, which is in turn underlain by a stiff to very stiff clay. The water table fluctuates around 2-m depth (Matasovic and Vucetic, 1993). This site shows historically one direct in situ observation of nonlinearity in borehole data. The Wildlife Refuge liquefaction array recorded acceleration at the surface and 7.5-m depth, and pore pressure on six piezometers at various depths (Holzer et al., 1989). The acceleration time histories for the Superstition Hills events at GL-0 m and GL-7.5 m, respectively, are shown in Figure 4.5 (left). Note how the acceleration changes abruptly for the record at GL-0 m after the S wave. Several sharp peaks are observed; they are very close to the peak acceleration for the whole record. In addition, these peaks have lower frequency than the previous part of the record (the beginning of the S wave, for instance). Zeghal and Elgamal (1994) used the Superstition Hills earthquakes to estimate the stress and strain from borehole acceleration recordings. They approximated the shear stress (ℎ, ) at depth ℎ, and the mean shear strain ̅ between the two sensors as follows, where (0, ) is the horizontal acceleration at the ground surface, (ℎ, ) is the acceleration at depth ℎ (evaluated through linear interpolation); ( , ) is the acceleration at the bottom of the layer; ( , ) and (0, ) are the displacement histories obtained by integrating twice the corresponding acceleration histories; is the thickness of the layer; and is the density. Using this method, the stress and strain at GL-2.9 m were computed (Figure 4.5). This figure clearly shows the large nonlinearity developed during the Superstition Hills event. The stress-strain loops form an S-shape and the strains are as large as 1.5%. At this depth, there is a piezometer (P5 according to Holzer et al., 1989). With this information it is also possible to reconstruct the stress path (bottom right of Figure 4.5). Note that some of the pore pressure pulses are correlated with episodes of high shear stress development. The stress path shows a strong contractive phase followed by dilatancy when the effective mean stress is close to 15 kPa.  Zeghal and Elgamal (1994), stress-strain loops and stress path history reconstitution (right).
Using the stress and stress time histories at GL-2.9 m computed earlier, Bonilla et al. (2005) performed a trial-and-error procedure in order to obtain the dilatancy parameters that best reproduce such observations. Figure 4.6 compares the computed shear stress time history with the observed shear strain at GL-2.9 m. The stress-strain hysteresis loops are also shown. We observe that the computed shear stress is well simulated; the stress-strain space also shows the same dilatant behavior (S-shape hysteresis loops) as the observed data.
Once the model parameters were determined, they proceed to compute the acceleration time history at GL-0 m using the north-south record at GL-7.5 m as input motion. Figure 4.7 shows the accelerograms (left) and the corresponding response spectra (right). The observed data are shown with no filtering, whereas the computed data are low-pass filtered at 10 Hz.
The computed accelerogram shows the transition from high-frequency content between 0 and 15 sec to the intermittent spiky behavior after 15 sec. The response spectra show that the computed accelerogram accurately represents the long periods; yet, the short periods are still difficult to model accurately. This is the challenge of nonlinear simulations; the fit should be as broadband as possible.

NOAH2D 2D P-SV analyses of the maximum observed peak acceleration
Current nonlinear formulations generally reproduce all first-order aspects of nonlinear soil responses. To illustrate this point, we present a nonlinear analyses of the largest peak ground acceleration recorded to date of > 4 g, that has a peak vertical acceleration of 3.8 g (Aoi et al., 2008). Aoi et al. (2008) analyzed ground motions recorded by the Kyoshin Network (Kik-net) during the M 6.9 2008 Iwate-Miyagi earthquake that included one soilsurface site that recorded a vertical acceleration of 3.8g (station IWTH25). The horizontal borehole and surface motions reported in Aoi et al. (2008) for station IWTH25 are generally consistent with the soil reducing surface horizontal accelerations at high frequencies, as is widely observed at soil sites (Field et al., 1997;1998;Seed and Idriss, 1970;Beresnev et al., 1995Beresnev et al., , 2002 and reproducible using 1D nonlinear site response modeling (O'Connell, 2008). However, the surface vertical peak acceleration exceeded 3.8g, exceeding the maximum expected amplification, based on the site velocity profile between the borehole and the surface accelerometers, and current 1D linear or nonlinear theories of soil behavior (O'Connell, 2008). In particular, application of the nonlinear approach of shear-modulus reduction advocated and tested by Bersenev et al. (2002) to predict nonlinear vertical responses, failed to predict peak vertical accelerations in excess of 2g (O'Connell, 2008). Further, Aoi et al. (2008) observed largest upward accelerations at the surface that were 2.27 times larger than the largest downward accelerations, a result not reproduced using 1D approaches to approximate soil nonlinearity. The 2D nonlinear wave propagation implementation of Bonilla et al. (2006) uses a plane-strain model (Iai et al., 1990a(Iai et al., , 1990b. In this section we show that this model could explain the first-order soil responses observed at station IWTH25 using fairly generic approximation to the site's nonlinear soil properties. The P-SV nonlinear rheology developed by Iai et al. (1990aIai et al. ( , 1990b was used in the Bonilla et al. (2006) implementation of 2D nonlinear wave propagation. The constitutive equation implemented corresponds to the strain space multishear mechanism model developed by Towhata and Ishihara (1985) and Iai et al. (1990aIai et al. ( , 1990b with its backbone characterized by the hyperbolic equation (Hardin and Drnevich, 1972). The multishear mechanism model is a plane strain formulation to simulate cyclic mobility of sands under undrained conditions. In the calculations of this study, a total stress rheology (pore pressure was ignored) was used in the second order staggered grid P-SV plane-strain finite difference code. Perfectly matched layer (PML) absorbing boundary conditions were used to approximate elastic (transmitting) boundary conditions at the bottom and side edges using an implementation adapted for finite differences from Ma and . Linear hysteretic damping (Q) was implemented using the method of Liu and Archuleta (2006). The geologic environment at station IWTH25 will clearly produce lateral changes in shallow velocity structure. In particular, the hangingwall uplift associated with repeated faulting similar to the 2008 earthquake will produce a series of uplift terraces adjacent to the stream next to station IWTH25, with the lowest shallow velocities being found on the lowest terrace adjacent to the stream, where station IWTH25 is located. The width of the stream and lowest terrace is about 100m near station IWTH25. We constructed a 2D velocity model by including a region 100 m wide with surface Vs=300 m/s layer 2-m deep and then extended Vs=500 m/s to the free surface in the region surrounding the 100-m-wide low-velocity surface layer. Station IWTH25 is assumed to be located relatively close (4-5 m) to the lateral velocity change within the lowst-velocity portion of the 2D velocity model because the geologic log from station IWTH25 indicates a only 1-2 m of young terrace deposits (Aoi et al., 2008), but the youngest terrace probably extends across and encompasses the stream channels and their margins. The dominant large-amplitude arrivals in the borehole motions are associated with large slip regions below and just south of station IWTH25. Consequently, a planewave incident at 80 degrees from the south was used to propagate the borehole motion to the surface in the 2D model.
The nonlinear properties were simplified to a depth-independent plasticity index (PI) of 20% for the NOAH2D calculations. Overall the 2D synthetic nonlinear horizontal motions provide a good fit to the acceleration response spectra (Figs. 4.8a and 4.8d) and acceleration seismograms (Figs. 4.8b and 4.8e). The 2D synthetic horizontal velocities match the observed velocity seismograms well, except in the early portion of the record where the translation ("fling") associated with permanent displacement that dominates early portions of the observed seismograms (Figs. 4.8c and 4.8f). Synthetic vertical responses were calculated for each horizontal-vertical component pair which is a crude approximation to total 3D wavefield. The east component is nearly fault-normal and has the largest peak accelerations and velocities of the two horizontal components, so the eastvertical combination probably best corresponds to the dominant P-SV responses. Except for the obvious asymmetry in both the acceleration and velocity vertical seismograms, both the north-vertical and east-vertical 2D nonlinear synthetic vertical surface motions provide a good fit to the observed acceleration response spectra (Figures 4.9a and   It is important to mention some factors that are not explicitly accounted for in the approach of Bonilla et al. (2006). Goldberg (1960) was among the first to theoretically show the interaction between P and S waves in an elastic medium for large-amplitude seismic waves. His solution yielded the following results: (1) P-and S-waves couple, (2) S waves induce P waves, (3) the induced waves have a dominant frequency twice the Swave frequency, (4) the induced P waves propagate ahead with P-wave velocity. Loukachev et al. (2002) develop a mechanism to explain the frequency dispersion in the Pwave spectra relative to S-wave spectra due to the interaction between compressional (P) and shear (S) waves in granular materials. Shear waves induce dilatancy and contractancy in granular materials which produces longitudinal dilatancy waves (so-called D waves) with approximately double frequency. The results of Goldberg (1960) are sufficient to partially explain the frequency doubling of the dominant vertical resonance frequency even if the IWTH25 material does not exhibit granular dilatancy and is dominated by relatively stiff clay responses. However, it appears more likely that the IWTH25 responses reflect a combination of large-amplitude P-S-wave linear coupling combined with the nonlinear dilatant material responses predicted by Loukachev et al. (2002) for granular soils. This may explain the broad frequency bandwidth (9-18 Hz) of large vertical spectral response. Even with these caveats, the synthetic results illustrated in Figures 9-12 show that the approach of Bonilla et al. (2006) clearly provides a sound basis to evaluate firstorder nonlinear horizontal and vertical response at shallow soil sites and deeper soil basins.

Ground motion prediction equations based on empirical data
Ground motion observations are the result of a long history of instrument development and deployment, instigated primarily by earthquake engineers, to acquire data to develop an empirical foundation to understand and predict earthquake ground motions for use in the design of engineered structures. Strong motion instruments usually record time histories of ground acceleration that can be post-processed to estimate ground velocities and displacements. A particularly useful derived quantity for engineering analyses are response spectra, which are the maximum amplitudes of modestly damped resonant responses of single-degree-of-freedom oscillators (an idealization of simple building responses) to a particular ground motion time history, as a function of natural period or natural frequency. While peak accelerations are always of concern for engineering analyses, peak ground velocity is now recognized as a better indicator of damage potential for large structures than is peak ground acceleration (EERI, 1994). Engineering analyses often consist of linear approaches to determine if structures reach their linear strength limits. Ground motion estimation quantities required for linear analyses are peak accelerations and velocities and associated response spectra. Nonlinear engineering analyses require estimates of future acceleration time histories. The discussion presented in this section focuses on empirical ground motion parameter estimation methods. Ground motion estimation methods required for nonlinear engineering analyses are presented in subsequent sections.
Historically the estimation of ground motion parameters such as peak acceleration, velocity, and displacement, and response spectral ordinates, and duration has been based on regression relationships developed using strong motion observations. These ground motion prediction equations strive to interpolate and extrapolate existing ground motion measurements to serve the needs to design for seismic loads.

Function form of GMPEs for regression
In their simplest form, these empirical GMPEs predict peak ground motions based on a limited parametric description of earthquake and site characteristics. Peak ground motion amplitudes generally increase with increasing magnitude up to a threshold magnitude range where peak accelerations saturate, i.e., only slightly increase or stay nearly constant above the threshold magnitude range (Campbell, 1981;Boore et al., 1994). Similarly, observed peak ground motion amplitudes decrease with increasing distance from the earthquake fault, but saturate at close distances to faults such that the decrease in amplitudes with increasing distance is small within several km of faults. These GMPEs relate specific ground motion parameters to earthquake magnitude, reduction (attenuation) of ground motion amplitudes with increasing distance from the fault (geometric spreading), and local site characteristics using either site classification schemes or a range of quantitative measures of shallow to deeper velocity averages or thresholds. The 30-m-average shear-wave velocity (Vs30) is most commonly used to account for firstorder influences of shallow site conditions. Depths to shear-wave velocities of 1.0, 1.5, and 2.5 km/s (Z1.0 in  and Chiou and Youngs (2008), Z1.5 in  and Day et al. (2008), and Z2.5 in Campbell and Bozorgnia (2008), respectively) are sometimes used to account for influences of larger scale crustal velocity structure on ground motions. The "Next Generation Attenuation" (NGA) Project was a collaborative research program with the objective of developing updated GMPEs (attenuation relationships) for the western U.S. and other worldwide active shallow tectonic regions. These relationships have been widely reviewed and applied in a number of settings (Stafford et al., 2008;Shoja-Taheri et al., 2010). Five sets of updated GMPEs were developed by teams working independently but interacting throughout the NGA development process. The individual teams all had previous experience in the development of GMPEs. The individual teams all had access to a comprehensive, updated ground motion database that had been consistently processed . Each team was free to identify portions of the database to either include or exclude from the development process. A total of 3551 recordings were included in the PEER-NGA database. The number of records actually used by the developers varied from 942 to 2754. The individual GMPEs are described in , Boore and Atkinson (2008), Campbell and Bozorgnia (2008), Chiou and Youngs (2008), and Idriss (2008). These models are referred to as AS08, BA08, CB08, CY08, and I08, respectively, below. The NGA GMPEs developed equations for the orientation-independent average horizontal component of ground motions (Boore et. al., 2006). The NGA account for these ground motion factors using the general form, and, where Y is the ground motion parameter of interest (peak acceleration, velocity, displacement, response spectral ordinate, etc.), M is magnitude, R is a distance measure, M REF and C SOURCE are magnitude and distance terms that define the change in amplitude scaling and the F [source, site, HW, main] are indicator variables of source type, site type, hanging wall geometry and main shock discriminator. The A i are coefficients to be determined by the regression. Not all of the five NGA GMPEs utilize all of these F indicator variables. The  lnY term represents the estimate of the period dependent standard deviation in the parameter lnY at the magnitude and distance of interest. The NGA models use different source parameters and distance measures. Some of the models include the depth to top of rupture (TOR) as a source parameter. This choice was partially motivated by research  that suggested a systematic difference in the ground motion for earthquakes with buried ruptures producing larger short period ground motions as compared to earthquakes with surface rupture. Large reverse-slip earthquakes tend to be buried ruptures more often than large strike-slip earthquakes so the effect of buried ruptures may be partially incorporated in the style-offaulting factor. Not all the NGA developers found the inclusion of TOR to be a statistically significant factor. All of the models except for I08 use the time-averaged S-wave velocity in the top 30 m of a site, Vs30, as the primary site response parameter. I08 is defined only for a reference rock outcrop with Vs30 = 450-900 m/s. Approximately two thirds of the recordings in the PEER-NGA database were obtained at sites without measured values of shear-wave velocity. Empirical correlations between the surface geology and Vs30, were developed (Chiou and others, 2008) and used with assessments of the surface geology to estimate values of Vs30 at the sites without measured velocities. The implications of the use of estimated Vs30 on the standard deviation ( T ) was evaluated and included by AS08.
All of the relationships that model site response incorporate nonlinear site effects. Two different metrics for the strength of the shaking are used to quantify nonlinear site response effects. AS08, BA08, and CB08 use the median estimate of PGA on a reference rock outcrop in the nonlinear site response term. CY08 uses the median estimate of spectral acceleration on a reference rock outcrop at the period of interest. The definition of "reference rock" varies from Vs30=535 m/s (I08) to Vs30=1130 m/s (CY08). A very small fraction of the strong-motion data in the PEER-NGA data set was obtained at sites with Vs30> 900m/s. Depths to shear-wave velocities of 1.0, 1.5, and 2.5 km/s (Z1.0 in  and Chiou and Youngs (2008), Z1.5 in  and Day et al. (2008), and Z2.5 in Campbell and Bozorgnia (2008), respectively) are sometimes used to account for influences of larger scale crustal velocity structure on ground motions. The implications of the methodology chosen to represent larger scale crustal velocity structure on ground motions is discussed in more detail below. The standard deviation or aleatory variability, often denoted as sigma (σ T ), exerts a very strong influence on the results of probabilistic seismic hazard analysis (PSHA) (Bommer and Abrahamson 2006). For this reason it is important to note that the total aleatory uncertainties, as well as the intra-and inter-event uncertainties are systematically larger for the new NGA equations relative to previous relationships Sadigh et al., 1997;Campbell, 1997). Three of the NGA models incorporate a magnitude dependence in the standard deviation. For magnitudes near 7, the five NGA models have similar standard deviations. However, for M < 5.5, there is a large difference in the standard deviations, with the three magnitude-dependent models exhibiting much larger standard deviations ( T > 0.7) than the magnitude-independent models ( T ~0.54). The three models that include a magnitude-dependent standard deviation (AS08, CY08 and I08) all included aftershocks, whereas the two models that used a magnitude-independent standard deviation (BA08 and CB08) excluded them. Including aftershocks greatly increases the number of small-magnitude earthquakes. However, there is a resulting trade-off of significantly larger variability in predicted ground motions than if only large magnitude mainshocks are used. Significant differences in the standard deviations are also noted for soil sites at short distances, this is most likely due to inclusion or exclusion of nonlinear site effects on the standard deviation. In general, the NGA models predict similar median values (within a factor of ~1.5) for vertical strike-slip earthquakes with 5.5 < M < 7.5. The largest differences are for small magnitudes (M < 5.5), for very large magnitudes (M = 8), and for sites located over the hanging wall of dipping faults . As more data has become available to the GMPE developers the number of coefficients in the relationships has increased significantly (>20 in some cases). However, the aleatory variability values ( T ) have not decreased through time (J. Bommer, pers. comm.). Since empirical GMPEs, including NGA GMPEs, are by necessity somewhat generic compared to the wide range of seismic source, crustal velocity structure, and site conditions encountered in engineering applications, there are cases when application of empirical GMPEs is difficult and most importantly, more uncertain. In the context of PSHA, these additional epistemic (knowledge) uncertainties, when quantified, are naturally incorporated into the probabilistic estimation of ground motion parameters. We present two situations of engineering interest, where the application of empirical GMPEs is challenging, to illustrate the difficulties and suggest a path forward in the ongoing process to update and improve empirical GMPEs.

Application of NGA GMPEs for near-fault Vs30 > 900 m/s sites
Independent analysis of the performance of the NGA GMPEs with post-NGA earthquake ground motion recordings demonstrate that use of measured site Vs30 characteristics leads to greatly improved ground motion predictions, with lower performance for sites where Vs30 is inferred instead of directly measured . Thus, the use of Vs30 represents a significant improvement over previous generations of GMPEs that use a simple qualitative site classification scheme. Kaklamanos and Baise (2011) suggest that development of of better site characteristics than Vs30 may also improve the prediction accuracy of GMPEs. In this section we illustrate the challenges presented in the use and application of Vs30 in the NGA GMPE regressions and application of the NGA GMPEs for "rock" sites. It is becoming more common to need ground motion estimates for "rock" site conditions to specify inputs for engineering analyses that include both structures and shallow lowervelocity materials within the analysis model. In this section we consider the challenges in estimating ground motions for site conditions of Vs30 > 900 m/s close to strike slip faults. The problem is challenging for empirical GMPEs because most of the available recordings of near-fault strike-slip ground motions are from sites with Vs30 on the order of 300 m/s. The NGA GMPEs that implement Vs30 used empirical and/or synthetic amplification functions that involve modifying the observed ground motions prior to regression. In this section we discuss some of the challenges of this approach as it applies to estimating ground motions at rock (Vs30 > 900 m/s) sites that are typical of foundation conditions for many large and/or deeply embedded structures. The four NGA GMPEs that implement Vs30 using deterministic ("constrained") amplification coefficients to remap the observed near-fault strike-slip strong motion data that have an average Vs30=299 m/s PSA (Table 5) prior to regression. In contrast,  applied non-linear multi-stage regression using the observed data directly; the observed ground motion values were directly employed in their regression with no remapping of values due to site characteristics). Boore and Atkinson (2008) used the  linear amplification coefficients to remap observed response spectra to a reference Vs30=760 m/s. Campbell and Bozorgnia (2008) used 1D nonlinear soil amplification simulation results of Walling et al. (2008) to deterministically fix nonlinear amplification and remap all response spectra with Vs30 < 400-1086 m/s, depending on period, to create the response spectral "data" input into the non-linear multi-stage regression.  use an approach similar to Campbell and Bozorgnia (2008). Chiou and Youngs (2008) do not explicitly specify how the coefficients for linear and nonlinear amplification were constrained or obtained. Thus, Boore and Atkinson (2008) remap observed response spectra prior to regression using the linear coefficients from , Campbell and Bozorgnia (2008) and  remap observed response spectra prior to regression using the nonlinear coefficients from Walling et al. (2008), and it is not clear what Chiou and Youngs (2008) did. We use Boore and Atkinson (2008) and Campbell and Bozorgnia (2008) to illustrate how the observed response spectral data for sites with Vs30=300 m/s are changed to create the actual "data" used in the regression to estimate Vs30=915 m/s ground motions. It is instructive to compare the approaches and resulting near-fault ground motion predictions. For Boore and Atkinson (2008) amplification normalized by the GMPEs longest period amplification (10 s) is used in Fig. 5.1, to clearly illustrate the scale of the a priori deterministic linear amplification as a function of period. The a priori deterministic linearamplification normalization (Fig. 5.1a) takes the original median near-fault response spectra that have a peak amplitude at about 0.65 s ( Figure 5.1) and create response spectra with peak amplitude at 0.2 s that is used as the "observed data" (red curve in Figure 5.1b) in the nonlinear multi-stage GMPE regression. For Campbell and Bozorgnia (2008), the nonlinear Vs30 amplification coefficients are fixed and create the deterministic nonlinear amplification function ( Figure 5.2a) that is always applied to Vs30 < 400 m/s PSA at all periods to create the "data" (red curve in Figure 5.2b) used in the non-linear multi-stage GMPE regression. In the case of nonlinear deterministic amplification it is necessary to specify a reference PGA. We use 0.45 g for the reference PGA for illustration since this is close to the median ground motion case for sites about 2 km from strike-slip faults and M > 6.; use of a higher reference PGA would increase the nonlinear amplification in Fig. 5.2a.

Earthquake
The use of a single deterministic amplification function for Vs30, whether linear or nonlinear, assumes that there is a one-to-one deterministic mapping of period-dependent amplification to Vs30, which Idriss (2008) suggests is not likely; a single Vs30 can be associated with a wide variety of amplification functions. Further, in the case of nonlinear amplification Abrahamson and, a single deterministic nonlinear amplification function used to account for modulus reduction and damping that vary widely as a function of soil materials, as discussed in Section 4 and Bonilla et al. (2005).    to have a peak acceleration response at about 0.2 s, prior to regression. Thus, in hindsight it may not be a surprise that the NGA response spectra maintain a strong bias to peak at 0.2 s period that in large part is the result of the deterministic amplification modifications to the observed data prior to non-linear multi-stage regression. What is remarkable is that all four NGA GMPEs that implement Vs30 and Idriss (2008) predict that spectral accelerations normalized by peak ground acceleration always peak at about 0.2 s, virtually independent of magnitude for M > 6; the overall shape of Boore and Atkinson (2008) response spectra normalized by peak ground acceleration in Fig. 5.3a are representative of all NGA GMPE response spectral shapes in terms of overall spectra shape and the 0.2 s period of maximum response.  obtained a quite different result, with the period of peak spectral amplitude shifting to longer periods as magnitude increases above M 6.6 ( Fig. 5.3b). The few near-fault data from sites with Vs30 > 900 m/s (Table 6 and Fig. 5.3c) are more consistent with the  results, with the median period of maximum spectral acceleration of 0.38 s being nearly twice that predicted by the four NGA GMPEs (Fig. 5.3c). Mostly, important for a wide class of larger buildings and critical structures, the broader longer-period spectral shift of  suggests much larger (> 50% larger spectral responses at periods > 0.4 s (Fig. 5.3c) Based on current seismological knowledge and theory, what should we expect?  Table 6 is compared to the mean Boore and Atkinson (2008) and  estimates in (c).
Ground motion acceleration at high frequency scales in proportion to dynamic stress drop (Boore, 1983). Average slip is proportional to the product of average dynamic stress drop and average rise time. Dynamic stress drop averaged over the entire fault plane is generally found to remain relatively constant with magnitude (Aki, 1983;Shaw, 2009). Thus, as average slip increases with magnitude (Somerville et al., 1999;Mai and Beroza, 2000;Mai et al., 2006) average rise time must also increase with increasing magnitude. Somerville (2003) notes that the period of the dominant amplitude near-fault motions is related to source parameters such as the rise time and the fault dimensions, which generally increase with magnitude. Mai et al. (2006) present an analysis of scaling of stress drop with seismic moment and find a strong increase of maximum stress drop on the fault plane as a function of increasing moment. In contrast, average stress drop over the entire fault plane at most only slightly increases with increasing moment; the substantial scatter of average stress drop values in Figure 1 of Mai et al. (2006) is consistent with average stress drop that is constant with moment. The Mai et al. (2006) results for maximum stress drop are consistent with first-order constraints on stochastic aspects of seismic source properties (Andrews, 1981;Boore, 1983;Frankel, 1991). As fault area increases, the probability of observing a larger stress drop somewhere on the fault plane increases since stress drop must exhibit correlatedrandom variability over the fault to explain the first-order observations of seismic source properties inferred from ground motion recordings, such as the  2 spectral shape (Andrews, 1981;Frankel, 1991). However, for the moment range (6.5 < M < 7.5) that dominates the hazard at many sites, the stress drop averaged over the entire fault plane is generally found to remain relatively constant with magnitude (Aki, 1983;Shaw, 2009), thus requiring average rise time to increase with increasing magnitude. These fundamental seismological constraints derived from analyses of many earthquakes require that the period that experiences peak response spectral amplitudes should increase with magnitude for some threshold magnitude. The results of the  GMPEs suggest the threshold magnitude is about M 6.6 ( Fig. 5.3b). That all five NGA GMPEs predict invariance of the period of peak spectral response amplitude for M > 6.6 to M 8.0 (example from M 6.6 to M 7.4 are shown in Fig. 5.3a) implies that stress drop increases strongly with increasing magnitude, which is inconsistent with current knowledge of seismic source properties. In contrast, the  response-spectral-magnitude-period-dependent results are more consistent with available seismological constraints. It is important to understand why.  implements Vs30 site factors in a quite different manner than the four NGA GMPEs that use Vs30.  applied non-linear multi-stage regression using the observed data directly, with no deterministic remapping of data by Vs30 prior to regression. Except for their deterministic treatment of Vs30, Boore and Atkinson (2008) use a similar regression approach for . Since Boore and Atkinson (2008) regress period-by-period, the linear site-response remapping (Fig. 5.1a) effectively swamps any signal associated with a period shift with increasing magnitude observed by ; a non-linear regression will operate on the largest signals. The deterministic linear amplification function in Boore and Atkinson (2008) becomes a very large signal ( Figure  5.1a) when operating on data from Vs30=300 m/s sites. The other NGA GMPE regressions normalize response spectra by peak ground acceleration prior to regression, which  suggest tends to reduce resolution of the period-amplitude response-spectra variations in multi-stage regression. Figs 5.1 and 5.2 illustrate why the NGA GMPEs predict PSA shapes that barely change with magnitude ( Fig. 5.3a) and why the NGA GMPEs do not match the first-order characteristics of M > 6.6 near-fault PSA (Fig. 5.3c). It simply might be true that once nonlinear amplification occurs it is impossible to resolve differences between period shifts associated with source processes and site responses. Yet, implicitly the NGA GMPE non-linear regressions assume resolution of all possible response-spectral shape changes as a function of magnitude using deterministic site response amplification functions, an assumption Idriss (2008) does not find credible. In contrast,  used the actual unmodified response spectral data in their multi-stage regression and obtained results compatible with existing seismological constraints. Unfortunately, this leaves us in a bit of a conundrum based on GMPE grading criteria suggested by Bommer et al. (2010) and Kaklamanos and Baise (2011), which clearly establish that NGA is a significant improvement for a wide range of applications than previous generation GMPEs, including . A primary contributor to this conundrum about appropriate spectral behaviour for near-fault Vs30 > 900 m/s sites is the lack of near-fault ground motion data for Vs30 > 900 m/s (Table 6 and Figure 5.3c), providing a vivid real-world example of epistemic uncertainty. The site amplification approach used in NGA is discussed by Boore and Atkinson (2008), "The rationale for pre-specifying the site amplifications is that the NGA database may be insufficient to determine simultaneously all coefficients for the nonlinear soil equations and the magnitude-distance scaling, due to trade-offs that occur between parameters, particularly when soil nonlinearity is introduced. It was therefore deemed preferable to "hard-wire" the soil response based on the best-available empirical analysis in the literature, and allow the regression to determine the remaining magnitude and distance scaling factors. It is recognized that there are implicit trade-offs involved, and that a change in the prescribed soil response equations would lead to a change in the derived magnitude and distance scaling. Note, however, that our prescribed soil response terms are similar to those adopted by other NGA developers who used different approaches; thus there appears to be consensus as to the appropriate level for the soil response factors." This consensus is both a strength and weakness of the NGA results. The weakness is that if there is a flaw in the deterministic site response approach, then all the NGA GMPEs that use Vs30 are adversely impacted. Ultimately, three data points (Table 6 and Fig. 5.3c) are insufficient for the data to significantly speak for themselves in this particular case. Consequently, one can argue for one interpretation (invariant spectral shape) or the other (spectral peaks shift to longer periods at M > 6.6), and while a Bayesian evidence analysis shows that limited available data support a spectral shift with increasing magnitude, without data from more earthquakes, an honest result is that large epistemic uncertainty remains a real issue for Vs30 > 900 m/s near-fault sites. Epistemic uncertainties can be rigorously accounted for in probabilistic ground motion analyses. However, it is necessary to develop a quantitative description of the epistemic uncertainties to accomplish this. Uncertainty in spectral shape as a function of magnitude, particular the period band of maximum acceleration response are important issues because many structures have fundamental modes of vibration at periods significantly longer than 0.2 s, the period the NGA GMPEs suggest that maximum acceleration responses will occur for M > 6.6 earthquakes at Vs30 > 900 m/s near-fault sites. We can reduce these site uncertainties and improve ground motion prediction, with the ground motion data that currently exist by collecting more quantitative information about site characteristics that more directly and robustly determine site amplification, like Vsdepth profiles. Kaklamanos and Baise (2011) showed through empirical statistical analyses that actual Vs30 measurements produced better performance than occurred at sites where Vs30 is postulated based on geology or other proxy data.  suggested that quarter-wavelength approximation of Joyner et al. (1981) would likely be a better predictor of site responses than Vs30. For a particular frequency, the quarter-wavelength approximation for amplification is given by the square root of the ratio between the seismic impedance (velocity times density) averaged over a depth corresponding to a quarter wavelength and the seismic impedance at the depth of the source. The analyses of this section suggest that the combination of Vs30 and its deterministic implementation in NGA is not the best approach. Thompson et al. (2011) show that the quarter-wavelength approximation more accurately estimates amplification than amplification estimated using Vs30. Given the rapid growth in low-cost, verified passive measurement methods to quickly estimate robust Vs-depth to 50-100 m or more (Stephenson et al., 2005;Boore and Asten, 2008;Miller et al., 2010;O'Connell and Turner, 2011), it would greatly improve the prospects for substantial improvements in future GMPEs to acquire Vs-depth data for as much of the empirical ground motion database as possible to improve resolution of site amplification.
These results illustrate how difficult it is to formulate a GMPE formulation and regression strategy a-priori, for a "single" parameter like Vs30. This analysis does not show that the NGA GMPEs are incorrect. Instead, it demonstrates some of the trade-offs, dependencies, and uncertainties that occur in the NGA GMPEs between Vs30 and spectral shape. This near-fault high Vs30 example illustrates that it is important to conduct independent analyses to determine which GMPEs are best suited for a particular application and to use multiple GMPEs, preferably with some measure of independence in their development to account for realistic epistemic GMPE uncertainties.

Near-fault application of NGA GMPEs and site-specific 3D ground motion simulations: Source and site within the basin
In tectonically active regions near plate boundaries active faults are often located within or along the margins of sedimentary basins. Basins are defined by spatially persistent strong lateral and vertical velocity contrasts that trap seismic waves within the basin. Trapped seismic waves interact to amplify ground shaking and sometimes substantially increase the duration of strong shaking. Basin amplification effect is the result the combination of lateral and vertical variations in velocity that make the basin problem truly three-dimensional in nature and difficult to quantify empirically with currently available strong motion data. The basin problem is particularly challenging for estimating amplifications for periods longer than 1 s and sedimentary basin thicknesses exceeding about 3 km . Unfortunately, some of the largest urban populations in the world are located within basins containing active faults, including many parts of Japan (Kawase, 1996;Pitarka et al., 1998, NIED, 2011, the Los Angeles and other basins in southern California (Day et al., 2008), and Seattle, Washington, (Frankel et al., 2009). Consequently, estimation of long-period ground motions in sedimentary basins associated with near-fault faulting is an important practical need.  used empirical and synthetic analyses to consider the effects of two types of basin situations. They denoted sites located in a basin overlying the source as having coincident source and site basin locations (CBL) and differentiated them from distinct source and site basin locations (DBL). They used pre-NGA GMPEs for "stiff-soil/rock", but modified to account for Vs30 using  to regress for additional basin amplification factors as a function a scalar measure of basin depth, Z1.5, the depth to a shear-wave velocity of 1.5 km/s. Using ground motion data from southern and northern California basins,  found strong empirical evidence that ground-motion amplification in coincident source and site basin locations (CBL) is significantly depthdependent at medium to long periods (T > 0.3 s). In contrast, They found that when the seismic source lies outside the basin margin (DBL), there is a much lower to negligible empirical evidence for depth-dependent basin amplification.
In support of NGA GMPE development Day et al. (2008) proposed a model for the effect of sedimentary basin depth on long-period response spectra. The model was based on the analysis of 3D numerical simulations (finite element and finite difference) of long-period 2-10 s ground motions for a suite of sixty scenario earthquakes (M 6.3 to M 7.1) within the Los Angeles basin region. Day et al. (2005) used a deterministic 3D velocity model for southern California (Magistrale et al., 2000) to calculate the wave responses on a grid and determine the amplification of basin sites as a function of Z1.5 in the 3D model. Being a purely synthetic model, but primarily concerned with ratios (amplification), it is relatively unimportant to consider the correlated-random effects on wave amplitude (Table 3) and phase (Table 4) to calculate to first-order amplification effects for shallow (< 2 km) and/or relatively fast basins. In shallow and/or fast basins the additional stochastic basin path length difference between the shallow basin and bedrock paths is less than a couple wavelengths at periods > 1 s, so the effects of differential correlated-random path lengths on S-wave amplification are negligible (O'Connell, 1999a). For typical southern California lower-velocity basins deeper than 3 km both the 3D viscoelastic finite-difference simulations of O'Connell (1999a) and phase-screen calculations of Hartzell et al. (2005) show correlated-random velocity variations will significantly reduce estimated basin amplification relative to deterministic 3D models. The primary purpose of O'Connell's (1999a) investigations was to determine the likely amplification of higher-velocity rock sites where few empirical data exist (see Section 5.2), relative to the abundant ground motion recordings obtained from stiff soil sites. O'Connell (1999a) showed that basin amplification in > 3 km deep basins is reduced relative to rock as the standard deviation of correlated-random velocity variations increases because the mean-freepath scattering in the basins significantly increases relative to rock at periods of 1-4 s. Consequently, because Day et al. (2008) use a deterministic 3D velocity model, we expect that their estimated basin amplifications will generally correspond to upper bounds of possible mean amplifications for southern California basins deeper than 3 km, but provide accurate first-order estimates of basin amplification for shallower (< 2 km) basins. Several NGA GMPEs worked to empirically evaluate and incorporate "basin effects" in some way, but it is important to note that none of the empirical NGA GMPEs explicitly consider 3D basin effects by separately considering data in coincident source and site basin locations (CBL) from other data as  showed is necessary to empirically estimate 3D basin effects for coincident source and site basin locations. NGA GMPEs lack sufficient parameterization to make this distinction, thus lumping all sites, CBL, DBL, and sites not located in basins into common Z1.0 or Z2.5 velocity-depth bins. All these sites, no matter what their actual location relative to basins and sources are apportioned some "basin-effect" through their Vs30 site velocity, Z1.0, and Z2.5 "basin-depth" terms (Day et al., 2008). It is important to understand that Z1.0 and Z2.5 are not empirically "basin-depth" terms, but "velocity-depth" terms. We use "velocity-depth" to refer to Z1.0 and Z2.5 instead of "basin-depth" because empirically, the NGA empirical GMPEs do not make the necessary distinctions in their GMPE formulations for these terms to actually apply to the problem of estimating 3D CBL amplification effects, the only basin case where a statistically significant empirical basin signal has been detected . Campbell and Bozorgnia (2008) found empirical support for significant "velocity-depth" Z2.5 term after application of their Vs30 term, but only for sites where Z2.5 < 3 km, which roughly correspond to Z1.5 < 1.5. For Z2.5 > 3 km, Campbell and Bozorgnia (2008) used the parametric 3D synthetic basin-depth model from Day et al. (2008). Day et al. (2008) note that correlation between Vs30 and "basin" depth is sufficiently strong to complicate the identification of a basin effect in the residuals after having fit a regression model to Vs30. Chiou and Youngs (2008) found that to implement a velocity-depth term using Z1.5 would require removing the Vs30 site term from their GMPE because of the Z1.5-Vs30 correlation. Instead, Chiou and Youngs (2008) retained Vs30 at all periods and included a "velocitydepth" Z1.0 term to empirically capture the portion of velocity-depth amplification not fully accounted for by the correlation between Vs30 and Z1.0.  used a similar Vs30 and Z1.0 parameterization approach for their GMPE. Since none of the NGA implementations of Z1.0 make the distinction whether a site is actually contained in a CBL or is not even in a basin, it is useful to evaluate the predictions of the four NGA GMPEs that implement Vs30, including the three NGA GMPEs that incorporate Z1.0 and Z2.5 velocitydepth terms, for four CBL sites along a portion of the North Anatolia Fault, where the fault is embedded below a series of connected 3D basins (Fig. 5.4). The North Anatolia Fault is a major right-lateral strike slip fault forming the tectonic boundary between the Eurasian plate and Anatolian plate. Over much of its length in eastern and central Turkey, the NAF is a relatively simple, localized structure (Hubert-Ferrari et al. 2002). In the Marmara region, however, the North Anatolia Fault splays into several strands and forms a wide zone of deformation between the Anatolia block and Eurasian Plate. Three hypocenters positions are used to evaluate forward, bilateral, and reverse rupture directivity (Fig. 5.4). These simulated ground motions are compared to NGA GMPE response spectra predictions from the four NGA GMPEs with Vs30, including the three with velocity-depth terms Campbell and Bozorgnia, 2008;Chiou and Youngs, 2008), and Boore and Atkinson (2008) for periods of 1 s. The NGA results are modified to account for rupture directivity using Abrahamson (2001), Spudich and Chiou (2008), and Rowshandel (2010) to isolate residual 3D basin and directivity effects relative to the NGA-based empirical predictions. A 3D velocity model encompassing the eastern Marmara Sea and Izmit Bay regions was constructed to span a region including the fault segments of interest, the ground motion estimation sites, and local earthquakes and recording stations (Fig. 5.5). Synthetic waveform modeling of local earthquake ground motions was used to iteratively improve and update the 3D model. The initial 3D velocity model was constructed using published 1-D velocity model data (Bécel et al., 2009;Bayrakci, 2009), tomographically assessed top of basement contours (Bayrakci, 2009), seismic reflection profiles (Carton et al., 2007;Kurt and Yucesoy, 2009), Boguer gravity profiles (Ates et al., 2003), geologic mapping (Okyar et al. 2008) and fault mapping (Armijo et al., 2002). Additional understanding of the basin-basement contact was gained by assessment of seismic reflection data collected by the SEISMARMARA cruise and made available at http://www.ipgp.fr/~singh/DATA-SEISMARMARA/.
The empirical wavespeed and density relations from Brocher (2005) were used to construct 3D shear-wave and density models based on the initial 3D acoustic-wave model. Shear-wave velocities were clipped so that they were not less than 600 m/s to ensure that simulated ground motions would be accurate for periods > 0.7 s for the 3D variable grid spacing used in the finite-difference calculations. This initial 3D velocity model was used to generate synthetic seismograms to compare with recordings of local M 3.2-4.3 earthquakes recorded on the margins of the Sea of Marmara, Izmit Bay, and inland locations north of Izmit Bay to assess the ground motion predictive performance of the initial 3D model. Several iterations of forward modeling were used to modify the 3D velocity model to obtain models that produce synthetic ground motions more consistent with locally-recorded local earthquake ground motions. The resulting shear-wave surface velocities mimic the pattern of acoustic-wave velocities that are consistent to first-order with the 3D acoustic-wave tomography results for the eastern Marmara Sea from Bayrakci (2009). Following O'Connell (1999a and Hartzell et al. (2005) the final 3D model incorporates 5% standard deviation correlated-random velocity variations to produce more realistic peak ground motion amplitudes than a purely deterministic model. Since there are three distinct geologic volumes in the 3D model, three independent correlatedrandomizations were used, one for the basin materials with a correlation length of 2.5 km, and one each for the basement north and south of the NAF that both used a correlation length of 5 km. Similar to Hartzell et al. (2010) we use a von Karman randomization with a Hurst coefficient close to zero and 5% standard deviation. Velocity variations are clipped to that shear-wave velocities are never smaller than 600 m/s to ensure a consistent dispersion limit for all calculations and randomized acoustic velocities were never larger than the maximum deterministic acoustic velocity to keep the same time step for all simulations. Realistic ground motion simulations require accounting for first-order anelastic attenuation, even at long periods (Olsen et al., 2003). The fourth-order finite-difference code employs the efficient and accurate viscoelastic formulation of Liu and Archuleta (2006) that accurately accounts for Q. The Vs-Qs relationship of Olsten et al. (2003) developed for the Los Angeles basin in southern California was modified slightly to ensure that Qs was bounded at the lowest shear-wave velocity by Qs=20 and at the highest shear-wave velocity by Qs=200. The simple relationship Qp=1.5*Qs was used to specify acoustic wave attenuation. These Qs-Vs and Qs-Qp relationships produced synthetic responses consistent with observed ground shaking durations and overall amplitude-frequency responses observed in local recordings of moderate-sized local earthquakes. The Karamürsel basin in Izmit Bay close to the ground motion sites is less than 2 km deep ( Fig. 5.5a). The western Cinarcik basin is deep and broad (Fig. 5.5b). The basins in Izmit Bay are separated from the Cinarcik basin and other basins in the Marmara Sea by a relatively shallow bedrock sill (Fig. 5.5c) that is well constrained by seismic reflection and refraction data from Bayrakci (2009 Liu and Archuleta (2002; was used to calculate ground motion responses from the kinematic finite fault rupture simulations. The kinematic rupture model mimics the spontaneous dynamic rupture behavior of a self-similar stress distribution model of Andrews and Boatwright (1998). The kinematic rupture model is also similar to the rupture model of Herrero and Benard (1994). Self-similar displacements are generated over the fault with rise times that are inversely proportional to effective stress. Peak rupture slip velocities evolve from ratios of 1:1 relative to the sliding (or healing peak) slip velocity at the hypocenter to a maximum ratio of 4:1. This form of slip velocity evolution is consistent with the dynamic rupture results of Andrews and Boatwright (1998) that show a subdued Kostrov-like growth of peak slip velocities as rupture grows over a fault. The kinematic model used here produces slip models with 1/k 2 (k is wavenumber) distributions consistent with estimates of earthquake slip distributions (Somerville et al., 1999) and  2 ( is angular frequency) displacement spectra in the far-field. Oglesby and Day (2002) and Schmedes et al. (2010) used numerical simulations of dynamic fault rupture to show that rupture velocity, rise time, and slip are correlated with fault strength and stress drop, as well as each other. The kinematic rupture model used here enforces correlations between these parameters by using a common fractal seed to specify relationships between all these fault rupture parameters. Oglesby and Day (2002), Guatteri et al. (2003), and Schmedes et al. (2010) used dynamic rupture simulations to demonstrate that rupture parameter correlation, as implemented in the stochastic kinematic rupture model outlined here, is necessary to produce realistic source parameters for ground motion estimation. The fault slip variability incorporates the natural log standard deviation of strike-slip displacement observed by Petersen et al. (2011) in their analyses of global measurements of strike-slip fault displacements. Consequently, although mean displacements are on the order of 1.5 m for the M 7.1 three-segment scenario earthquake, asperities within the overall rupture have displacements of up to 3-4 m. The Liu et al. (2006) slip velocity function is used with the specified fault slips and rise times to calculate slipvelocity time functions at each grid point. Three hypocenters were used to simulate forward, reverse, and bilateral ruptures relative to Izmit Bay sites (Fig. 5.4). To find an appropriate "median" randomization of the 3D velocity model, ten correlated-random 3D velocity models were created and then a single, threesegment randomized kinematic rupture model was used to simulate ten sets of ground motions. The randomized 3D model that most consistently produced nearly median motions across the five sites over the 1-10 s period band was used to calculate all the ground motion simulations for all two-segment and three-segment rupture scenarios. Ten kinematic randomizations were used for each case resulting in 60 rupture scenario ground motion simulations.
The simulated ground motions were post-processed to calculate acceleration response spectra for 5% damping. The geometric mean of Boore et al. (2006) (GMRotI50) was calculated from the two horizontal components to obtain GMRotI50 response spectra (SA). Response spectral results are interpreted for periods longer than 1 s, consistent with the fourth-order finite-difference accuracy for the variable grid spacing, minimum shear-wave velocity of 600 m/s, and broad period influence of oscillator response (Day et al., 2008). Fault-normal peak velocities decrease from sites 1-3 close to the fault and near the deeper portion of the basin (Fig. 5.5a) toward the shallow basin (site 4 in Figs. 5.4a and Fig. 5.6), and bedrock outside the basin (site 5 in Figs. 5.4a and Fig. 5.6). The four NGA ground motion prediction equations (GMPE) that implement Vs30 were used to calculate ground motion estimates at all four sites using the Z1.0 and Z2.5 below each site in the 3D synthetic velocity model, Vs30=600 m/s, the directivity corrections of Abrahamson (2001), Spudich and Chiou (2008), and Rowshandel (2010) equally weighted, and the three rupture hypocenters (forward, bilateral, and reverse directivity in Fig. 5.3) equally weighted. Site 4 was located away from basin-edge effects and in the shallow portion of the basin with the same Vs30=600 m/s, as sites 1-3 and consistent with a relatively linear site response making direct comparison of linear 3D simulated motions with empirical GMPE feasible. Site 4 horizontal spectra where estimated as the log-mean average of the set of six earthquake rupture scenarios (two-segment and three-segment rupture, and forward, bilateral, and reverse directivity) used in the 3D ground motion simulations. To obtain robust estimates of mean synthetic spectra, we omitted the two largest and smallest amplitudes at each period to estimate log-mean spectra for comparison. The Site 4 horizontal responses are comparable in amplitude to NGA predicted response spectra for periods > 1 s (Fig. 5.7a). The reduced synthetic responses between 1-2 s in Fig. 5.7 are an artifact of finite-difference grid dispersion similar to that noted by Day et al. (2008). The site 4 3D simulated responses are generally slightly less than the empirical GMPE median estimates over the 1-8 second period range, except for a small amplification at 3 seconds of < 10% (Fig. 5.7b). This confirmed that the 3D ground motion simulations and empirical NGA GMPE predict comparable spectral hazard at site 4 and establish site 4 as an appropriate reference point to compare to responses at sites 1-3 closer to the fault and within deeper portions of the basin. We use the empirical NGA GMPEs to estimate the amplitude effects of differential sourcesite distances on amplitudes relative to site 4 and changes in Z1.0 and Z2.5 between sites 1-3 and reference site 4. The three empirical directivity relations are used with equal weight to remove the differential directivity effects for two separate sets of rupture cases designed to determine if 3D basin amplification is dependent on rupture directivity. We consider in the first case, the two rupture scenarios away from the sites, to determine 3D basin amplification in the absences of forward rupture directivity. In the second case, we average all six rupture scenarios, four of which have strong forward rupture directivity, to see if any of the sites shows significantly 3D basin amplification in the case of solely reverse rupture direction ground motions. Empirical NGA distance, directivity, and Z1.0-Z2.5 sites 1-3 amplifications relative to reference site 4 are the lowest curves at the bottom of Fig. 5.8 and represent the sum total of the effects of all NGA GMPE terms related to differential distance, directivity, and Z1.0 and Z2.5 velocity-depth. Although sites 1-3 are much closer to the fault than site 4, the relative Fig. 5.8. Mean and reverse-rupture only residual 3D basin amplifications for sites 1-3 relative to reference site 4 with NGA differential site amplitude correction functions. changes in amplitudes are much smaller than the proportional differences in site-source distances as a result of saturation, the condition enforced in NGA that ground motion amplitudes cease to increase as distance to the fault approaches zero. The directivity amplitude reduction from Rowshandel (2010) for reverse rupture accounts for the dip in longer-period NGA differential responses at periods of about 5 s in Fig. 5.8. The most striking aspect of the NGA transfer functions is that although three of the four GMPE include Z1.0 or Z2.5 "basin-depth" terms, there is no hint of an empirical resonant 3D basin response, just slight steady increases of "amplification" with increasing period. The non 3Dbasin-like NGA differential amplification results are not surprising because the NGA basindepth formulation pools ground motion observations from all scales of basins and nonbasins in each Z1.0 and Z2.5 bin. Consequently, the NGA Vs30 and velocity-depth Z1.0 and Z2.5 basin terms do not capture any of the strongly period-dependent amplification associated with the site-specific basin of < 2 km total depth near the sites. The residual site-specific synthetic 3D amplifications at sites 1-3 relative to reference site 4 are essentially independent of rupture direction ( Fig. 5.8). Site 1 closest to the fault shows the largest amplification for case 2 with 2/3 forward rupture directivity, but the difference at site 1 between 2/3 forward rupture directivity basin amplification and reverse rupture basin amplification is < 10%. For sites 2 and 3 located slightly further from the fault, differences in case one and case two directivity 3D basin amplifications deviate < 4% from their mean peak amplifications. The remarkable result is that even in this case of a strike-slip fault embedded below the center of a basin and rupturing within basins continually along the entire rupture length, to first order 3D basin amplification is independent of rupture directivity/rupture direction. These 3D synthetic calculations shows that the three empirical directivity corrections applied with NGA GMPE effectively accounted for first-order directivity in this rather severe case of strike-slip fault rupture within a basin. The Izmit Bay basins are quite similar in width, depth, and velocity characteristics to the San Fernando Basin, one of the basins included in the Day et al. (2008) 3D synthetic calculations to represent basin amplification in younger shallower basins. Thus, it is interesting to compare the Day et al. (2008) synthetic amplification predictions calculated across a spectrum of shallow and deeper basins using a deterministic 3D velocity model with these simulations using a site-specific weakly-randomized 3D velocity model. We calculate the 3D simulation and Day et al. (2008) response ratios of sites 1-3 to site 4 using the Z1.5 values from the 3D simulation model in the Day et al. (2008) Z1.5-amplification relationships ( Fig. 5.9). Both 3D synthetic approaches predict comparable peak amplifications at comparable periods ( Fig. 5.9), with the site-specific 3D model predicting a more rapid decrease with in increasing period that reflects the details of the site-specific 3D model; Day et al. (2008) have a wider period range of larger amplification because they pooled basin amplifications from a wider range of basin configurations than representative of the site-specific 3D velocity structure. When soils are significantly less linear than clays with a plasticity index of 20, the fully nonlinear shear P-SV 2D investigations of O'Connell et al. (2010) suggest that combining the outputs of linear 3D simulations that omit the very-low-velocity basin with 1D nonlinear analyses to account for the very-low-velocity basins will produce comparable amplifications within the basin to full nonlinear 2D or 3D analyses. Linear 1D P-SV vertical analyses in the central portions of basins will typically provide appropriate vertical amplifications throughout most of the basin. Thus, it appears that it may be feasible in most of these cases to omit the shallow soft low-velocity regions from the top of basins from 3D linear or nonlinear analyses and use the outputs from linear 3D analyses with simplified 1D nonlinear SH and P-SV nonlinear amplification calculations to estimate realistic horizontal and vertical peak velocities and accelerations in the upper low-velocity soft soils. These results illustrate that at present, the NGA GMPE do not effectively estimate sitespecific 3D basin amplification for the most extreme case of a strike slip source and sites located within a closed basin. In such situations it is necessary to use site-specific 3D basin amplification calculations or compiled synthetic 3D generic basin amplification relations like Day et al. (2008) to estimate realistic site-specific 3D basin amplification effects. However, the NGA GMPEs and associated empirical directivity relations are shown to effectively account for geometric spreading and directivity in the demanding application of source and site located within a closed basin and provide a robust means to extract residual 3D basin amplification relative to NGA GMPE predictions. This approach requires a suitable reference site in shallow portions of the basin that are not strongly influenced by basin effects or a site outside the basin. In future GMPE development, the basin analyses of , Day et al. (2008), and this analysis suggest separate consideration and analysis of data that is within closed basins with faults beneath or adjacent to the basin is warranted to evaluate empirical evidence for systematic basin responses. Such analyses need to be done separately for ground motion observations outside of this specific basin configuration to discern the relative effects of velocity-depth versus basin-depth on parameters like Vs30, Z1.0, and Z1.5. We suggest that it is more appropriate and prudent to refer to Z1.0 and Z1.5 as velocity-depth terms, not basin-terms, since they will fail to account for significant systematic period-dependent 3D basin amplification in the cases of sources and sites located within low-velocity basins.

Conclusion and recommendations
Geologic seismic source characterization is the fundamental first step in strong ground motion estimation. Many of the largest peak ground motion amplitudes observed over the past 30 years have occurred in regions where the source faults were either unknown or major source characteristics were not recognized prior to the occurrence of earthquakes on them. The continued development of geologic tools to discern and quantify fundamental characteristics of active faulting remains a key strong ground motion estimation research need. As Jennings (1983) noted, by the early 1980s efforts to develop empirical ground motion prediction equations were hampered not only by the insufficient recordings of ground motions to constrain the relationships between magnitude, distance, and site conditions, but insufficient physical understanding of how to effectively formulate the problem. Strong ground motion estimation requires both strong motion observations and understanding of the physics of earthquake occurrence, earthquake rupture, seismic radiation, and linear and nonlinear wave propagation. In sections 2-4 we provided an overview of the physics of strong ground motions and forensic tools to understand their genesis. The physics are complex, requiring understanding of processes operating on scales of mm to thousands of km, most of the physical system is inaccessible, and the strong motion observations are sparse. As O'Connell (1999a) and Hartzell et al. (2005) showed, surface ground motion observations alone are insufficient to constrain linear and nonlinear amplification and seismic source properties. The observational requirements to understand the earthquake system and how ground motions are generated are immense, and require concurrent recording of ground motions at the surface and at depth. These observations have only recently been undertaken at a comprehensive large scale. In Japan, the National Research Institute for Earth Science and Disaster Prevention (NIED) operates K-NET (Kyoshin Network) with 660 strong motion stations. Each station records triaxial accelerations both at the surface and at sufficient depth in rock to understand the physics of earthquake fault rupture and to directly observe linear and nonlinear seismic wave propagation in the shallow crust. These borehole-surface data have provided fundamental new constraints on peak ground motions (Aoi et al., 2008), direct observation of nonlinear wave propagation, and new constraints on ground motion variability (Rodriguez-Marek et al., 2011). It will be necessary to expand the deployment of K-NET scale networks to other tectonically active regions like the western United States, to make real long-term progress understanding and significantly improving our ability to predict strong ground shaking. The synergy between earthquake physics research and strong ground motion estimation is based on ground motion observations and geologic knowledge. The need for new recordings of strong ground motions in new locations is clear, but there is immensely valuable information yet to be extracted from existing strong ground motion data. One of single biggest impediments to understanding strong ground motions is the lack of site velocity measurements for most of the current strong ground motion database Kaklamanos and Baise, 2011). The last 10 years has seen an explosion in the development and successful application of rapid, inexpensive, and non-invasive methods to measure site shear-wave velocities over depths of 50-1000 m that can provide site amplification estimates accurate to on the order of 10-20% (Stephenson et al., 2005;Boore and Asten, 2008). Using the large borehole-surface station network in Japan, Rodriguez-Marek et al. (2011) showed that the difference in the single-station standard deviation of surface and borehole data is consistently lower than the difference in ergodic standard deviations of surface and borehole data. This implies that the large difference in ergodic standard deviations can be attributed to a poor parameterization of site response. Vs30 does not constrain frequency-dependent site amplification because literally, an infinite number of different site velocity-depth profiles can have the same Vs30. Even given geologic constraints on near-surface material variability, the scope of distinct velocity profiles and amplification characteristics that share a common Vs30 is vast. Ironically, the implementation of Vs30 in four of the NGA GMPE produced significant uncertainties in spectral shape as a function of magnitude as illustrated in Section 5.1. Vs30 also trades off with other velocity-depth factors (Section 5.2). We propose that one of the most valuable new strong ground motion datasets that can be obtained now is measurement of site shearwave velocity profiles at the sites of existing strong ground motion recordings. These measurements would provide a sound quantitative basis to constrain frequency-dependent linear-site amplification prior to regression and reduce uncertainties in ground motion estimations, particularly spectral shape as a function of site conditions. As Rodriguez-Marek et al. (2011) note, reduction of exaggerated ground motion variability results in more realistic ground motion estimates across widely differing sites in probabilistic analyses. The analyses of  and section 5.2 suggest that accounting for positions of ground motion recordings and earthquake inside or outside of closed basins may provide a path forward to improve the ability of future empirical GMPE to accurately estimate responses within basins.

Acknowledgments
This paper is dedicated to the memory of William Joyner, who generously participated in discussions of directivity, wave propagation, site response, and nonlinear soil response, and who encouraged us to pursue many of the investigations presented here. David Boore kindly read the initial draft and provided suggestions that improved it.