Daniel R H O'connell

Jon P Ake

Fabian Bonilla

Pengcheng Liu

Roland Laforge

Dean Ostenaa

Strong Ground Motion Estimation

Introduction

At the time of its founding, only a few months after the great 1906 M 7.7 San Francisco Earthquake, the Seismological Society of America noted in their timeless statement of purpose "that earthquakes are dangerous chiefly because we do not take adequate precautions against their effects, whereas it is possible to insure ourselves against damage by proper studies of their geographic distribution, historical sequence, activities, and effects on buildings." Seismic source characterization, strong ground motion recordings of past earthquakes, and physical understanding of the radiation and propagation of seismic waves from earthquakes provide the basis to estimate strong ground motions to support engineering analyses and design to reduce risks to life, property, and economic health associated with earthquakes. When a building is subjected to ground shaking from an earthquake, elastic waves travel through the structure and the building begins to vibrate at various frequencies characteristic of the stiffness and shape of the building. Earthquakes generate ground motions over a wide range of frequencies, from static displacements to tens of cycles per second [Hertz (Hz)]. Most structures have resonant vibration frequencies in the 0.1 Hz to 10 Hz range. A structure is most sensitive to ground motions with frequencies near its natural resonant frequency. Damage to a building thus depends on its properties and the character of the earthquake ground motions, such as peak acceleration and velocity, duration, frequency content, kinetic energy, phasing, and spatial coherence. Strong ground motion estimation must provide estimates of all these ground motion parameters as well as realistic ground motion time histories needed for nonlinear dynamic analysis of structures to engineer earthquake-resistant buildings and critical structures, such as dams, bridges, and lifelines. Strong ground motion estimation is a relatively new science. Virtually every M > 6 earthquake in the past 35 years that provided new strong ground motion recordings produced a paradigm shift in strong motion seismology. The 1979 M 6.9 Imperial Valley, California, earthquake showed that rupture velocities could exceed shear-wave velocities over a significant portion of a fault, and produced a peak vertical acceleration > 1.5 g [START_REF] Spudich | Direct observation of rupture propagation during the 1979 Imperial Valley earthquake using a short baseline accelerometer array[END_REF]Archuleta;[START_REF] Archuleta | A faulting model for the 1979 Imperial Valley earthquake[END_REF]. The 1983 M 6.5 Coalinga, California, earthquake revealed a new class of seismic sources, blind thrust faults [START_REF] Stein | Seismicity and geometry of a 110-km-long blind thrust fault: 2. Synthesis of the 1982-1985 California earthquake sequence[END_REF]. The 1985 M 6.9 Nahanni earthquake produced horizontal accelerations of 1.2 g and a peak vertical acceleration > 2 g (Weichert et al., 1986). The 1989 M 7.0 Loma Prieta, California, earthquake occurred on an unidentified steeply-dipping fault adjacent to the San Andreas fault, with reverse-slip on half of the fault [START_REF] Hanks | The 1989 Loma Prieta, California, earthquake and its effects: Introduction to the Special Issue[END_REF], and produced significant damage > 100 km away related to critical reflections of shear-waves off the Moho [START_REF] Somerville | The influence of critical Moho reflections on strong ground motions recorded in San Francisco and Oakland during the 1989 Loma Prieta earthquake[END_REF][START_REF] Catchings | Reflected seismic waves and their effect on strong shaking during the 1989 Loma Prieta, California, earthquake[END_REF]. The 1992 M 7.0 Petrolia, California, earthquake produced peak horizontal accelerations > 1.4 g [START_REF] Oglesby | A faulting model for the 1992 Petrolia earthquake: Can extreme ground acceleration be a source effect?[END_REF]. The 1992 M 7.4 Landers, California, earthquake demonstrated that multisegment fault rupture could occur on fault segments with substantially different orientations that are separated by several km [START_REF] Li | Fine structure of the Landers fault zone; segmentation and the rupture process[END_REF]. The 1994 M 6.7 Northridge, California, earthquake produced a then world-record peak horizontal velocity (> 1.8 m/s) associated with rupture directivity (O'Connell, 1999a), widespread nonlinear soil responses [START_REF] Field | Nonlinear ground-motion amplification by sediments during the 1994 Northridge earthquake[END_REF][START_REF] Cultera | Nonlinear soil response in the vicinity of the Van Norman Complex following the 1994 Northridge, California, earthquake[END_REF], and resulted in substantial revision of existing ground motion-attenuation relationships [START_REF] Abrahamson | Overview[END_REF]. The 1995 M 6.9 Hyogoken Nanbu (Kobe) earthquake revealed that basin-edge generated waves can strongly amplify strong ground motions [START_REF] Kawase | The cause of the damage belt in Kobe: "The basin-edge effect," constructive interference of the direct S-wave with the basin-induced diffracted/Rayleigh waves[END_REF][START_REF] Pitarka | Three-dimensional simulation of the near-fault ground motions for the 1995 Hyogo-Nanbu (Kobe), Japan, earthquake[END_REF] and provided ground motion recordings demonstrating time-dependent nonlinear soil responses that amplified and extended the durations of strong ground motions [START_REF] Archuleta | Nonlinearity in observed and computed accelerograms[END_REF]. The 1999 M > 7.5 Izmit, Turkey, earthquakes produced asymmetric rupture velocities, including rupture velocities ~40% faster than shear-wave velocities, which may be associated with a strong velocity contrast across the faults [START_REF] Bouchon | How Fast is Rupture during an Earthquake? New Insights from the 1999 Turkey Earthquakes[END_REF]. The 1999 M 7.6 Chi-Chi, Taiwan, earthquake produced a world-record peak velocity > 3 m/s with unusually low peak accelerations [START_REF] Shin | A preliminary report of the 1999 Chi-Chi (Taiwan) earthquake[END_REF]. The 2001 M 7.7 Bhuj India demonstrated that M > 7.5 blind thrust earthquakes can occur in intraplate regions. The M 6.9 2008 Iwate-Miyagi, Japan, earthquake produced a current world-record peak vector acceleration > 4 g, with a vertical acceleration > 3.8 g (Aoi et al., 2008). The 2011 M 9.1 Tohoku, Japan, earthquake had a world-record peak slip on the order of 60 m (Shao et al., 2011) and produced a world-record peak horizontal acceleration of 2.7 g at > 60 km from the fault [START_REF] Nied | Off the Pacific Coast of Tohoku Earthquake, Strong Ground Motion[END_REF]. 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. [START_REF] Lettis | Empirical observations regarding reverse earthquakes, blind thrust faults, and Quaternary deformation: Are blind thrust faults truly blind[END_REF] 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, 1999 M 7.6 Chi-Chi Taiwan, 2001 M 7.7 Bhuj, India, 2010 M 7.0 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 [START_REF] Nied | Off the Pacific Coast of Tohoku Earthquake, Strong Ground Motion[END_REF] 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 [START_REF] Spudich | Techniques for earthquake ground-motion calculation with applications to source parameterization to finite faults[END_REF]

        u t s t g t k i j k i j ij nm    (1)
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 [START_REF] Hutchings | Empirical Green's functions from small earthquakes -A waveform study of locally recorded aftershocks of the San Fernando earthquakes[END_REF] and [START_REF] Hutchings | Kinematic earthquake models and synthesized ground motions using empirical Green's functions[END_REF]

using,             U S e G e k i j i kij ij nm i ij kij              
(2)

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 [START_REF] Spudich | Use of ray theory to calculate high-frequency radiation from earthquake sources having spatially variable rupture velocity and stress drop[END_REF]1987) 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 [START_REF] Spudich | Use of ray theory to calculate high-frequency radiation from earthquake sources having spatially variable rupture velocity and stress drop[END_REF]Zeng et al., 1991;Smedes and Archuleta, 2008)

Isochrones analysis of rupture directivity

        v t f t dl y t x sGc x ,  ,     (3)         2 , , y t d d d a t f t dl dq dq dq                               2 2 x s G c x c G c s s G sGc ( 4
)
where c is isochrone velocity, s is slip velocity (either rupture or healing), G is a ray theory Green function, x are position vectors, y(t,x) are isochrones,  is the curvature of the isochrone, dl denotes isochrone line integral integration increment, and dq denotes a spatial derivative. Since isochrones are central to understanding ground motions, we provide explicit expressions for rupture and healing isochrones to illustrate how source and propagation factors can combine to affect ground motions. The arrival times of rupture at a specific receiver are

      T t t r r x x      , (5)
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

      T T R h r x x   ,  (6)
where R are the rise times (the durations of slip) at all fault positions. [START_REF] Archuleta | A faulting model for the 1979 Imperial Valley earthquake[END_REF] 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 [START_REF] Spudich | Techniques for earthquake ground-motion calculation with applications to source parameterization to finite faults[END_REF]. This is not the case with rupture velocity or rise time; these quantities influence ground motions at all sites. However, as [START_REF] Anderson | Comparison of strong ground motion from several dislocation models[END_REF] showed, it takes a 300% change in rise time to compensate for a 17% change in rupture time. [START_REF] Spudich | Dense seismograph array observations of earthquake rupture dynamics[END_REF] 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) is caused by variations of isochrone velocity c, where

c T r    s 1 (7)
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 [START_REF] Farra | Fast near source evaluation of strong ground motion for complex source models[END_REF]. [START_REF] Spudich | Use of ray theory to calculate high-frequency radiation from earthquake sources having spatially variable rupture velocity and stress drop[END_REF] 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.

3.1.1

The fundamental difference between strike-slip and dip-slip directivity [START_REF] Boore | The effect of directivity on the stress parameter determined from ground motion observations[END_REF] and [START_REF] Joyner | Directivity for nonuniform ruptures[END_REF] 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, [START_REF] Schmedes | Near-source ground motion along strike-slip faults: Insights into magnitude saturation of PGV and PGA[END_REF] use kinematic rupture simulations and isochrone analyses to show why directivity becomes bounded during strike-slip fault along long faults. [START_REF] Schmedes | Near-source ground motion along strike-slip faults: Insights into magnitude saturation of PGV and PGA[END_REF] 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;Abrahamson and Silva, 2008;[START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF][START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF][START_REF] Chiou | An NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF]. 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 [START_REF] O'connell | Influence of dip and velocity heterogeneity on reverse-and normal-faulting rupture dynamics and near-fault ground motions[END_REF]. 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 [START_REF] Schmedes | Near-source ground motion along strike-slip faults: Insights into magnitude saturation of PGV and PGA[END_REF]. When the two-dimensional nature of finite-fault rupture is considered, rupture directivity is not as strong as suggested by this one-dimensional analysis [START_REF] Bernard | Modeling directivity of heterogeneous earthquake ruptures[END_REF], 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

    ij .
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 [START_REF] Brune | Tectonic stress and the spectra of seismic shear waves from earthquakes[END_REF]. 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 [START_REF] Aki | Strong-motion seismology[END_REF]. 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 [START_REF] Wells | New empirical relationships amoung magnitude, rupture length, rupture width, rupture area, and surface displacement[END_REF][START_REF] Somerville | Characterizing crustal earthquake slip models for the prediction of strong ground motion[END_REF]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. [START_REF] Andrews | A stochastic fault model, 2, Time-dependent case[END_REF] and [START_REF] Frankel | High-frequency spectral falloff of earthquakes, fractal dimension of complex rupture, b value, and the scaling strength on faults[END_REF] 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. [START_REF] Somerville | Characterizing crustal earthquake slip models for the prediction of strong ground motion[END_REF] 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, [START_REF] Somerville | Characterizing crustal earthquake slip models for the prediction of strong ground motion[END_REF] provide recipes for specifying fault rupture parameters such as slip, rise times, and asperity dimensions as a function of moment. [START_REF] Mai | Source scaling properties from finite-fault-rupture models[END_REF] showed that 5.3 < M < 8.1 magnitude range dip-slip earthquakes follow self-similar scaling as suggest by [START_REF] Somerville | Characterizing crustal earthquake slip models for the prediction of strong ground motion[END_REF]. 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 [START_REF] Oglesby | The three-dimensional dynamics of dipping faults[END_REF] suggest that self-similar slip scaling relations may also breakdown at larger moments for dip-slip events.

Factor

Influence Moment rate,

   S ij M 0 
Moment rate scales peak velocities and accelerations. Moment determines the average slip for a fixed fault area and known shear moduli.

Stress drop,

   S ij   Since    S ij     ,
   S ij C 
Diffraction at the crack tip introduces a frequency dependent amplitude to the radiation pattern [START_REF] Madariaga | High-frequency radiation from crack (stress-drop) models of earthquake faulting[END_REF][START_REF] Boatwright | A dynamic model for far-field acceleration[END_REF][START_REF] Fukuyama | Integral equation method for plane crack with arbitary shape in 3D elastic medium[END_REF].

Dynamics,    S ij D 
Fault rupture in heterogeneous velocity structure can produce anisotropic slip velocities relative to rupture direction [START_REF] Harris | Effects of a low-velocity zone on dynamic rupture[END_REF] and slip velocities and directivity are a function of rake and dip for dip-slip faults [START_REF] Oglesby | Earthquakes on dipping faults:The effects of broken symmetry[END_REF]2000;[START_REF] O'connell | Influence of dip and velocity heterogeneity on reverse-and normal-faulting rupture dynamics and near-fault ground motions[END_REF]. Frictional heating, fault zone fluids, and melting may also influence radiated energy (Kanamori and Brodsky, 2001;Andrews, 200X)

   S ij  ) Factor Influence Rupture velocity,     ij V r
High rupture velocities increase directivity. Rupture velocities interact with stress drops and rise times to modify the amplitude spectrum. Supershear rupture velocities can increase directivity far from the fault (Andrews, 2010). Healing velocity,

    ij V h
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,

    ij A
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 [START_REF] Madariaga | High-frequency radiation from crack (stress-drop) models of earthquake faulting[END_REF][START_REF] Boatwright | A dynamic model for far-field acceleration[END_REF][START_REF] Fukuyama | Integral equation method for plane crack with arbitary shape in 3D elastic medium[END_REF].

    ij R Since  c R  1 ,
Dynamics,     ij D
The same dynamic processes identified in Table 1 produce corresponding phase variability. (1998; 2000) 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. [START_REF] O'connell | Influence of dip and velocity heterogeneity on reverse-and normal-faulting rupture dynamics and near-fault ground motions[END_REF] 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. [START_REF] Day | Three-dimensional simulation of spontaneous rupture: The effect of nonuniform prestress[END_REF] 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;[START_REF] Spudich | Direct observation of rupture propagation during the 1979 Imperial Valley earthquake using a short baseline accelerometer array[END_REF][START_REF] Archuleta | A faulting model for the 1979 Imperial Valley earthquake[END_REF], the M 6.9 1980 Irpinia normal-faulting earthquake [START_REF] Belardinelli | Redistribution of dynamic stress during coseismic ruptures: Evidence for fault interaction and earthquake triggering[END_REF], the M 7.0 1992 Petrolia thrust-faulting earthquake [START_REF] Oglesby | A faulting model for the 1992 Petrolia earthquake: Can extreme ground acceleration be a source effect?[END_REF], the M 7.3 Landers strike-slip earthquake [START_REF] Olsen | Three-dimensional dynamic simulation of the 1992 Landers earthquake[END_REF][START_REF] Bouchon | Stress field associated with the rupture of the 1992 Landers, California, earthquake and its implications concerning the fault strenght at the onset of the earthquake[END_REF][START_REF] Hernandez | Contribution of radar interfermetry to a two-step inversion of the kinematic process of the 1992 Landers earthquake[END_REF] the M 6.7 1994 Northridge thrust-faulting earthquake [START_REF] O'connell | Possible super-shear rupture velocities during the 1994 Northridge earthquake[END_REF], and the 1999 M 7.5 Izmit and M 7.3 Duzce Turkey strike-slip earthquakes [START_REF] Bouchon | How Fast is Rupture during an Earthquake? New Insights from the 1999 Turkey Earthquakes[END_REF]. 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. [START_REF] O'connell | Possible super-shear rupture velocities during the 1994 Northridge earthquake[END_REF] 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. [START_REF] Harris | Effects of a low-velocity zone on dynamic rupture[END_REF] 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 [START_REF] Bouchon | How Fast is Rupture during an Earthquake? New Insights from the 1999 Turkey Earthquakes[END_REF]. 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. [START_REF] Somerville | Simulations of strong ground motions recorded during the Michoacan, Mexico and Valparaiso, Chile, earthquakes[END_REF][START_REF] Schneider | Ground motion model for the 1989 M 6.9 Loma Prieta earthquake including effects of source, path, and site[END_REF][START_REF] Hutchings | Kinematic earthquake models and synthesized ground motions using empirical Green's functions[END_REF][START_REF] Tumarkin | Scaling relations for composite earthquake models[END_REF][START_REF] Zeng | A composite source model for computing realistic strong ground motions[END_REF][START_REF] Beresnev | Modeling finite-fault radiation from the n spectrum[END_REF]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 [START_REF] Beresnev | Modeling finite-fault radiation from the n spectrum[END_REF], and larger uncertainties for sites subject to directivity. Rupture direction may change due to strength or stress heterogeneities on a fault. [START_REF] Beroza | Linearized inversion for fault rupture behavior: Application to the 1984 Morgan Hill, California, earthquake[END_REF] inferred that rupture was delayed and then progressed back toward the hypocenter during the M 6.2 1984 Morgan Hill earthquake. [START_REF] Oglesby | A faulting model for the 1992 Petrolia earthquake: Can extreme ground acceleration be a source effect?[END_REF] 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. [START_REF] Das | A numerical study of two-dimensional spontaneous rupture propagation[END_REF] 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. [START_REF] Day | Three-dimensional simulation of spontaneous rupture: The effect of nonuniform prestress[END_REF] 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, [START_REF] Das | Breaking of a single asperity: Rupture process and seismic radiation[END_REF] 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. [START_REF] Guatteri | Coseismic temporal changes of slip direction: The effect of absolute stress on dynamic rupture[END_REF] 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. [START_REF] Pitarka | Simulation of near-fault strong-ground motion using hybrid Green's functions[END_REF] 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. [START_REF] Spudich | Use of fault striations and dislocation models to infer tectonic shear stress during the 1995 Hyogo-ken Nanbu (Kobe) earthquake[END_REF] 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. [START_REF] Vidale | Influence of focal mechanism on peak accelerations of strong motions of the Whittier Narrows, California, earthquake and an aftershock[END_REF] 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, [START_REF] Liu | The 23:19 aftershock of the 15 October 1979 Imperial Valley earthquake: More evidence for an asperity[END_REF] 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. [START_REF] Bent | Source complexity of the October 1, 1987, Whittier Narrows earthquake[END_REF] 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 [START_REF] Vidale | Influence of focal mechanism on peak accelerations of strong motions of the Whittier Narrows, California, earthquake and an aftershock[END_REF]. 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. [START_REF] Frankel | A three-dimensional simulation of seimic waves in the Santa Clara Valley, California, from a Loma Prieta aftershock[END_REF][START_REF] Frankel | Three-dimensional simulations of ground motions in the San Bernardino Valley, California, for hypothetical earthquakes on the San Andreas fault[END_REF][START_REF] Olsen | Three-dimensional simulation of earthquakes on the Los Angeles fault system[END_REF][START_REF] Wald | The seismic response of the Los Angeles Basin, California[END_REF][START_REF] Archuleta | Direct observation of nonlinear soil response in acceleration time histories[END_REF][START_REF] Frankel | Three-dimensional simulations of ground motins in the Seattle region for earthquakes in the Seattle fault zone[END_REF][START_REF] Koketsu | Propagation of seismic ground motion in the Kanto Basin, Japan[END_REF][START_REF] Frankel | Observations of basin ground motions from a dense seismic array in San Jose, California[END_REF]. Basin-edge waves can substantially amplify strong ground motions in basins [START_REF] Liu | Array analysis of the ground velocities and accelerations from the 1971 San Fernando, California, earthquake[END_REF][START_REF] Frankel | High-frequency spectral falloff of earthquakes, fractal dimension of complex rupture, b value, and the scaling strength on faults[END_REF][START_REF] Phillips | Basin-induced Love waves observed using the strong motion array at Fuchu, Japan[END_REF][START_REF] Spudich | The seismic coda, site effects, and scattering in alluvial basins studied using aftershocks of the 1986 North Palm Springs, California, earthquakes as source arrays[END_REF][START_REF] Kawase | The cause of the damage belt in Kobe: "The basin-edge effect," constructive interference of the direct S-wave with the basin-induced diffracted/Rayleigh waves[END_REF][START_REF] Pitarka | Three-dimensional simulation of the near-fault ground motions for the 1995 Hyogo-Nanbu (Kobe), Japan, earthquake[END_REF][START_REF] Frankel | Observations of basin ground motions from a dense seismic array in San Jose, California[END_REF]. 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 [START_REF] Kawase | The cause of the damage belt in Kobe: "The basin-edge effect," constructive interference of the direct S-wave with the basin-induced diffracted/Rayleigh waves[END_REF][START_REF] Pitarka | Three-dimensional simulation of the near-fault ground motions for the 1995 Hyogo-Nanbu (Kobe), Japan, earthquake[END_REF], 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 [START_REF] Alex | Lens-effect in Santa Monica?[END_REF][START_REF] Graves | Ground motion amplification in the Santa Monica area: Effects of shallow basin structure[END_REF][START_REF] Davis | Northridge earthquake damage caused by geologic focusing of seismic waves[END_REF]. Basin-edge waves can be composed of both body and surface waves [START_REF] Spudich | The seismic coda, site effects, and scattering in alluvial basins studied using aftershocks of the 1986 North Palm Springs, California, earthquakes as source arrays[END_REF][START_REF] Meremonte | Urban seismology: Northridge aftershocks recorded by multiscale arrays of portable digital seismographs[END_REF][START_REF] Frankel | Observations of basin ground motions from a dense seismic array in San Jose, California[END_REF]) 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 [START_REF] Somerville | The influence of critical Moho reflections on strong ground motions recorded in San Francisco and Oakland during the 1989 Loma Prieta earthquake[END_REF][START_REF] Catchings | Reflected seismic waves and their effect on strong shaking during the 1989 Loma Prieta, California, earthquake[END_REF]. 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, [START_REF] Catchings | New Madrid and central California apparent Q values as determined from seismic refraction data[END_REF] 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

Factor Influence Geometric spreading,   G kij r 
Amplitudes decrease with distance at 1/r, 1/r 2 , and 1/r 4 for body waves and 1/ r for surface waves. The 1/r term has the strongest influence on high-frequency ground motions. The 1/ r term can be significant for locally generated surface waves. Large-scale velocity structure,

  G kij V D 3 
Horizontal and vertical velocity gradients and velocity discontinuities can increase or decrease amplitudes and durations. Low-velocity basins can amplify and extend ground motion durations. Abrupt changes in lateral velocity structure can induce basin-edge-waves in the lower velocity material that amplify ground motions. Near-surface resonant responses,

  G kij L  Low-
  G kij Q 
Linear hysteretic behavior that reduces amplitudes of the form

e f r Q       . High-frequency atten uation,   G kij  
Strong attenuation of high-frequencies in the shallow crust of the form

  e r f      . Scattering,   G kij S 
Scattering tends to reduce amplitudes on average, but introduces high amplitude caustics and low-amplitude shadow zones and produces nearly log-normal distributions of amplitudes (O'Connell, 1999a).

Anisotropy,

 

G kij A  Complicates shear-wave amplitudes and modifies radiation pattern amplitudes and can introduce frequency-dependent amplification based on direction of polarization.

Topography,

 

G kij T  Can produce amplification near topographic highs and introduces an additional sources of scattering.

Table 3. Seismic Wave Propagation Amplitude Factors (

  G kij  )
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. [START_REF] Borcherdt | Progress on ground motion predictions for the San Francisco Bay region, California, in Progress on Seismic Zonation in the San Francisco Bay Region[END_REF][START_REF] Joyner | The effect of Quaternary alluvium on strong ground motion in the Coyote Lake, California earthquake of 1979[END_REF][START_REF] Seed | The Mexico earthquake of September 19, 1985 -relationships between soil conditions and earthquake ground motions[END_REF]. [START_REF] Williams | Surface seismic measurements of near-surface P-and S-wave seismic velocities at earthquake recording stations[END_REF] showed that significant resonances can occur for impedance boundaries as shallow as 7-m depth. Boore and Joyner (1997) 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. Boore and Joyner (1997) 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 (Boore and Joyner, 1997).

Factor Influence Geometric spreading,     kij r
Introduces frequency dependent propagation delays.

Large-scale velocity structure,

    kij V D 3
Horizontal and vertical velocity and density gradients and velocity and density discontinuities produce frequency dependent phase shifts. Near-surface resonant responses,

    kij L
Interactions of shear-wave arrivals of varying angles of incidence and directions produce frequency dependent phase shifts.

Nonlinear soil responses,

    kij N u , (equivalent linear),    kij N u t , ( fully nonlinear)
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,

    kij Q
Linear hysteretic behavior produces frequency-dependent velocity dispersion that produces frequency dependent phase variations.

Scattering,

    kij S
The scattering strength and scattering characteristics determine propagation distances required to randomize the phase of shear waves as a function of frequency.

Anisotropy,

    kij A
Complicates shear-wave polarizations and modifies radiation pattern polarizations.

Topography,

    kij T
Complicates phase as a function of topographic length scale and near-surface velocities.

    kij )
A common site-response estimation method is to use the horizontal-to-vertical (H/V) spectral ratio method with shear waves [START_REF] Lermo | Site effect evaluation using spectral ratios with only one station[END_REF] to test for site resonances. The H/V method is similar to the receiver-function method of [START_REF] Langston | Structure under Mount Ranier, Washington, inferred from teleseismic body waves[END_REF].

Several investigations have shown the H/V approach provides robust estimates of resonant frequencies (e.g., [START_REF] Field | A comparison and test of various site response estimation techniques including threee that are not reference site dependent[END_REF][START_REF] Castro | S-wave site-response estimates using horizontal-to-vertical spectra ratios[END_REF][START_REF] Tsubio | Verification of horizontal-to-vertical spectral-ratio technique for estimate of site response using borehole seismographs[END_REF] although absolute amplification factors are less well resolved [START_REF] Castro | S-wave site-response estimates using horizontal-to-vertical spectra ratios[END_REF]Bonilla 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, [START_REF] Field | A modified ground-motion attenaution relationship for southern California that accounts for detailed site classification and a basin-depth effect[END_REF] found that the basin effect had a stronger influence on peak acceleration than detailed geology used to classify site responses. [START_REF] Hartzell | Variability of site response in Seattle, Washington[END_REF] 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. [START_REF] Chavez-Garcia | Lateral propagation effects observed at Parkway, New Zealand. A case history to compare 1D versus 2D site effects[END_REF] 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 [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogenous Earth[END_REF] 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 [START_REF] Frankel | Finite difference simulations of seismic scattering: implications for the propagation of short-period seismic waves in the crust and models of crustal heterogeneity[END_REF], most well-log data [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogenous Earth[END_REF], the frequency dependence of shear-wave attenuation [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogenous Earth[END_REF], and envelope broadening of shear waves with distance [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogenous Earth[END_REF]. 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 [START_REF] Frankel | Finite difference simulations of seismic scattering: implications for the propagation of short-period seismic waves in the crust and models of crustal heterogeneity[END_REF] fractal, self-similar velocity fluctuations are described with an autocorrelation function, P, of the form,

    P k a k a r n r n   1 (8)
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 [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogenous Earth[END_REF]. 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 [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogenous Earth[END_REF]. 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 [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogenous Earth[END_REF]. 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 [START_REF] Ishihara | Soil Behavior in Earthquake Geotechnics[END_REF]. While there has been evidence of nonlinear soil response in surface strong motion recordings [START_REF] Field | Nonlinear ground-motion amplification by sediments during the 1994 Northridge earthquake[END_REF][START_REF] Cultera | Nonlinear soil response in the vicinity of the Van Norman Complex following the 1994 Northridge, California, earthquake[END_REF], 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 [START_REF] Chang | Development of shear modulus reduction curves based on Lotung downhole ground motion data[END_REF]Wen et al., 1995, Ghayamghamain and[START_REF] Ghayamghamain | On the characteristics of non-linear soil response and dynamic soil properties using vertical array data in Japan[END_REF][START_REF] Satoh | Nonlinear behavior of soil sediments identified by using borehole records observed at the Ashigara Valley, Japan[END_REF][START_REF] Satoh | Nonlinear behavior of scoria evaluated from borehole records in eastern Shizuoka prefecture, Japan[END_REF][START_REF] Satoh | Inversion of strain-dependent nonlinear characteristics of soils using weak and strong motions observed by borehole sites in Japan[END_REF]. 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. [START_REF] Joyner | Calculation of nonlinear ground response in earthquakes[END_REF] 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. [START_REF] Archuleta | Nonlinearity in observed and computed accelerograms[END_REF] and [START_REF] Bonilla | Computation of linear and nonlinear response for near field ground motion[END_REF] 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 [START_REF] Joyner | Calculation of nonlinear ground response in earthquakes[END_REF]. 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 [START_REF] Archuleta | Nonlinearity in observed and computed accelerograms[END_REF]. [START_REF] Lomnitz | Seismic coupling of interface modes in sedimentary bains: A recipe for distaster[END_REF] 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. [START_REF] Lomnitz | Seismic coupling of interface modes in sedimentary bains: A recipe for distaster[END_REF] 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 [START_REF] O'connell | Influence of 2D Soil Nonlinearity on Basin and Site Responses[END_REF]. Anisotropy complicates polarizations of shear waves. [START_REF] Coutant | Observations of shallow anisotropy on local earthquake records at the Garner Valley, Southern California, downhole array[END_REF] showed that shallow (< 200 m) shear-wave anisotropy strongly influences surface polarizations of shear waves for frequencies < 30 Hz. [START_REF] Chapman | Ray tracing in azimuthally anisotropic media-II. Quasi-shear wave coupling[END_REF] 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. [START_REF] Chapman | Ray tracing in azimuthally anisotropic media-II. Quasi-shear wave coupling[END_REF] 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 [START_REF] Kawasaki | Radiation pattern of body waves due to the seismic dislocation occurring in an anisotropic source medium[END_REF].

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 [START_REF] Pitarka | Simulation of near-fault strong-ground motion using hybrid Green's functions[END_REF] or up to 10 Hz [START_REF] Zeng | Evaluation of numerical procedures for simulating nearfault long-period ground motions using the Zeng method[END_REF]. 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 [START_REF] Frankel | How does the ground shake[END_REF]. Topography can substantially influence peak ground motions [START_REF] Boore | A note of the effect of simple topography on seismic SH waves[END_REF][START_REF] Boore | The effect of simple topography on seismic waves: Implications for the acceleration recorded at Pacoima Dam, San Fernando Valley, California[END_REF]. [START_REF] Schultz | Enhanced backscattering of seismic waves from irregular interfaces[END_REF] 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. [START_REF] Bouchon | Effect of three-dimensional topography on seismic motion[END_REF] 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 [START_REF] Spudich | Directional topographic site response at Tarzana observed in aftershocks of the 1994 Northrige, Calfornia, earthquake: Implications for mainshock motions[END_REF]. 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 [START_REF] Hutchings | Ground-motion variability at the Highway 14 and I-5 interchange in the northern San Fernando Valley[END_REF]. Topography has a significant influence on longer-period amplification and groundshaking durations. [START_REF] O'connell | Influence of dip and velocity heterogeneity on reverse-and normal-faulting rupture dynamics and near-fault ground motions[END_REF] 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

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. [START_REF] Kramer | Geothechnical Earthquake Engineering[END_REF]. 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. [START_REF] Seed | Influence of soil conditions on ground motions during earthquakes[END_REF]), 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, [START_REF] Darragh | The site response of two rock and soil station pairs to strong and weak ground motion[END_REF] observed an amplification reduction at the Treasure Island soft soil site in San Francisco. [START_REF] Beresnev | Nonlinear site response -a reality?[END_REF] 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 and[START_REF] Aguirre | Nonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port Island, Kobe, during the Hyogo-ken Nanbu Earthquake[END_REF]. On the other hand, [START_REF] Darragh | The site response of two rock and soil station pairs to strong and weak ground motion[END_REF] 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 [START_REF] Field | Nonlinear sediment response during the 1994 Northridge earthquake: observations and finite-source simulations[END_REF]. 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;[START_REF] Darragh | The site response of two rock and soil station pairs to strong and weak ground motion[END_REF][START_REF] Field | Nonlinear ground-motion amplification by sediments during the 1994 Northridge earthquake[END_REF][START_REF] Beresnev | Nonlinearity at California generic soil sites from modeling recent strongmotion data[END_REF]. The accuracy of this approximation strongly depends on near surface rock weathering or topography complexity [START_REF] Steidl | What is A Reference Site? Bull[END_REF]. 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 [START_REF] Oglesby | A faulting model for the 1992 Petrolia earthquake: Can extreme ground acceleration be a source effect?[END_REF]. 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. [START_REF] Zeghal | Analysis of Site Liquefaction Using Earthquake Records[END_REF][START_REF] Aguirre | Nonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port Island, Kobe, during the Hyogo-ken Nanbu Earthquake[END_REF][START_REF] Satoh | Inversion of strain-dependent nonlinear characteristics of soils using weak and strong motions observed by borehole sites in Japan[END_REF][START_REF] Assimaki | Inverse analysis of weak and strong motion borehole array data from the Mw7.0 Sanriku-Minami earthquake[END_REF][START_REF] Assimaki | A Wavelet-based Seismogram Inversion Algorithm for the In Situ Characterization of Nonlinear Soil Behavior[END_REF][START_REF] Bonilla | Nonlinear site response evidence of K-NET and KiK-net records from the 2011 off the Pacific coast of Tohoku Earthquake[END_REF].

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. 4.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 [START_REF] Ishihara | Soil Behavior in Earthquake Geotechnics[END_REF][START_REF] Kramer | Geothechnical Earthquake Engineering[END_REF][START_REF] Beresnev | Nonlinear site response -a reality?[END_REF]

) = 1 + | |
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,

- 2 = - 2
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 [START_REF] Aguirre | Nonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port Island, Kobe, during the Hyogo-ken Nanbu Earthquake[END_REF]. 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. [START_REF] Iai | Response of a dense sand deposit during 1993 Kushiro-Oki Earthquake[END_REF], [START_REF] Archuleta | Direct observation of nonlinear soil response in acceleration time histories[END_REF][START_REF] Bonilla | Hysteretic and Dilatant Behavior of Cohesionless Soils and Their Effects on Nonlinear Site Response: Field Data Observations and Modeling[END_REF] 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 [START_REF] Towhata | Modeling Soil Behavior Under Principal Axes Rotation[END_REF] 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

[ ] = K ( ) ( ) + ( ) ( ) ( )
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 [START_REF] Konder | A hyperbolic stress-strain formulation for sands[END_REF] 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,

( ) = {1 1 0} ( ) = {cos -cos sin } = ( -1)
where ∆θ = π I ⁄ is the angle between each spring as shown in Figure 4.4. [START_REF] Towhata | Modeling Soil Behavior Under Principal Axes Rotation[END_REF] 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)). 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 [START_REF] Bonilla | Computation of linear and nonlinear response for near field ground motion[END_REF]. They are called generalized Masing rules because its formulation contains [START_REF] Pyke | Nonlinear soil model for irregular cyclic loadings[END_REF] 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 [START_REF] Ishihara | Modelling of stress-strain relations of soils in cyclic loading[END_REF] 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

= 1 1 + | | ⁄
where γ = τ G ⁄ 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

= ( ) = ⁄ 1 + | ⁄ |
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 [START_REF] Vucetic | Normalized behavior of clay under irregular cylic loading[END_REF][START_REF] Kramer | Geothechnical Earthquake Engineering[END_REF]. 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. [START_REF] Pyke | Nonlinear soil model for irregular cyclic loadings[END_REF][START_REF] Li | Dynamic skeleton curve of soil stress-strain relation under irregular cyclic loading Earthquake research in China[END_REF]. 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 [START_REF] Archuleta | Nonlinearity in observed and computed accelerograms[END_REF]. 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,

lim ⟶ ( )| | - ≤ ≤ |∞|
where is the specified finite or infinite strain condition, and correspond to the turning point and the hysteresis shape factor at the jth unloading or reloading; and

( ) is the sign of the strain rate. Thus,

= lim ⟶ ( )| | - + where = (| |)
, and ( , ) is the turning point pair at the jth reversal. Replacing the functional form of the backbone (the hyperbolic model) and after some algebra we have,

= - | -| | - - ( -) ≤ ≤ |∞|
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 [START_REF] Pyke | Nonlinear soil model for irregular cyclic loadings[END_REF], while → is similar to some extent to a method discussed in [START_REF] Li | Dynamic skeleton curve of soil stress-strain relation under irregular cyclic loading Earthquake research in China[END_REF]. In the following section, we will see an example of application of this soil constitutive model [START_REF] Towhata | Modeling Soil Behavior Under Principal Axes Rotation[END_REF]Iai et al., 1990aIai et al., , 1990b) ) together with the Generalized Masing hysteresis operator [START_REF] Bonilla | Computation of linear and nonlinear response for near field ground motion[END_REF].

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 [START_REF] Holzer | Dynamics of liquefaction during the 1987 Superstition Hills, California, earthquake[END_REF]. 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 [START_REF] Matasovic | Analysis of seismic records obtained on November 24, 1987 at the Wildlife Liquefaction Array[END_REF]. 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 [START_REF] Holzer | Dynamics of liquefaction during the 1987 Superstition Hills, California, earthquake[END_REF]. 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). [START_REF] Zeghal | Analysis of Site Liquefaction Using Earthquake Records[END_REF] 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,

(ℎ, ) = 1 2 ℎ[ (0, ) + (ℎ, )] (ℎ, ) = (0, ) + ℎ [ ( , ) -(0, )] ̅ ( ) = (ℎ, ) -(0, ) 2
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 [START_REF] Holzer | Dynamics of liquefaction during the 1987 Superstition Hills, California, earthquake[END_REF]. 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. Fig. 4.5. Wildlife Refuge station that recorded the 1987 Superstition Hills earthquake both acceleration and pores pressure time histories (left). Computed stress and strain time histories according to [START_REF] Zeghal | Analysis of Site Liquefaction Using Earthquake Records[END_REF], stress-strain loops and stress path history reconstitution (right).

Using the stress and stress time histories at GL-2.9 m computed earlier, [START_REF] Bonilla | Hysteretic and Dilatant Behavior of Cohesionless Soils and Their Effects on Nonlinear Site Response: Field Data Observations and Modeling[END_REF] 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.

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 [START_REF] Field | Nonlinear ground-motion amplification by sediments during the 1994 Northridge earthquake[END_REF][START_REF] Archuleta | Direct observation of nonlinear soil response in acceleration time histories[END_REF][START_REF] Seed | Analyses of ground motions at Union Bay, Seattle, during earthquakes and distant nuclear blasts[END_REF][START_REF] Beresnev | Nonlinear soil amplification: Its corroboration in Taiwan[END_REF][START_REF] Beresnev | Properties of vertical ground motions[END_REF] 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 [START_REF] Bonilla | 1D and 2D linear and non linear site response in the Grenoble area[END_REF] implementation of 2D nonlinear wave propagation. The constitutive equation implemented corresponds to the strain space multishear mechanism model developed by [START_REF] Towhata | Modeling Soil Behavior Under Principal Axes Rotation[END_REF] 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 Liu (2006). Linear hysteretic damping (Q) was implemented using the method of [START_REF] Liu | Efficient modeling of Q for 3D numerical simulation of wave propagation[END_REF]. The horizontal-and vertical component plane waves are inserted in the linear viscoelastic portion of the 2D with a userselectable range of incident angles. The Kyoshin Network, or Kik-net, in Japan (Fujiwara et al., 2005), has recorded numerous earthquakes with ground motion data recorded at surface and at depth in underlying rock and soil. We use the recording at Kik-net station IWTH25, where a 3.8g peak vertical acceleration was recorded (Aoi et al., 2008). Analyses of the combined downhole and surface ground motions from IWTH25 provide an opportunity to evaluate several strategies to estimate vertical ground motions since a P-and S-wave velocity profile is available to the bottom of the borehole at 260m (Aoi et al., 2008). Station IWTH25 is located in a region of rugged topography adjacent to a high-gradient stream channel on a fluvial terrace. The ruggest topography reflects the station's hangingwall location relative to the reverse fault . Station IWTH25 is located near a region of large slip along strike and updip of the hypocenter. Consequently, IWTH25 is subjected to significant rupture directivity and near-fault radiation associated with strong gradients of slip and rupture velocity on the portions of the fault close to the station (Miyazaki et al., 2009). The IWTH25 ground motion has been of particular interest because of the extreme peak vertical acceleration (3.8g) and peculiar asymmetric amplitudes distribution of the vertical accelerations (Aoi et al., 2008;[START_REF] O'connell | Assessing Ground Shaking[END_REF]Hada et al., 2009;Miyazaki et al., 2009;Yamada et al, 2009a); the upward vertical acceleration is much larger than the downward direction, although in the borehole record at a depth of 260 m at the same site, the upward and downward accelerations have symmetric amplitudes (Aoi et al., 2008) 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. It is important to mention some factors that are not explicitly accounted for in the approach of [START_REF] Bonilla | 1D and 2D linear and non linear site response in the Grenoble area[END_REF]. 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.

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 [START_REF] Campbell | Near-source attenuation of peak horizontal acceleration[END_REF][START_REF] Boore | Estimation of response spectra and peak accelerations from western North American earthquakes: An interim report[END_REF]. 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 Abrahamson and Silva (2008) and [START_REF] Chiou | An NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF], Z1.5 in Choi et al. (2005) and [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF]Z2.5 in Campbell andBozorgnia (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 [START_REF] Stafford | An Evaluation of the Applicability of the NGA Models to Ground Motion Prediction in the Euro-Mediterranean Region[END_REF][START_REF] Shoja-Taheri | A Test of the Applicability of the NGA Models to the Strong-Motion Data in the Iranian Plateau[END_REF]. 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 (Chiou et al., 2008). 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 Abrahamson and Silva (2008), [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF], [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF], Chiou andYoungs (2008), and[START_REF] Idriss | An NGA Empirical Model for Estimating the Horizontal Spectral Values Generated by Shallow Crustal Earthquakes[END_REF]. 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 [START_REF] Boore | GMRotD and GMRotI: Orientation-Independent Measures of Ground Motion[END_REF].

The NGA account for these ground motion factors using the general form,

    1 2 3 4 5 6 7 8 9 ln N N REF Source source site HW main Y A A M A M M A ln R C A R A F A F A F A F                         (9) and,   10 , 30 , lnY A M Vs   ( 10
)
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 [START_REF] Somerville | Differences in earthquake source and ground motion characteristics between surface and buried earthquakes[END_REF]) 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 Abrahamson and Silva (2008) and [START_REF] Chiou | An NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF], Z1.5 in Choi et al. (2005) and [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF], and Z2.5 in [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF], 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) [START_REF] Bommer | Why do modern probabilistic seismic-hazard analyses often lead to increased hazard estimates?[END_REF]. 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 (Boore et al., 1997;[START_REF] Sadigh | Atenuation relationships for shallow crustal earthquakes based on California strong motion data[END_REF][START_REF] Campbell | Empirical near-source attenuation relationships for horizontal and vertical components of peak ground acceleration, peak ground velocity, and pseudo-absolute acceleration response spectra[END_REF]. 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 (Abrahamson et al., 2008). 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 [START_REF] Kaklamanos | Model validations and comparisons of the next generation attenuation of ground motions (NGA-West) project[END_REF]. Thus, the use of Vs30 represents a significant improvement over previous generations of GMPEs that use a simple qualitative site classification scheme. [START_REF] Kaklamanos | Model validations and comparisons of the next generation attenuation of ground motions (NGA-West) project[END_REF] 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, Boore et al. (1997) 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). [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] used the Choi and Stewart (2005) linear amplification coefficients to remap observed response spectra to a reference Vs30=760 m/s. [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF] used 1D nonlinear soil amplification simulation results of [START_REF] Walling | Nonlinear site amplification factors for constraining the NGA models[END_REF] 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. Abrahamson and Silva (2008) use an approach similar to [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF]. [START_REF] Chiou | An NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF] do not explicitly specify how the coefficients for linear and nonlinear amplification were constrained or obtained. Thus, [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] remap observed response spectra prior to regression using the linear coefficients from Choi and Stewart (2005), [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF] and Abrahamson and Silva (2008) remap observed response spectra prior to regression using the nonlinear coefficients from [START_REF] Walling | Nonlinear site amplification factors for constraining the NGA models[END_REF], and it is not clear what [START_REF] Chiou | An NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF] did. We use [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] and [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF] 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 [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] 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 [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF], 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. 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 (Campbell andBozorgnia, 2008 and[START_REF] Abrahamson | Effects of rupture directivity on probabilistic seismic hazard analysis[END_REF]Silva, 2008), 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 [START_REF] Bonilla | Hysteretic and Dilatant Behavior of Cohesionless Soils and Their Effects on Nonlinear Site Response: Field Data Observations and Modeling[END_REF]. Chiou et al., 2008) 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 [START_REF] Idriss | An NGA Empirical Model for Estimating the Horizontal Spectral Values Generated by Shallow Crustal Earthquakes[END_REF] 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 [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] 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. Boore et al. (1997) 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 andFig Ground motion acceleration at high frequency scales in proportion to dynamic stress drop [START_REF] Boore | Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra[END_REF]. 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 [START_REF] Aki | Strong-motion seismology[END_REF][START_REF] Shaw | Constant stress drop from small to great earthquakes in magnitude-area scaling[END_REF]. Thus, as average slip increases with magnitude [START_REF] Somerville | Characterizing crustal earthquake slip models for the prediction of strong ground motion[END_REF][START_REF] Mai | Source scaling properties from finite-fault-rupture models[END_REF][START_REF] Mai | On scaling of fracture energy and stress drop in dynamic rupture models: Consequences for near-source ground motions, Earthquakes: Radiated Energy and the Physics of Faulting[END_REF] average rise time must also increase with increasing magnitude. [START_REF] Somerville | Magnitude scaling of the near fault rupture directivity pulse[END_REF] 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. [START_REF] Mai | On scaling of fracture energy and stress drop in dynamic rupture models: Consequences for near-source ground motions, Earthquakes: Radiated Energy and the Physics of Faulting[END_REF] 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 [START_REF] Mai | On scaling of fracture energy and stress drop in dynamic rupture models: Consequences for near-source ground motions, Earthquakes: Radiated Energy and the Physics of Faulting[END_REF] is consistent with average stress drop that is constant with moment. The [START_REF] Mai | On scaling of fracture energy and stress drop in dynamic rupture models: Consequences for near-source ground motions, Earthquakes: Radiated Energy and the Physics of Faulting[END_REF] results for maximum stress drop are consistent with first-order constraints on stochastic aspects of seismic source properties [START_REF] Andrews | A stochastic fault model, 2, Time-dependent case[END_REF][START_REF] Boore | Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra[END_REF][START_REF] Frankel | High-frequency spectral falloff of earthquakes, fractal dimension of complex rupture, b value, and the scaling strength on faults[END_REF]. 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 [START_REF] Andrews | A stochastic fault model, 2, Time-dependent case[END_REF][START_REF] Frankel | High-frequency spectral falloff of earthquakes, fractal dimension of complex rupture, b value, and the scaling strength on faults[END_REF]. 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 [START_REF] Aki | Strong-motion seismology[END_REF][START_REF] Shaw | Constant stress drop from small to great earthquakes in magnitude-area scaling[END_REF], 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 Boore et al. (1997) 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 Boore et al. (1997) response-spectral-magnitude-period-dependent results are more consistent with available seismological constraints. It is important to understand why. Boore et al. (1997) implements Vs30 site factors in a quite different manner than the four NGA GMPEs that use Vs30. Boore et al. (1997) 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, [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] use a similar regression approach for Boore et al. (1997). Since [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] 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 Boore et al. (1997); a non-linear regression will operate on the largest signals. The deterministic linear amplification function in [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] 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 Boore et al. (1997) 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 [START_REF] Idriss | An NGA Empirical Model for Estimating the Horizontal Spectral Values Generated by Shallow Crustal Earthquakes[END_REF] does not find credible. In contrast, Boore et al. (1997) 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 [START_REF] Bommer | On the selection of ground-motion prediction equations for seismic hazard analysis[END_REF] and [START_REF] Kaklamanos | Model validations and comparisons of the next generation attenuation of ground motions (NGA-West) project[END_REF], which clearly establish that NGA is a significant improvement for a wide range of applications than previous generation GMPEs, including Boore et al. (1997). 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 andFigure 5.3c), providing a vivid real-world example of epistemic uncertainty.

Earthquake

The site amplification approach used in NGA is discussed by [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF], "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 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. [START_REF] Kaklamanos | Model validations and comparisons of the next generation attenuation of ground motions (NGA-West) project[END_REF] 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. Boore and Joyner (1997) suggested that quarter-wavelength approximation of [START_REF] Joyner | The effect of Quaternary alluvium on strong ground motion in the Coyote Lake, California earthquake of 1979[END_REF] 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. [START_REF] Thompson | Multiscale site-response mapping: a case study of Parkfield, California[END_REF] 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 [START_REF] Stephenson | Blind shear-wave velocity comparison of ReMi and MASW results with boreholes to 200 m in Santa Clara Valley: Implications for earthquake ground-motion assessment[END_REF][START_REF] Boore | Comparisons of shear-wave slowness in the Santa Clara Valley, California, using blind interpretations of data from invasive and noninvasive methods[END_REF][START_REF] Miller | Advances in near-surface seismology and ground-penetrating radar: Introduction[END_REF][START_REF] O'connell | Interferometric multichannel analysis of surface waves (IMASW)[END_REF], 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.

5.3

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 [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF].

Unfortunately, some of the largest urban populations in the world are located within basins containing active faults, including many parts of Japan [START_REF] Kawase | The cause of the damage belt in Kobe: "The basin-edge effect," constructive interference of the direct S-wave with the basin-induced diffracted/Rayleigh waves[END_REF][START_REF] Pitarka | Three-dimensional simulation of the near-fault ground motions for the 1995 Hyogo-Nanbu (Kobe), Japan, earthquake[END_REF][START_REF] Nied | Off the Pacific Coast of Tohoku Earthquake, Strong Ground Motion[END_REF], the Los Angeles and other basins in southern California [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF], and Seattle, Washington, [START_REF] Frankel | Sedimentary basin effects in Seattle, Washington: Ground-motion observations and 3D simulations[END_REF]. Consequently, estimation of long-period ground motions in sedimentary basins associated with near-fault faulting is an important practical need. Choi et al. (2005) 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 Choi and Stewart (2005) 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, Choi et al. (2005) 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 [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF] 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. [START_REF] Day | 3-D ground motion simulation in basins: Final report for Lifelines Project 1A03[END_REF] 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 [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF] 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 Choi et al. (2005) 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 [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF]. 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 (Choi et al, 2005). [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF] 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, [START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF] used the parametric 3D synthetic basin-depth model from [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF]. [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF] 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. [START_REF] Chiou | An NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF] 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, [START_REF] Chiou | An NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF] 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. Abrahamson and Silva ( 2008) 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. 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 (Abrahamson and Silva, 2008;[START_REF] Campbell | NGA Ground Motion Model for the Geometric Mean Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response Spectra for Periods Ranging from 0.01 to 10 s[END_REF][START_REF] Chiou | An NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra[END_REF], and [START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] for periods of 1 s. The NGA results are modified to account for rupture directivity using [START_REF] Rodriguez-Marek | An empirical geotechnical seismic site response procedure[END_REF], Spudich andChiou (2008), and[START_REF] Rowshandel | Directivity correction for the next generation attenuation (NGA) relations[END_REF] 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 [START_REF] Bécel | Moho, crustal architecture and deep deformation under the North Marmara Trough, from the SEISMARMARA Leg 1 offshore-onshore reflection-refraction survey[END_REF][START_REF] Bayrakci | Approach to the complex 3D upper-crustal seismic structure of the Sea of Marmara by artificial source tomography on a 2D grid of OBS[END_REF], tomographically assessed top of basement contours [START_REF] Bayrakci | Approach to the complex 3D upper-crustal seismic structure of the Sea of Marmara by artificial source tomography on a 2D grid of OBS[END_REF], seismic reflection profiles [START_REF] Carton | Seismic imaging of the threedimensional architecture of the Çınarcık Basin along the North Anatolian Fault[END_REF]Kurt and Yucesoy, 2009), Boguer gravity profiles [START_REF] Ates | Structural interpretation of the Marmara region, NW Turkey, from aeromagnetic, seismic and gravity data[END_REF], geologic mapping [START_REF] Okyar | Late quaternary seismic stratigraphy and active faults of the Gulf of İzmit (NE Marmara Sea)[END_REF]) and fault mapping [START_REF] Armijo | Asymmetric slip partitioning in the Sea of Marmara pull-apart: a clue to propagation processes of the North Anatolian Fault?[END_REF]. 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 [START_REF] Bayrakci | Approach to the complex 3D upper-crustal seismic structure of the Sea of Marmara by artificial source tomography on a 2D grid of OBS[END_REF]. 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 [START_REF] Hartzell | Effects of 3D random correlated velocity perturbations on predicted ground motions[END_REF] 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 [START_REF] Olsen | Estimation of Q for long-period(>2 sec) waves in the Los Angeles basin[END_REF]. The fourth-order finite-difference code employs the efficient and accurate viscoelastic formulation of [START_REF] Liu | Efficient modeling of Q for 3D numerical simulation of wave propagation[END_REF] that accurately A kinematic representation of finite fault rupture is used where fault slip (displacement), rupture time, and rise time are specified at each finite-difference grid node intersected by the fault. The 3D viscoelastic fourth-order finite-difference method of [START_REF] Liu | The effect of a low-velocity surface layer on simulated ground motion[END_REF]2006) 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 [START_REF] Andrews | Dynamic simulation of spontaneous rupture with heterogeneous stress drop[END_REF]. 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 [START_REF] Andrews | Dynamic simulation of spontaneous rupture with heterogeneous stress drop[END_REF] 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 [START_REF] Somerville | Characterizing crustal earthquake slip models for the prediction of strong ground motion[END_REF] and  2 ( is angular frequency) displacement spectra in the far-field. [START_REF] Oglesby | Stochastic fault stress: Implications for fault dynamics and ground motion[END_REF] and [START_REF] Schmedes | Correlation of earthquake source parameters inferred from dynamic rupture simulations[END_REF] 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. [START_REF] Oglesby | Stochastic fault stress: Implications for fault dynamics and ground motion[END_REF], [START_REF] Guatteri | Strong ground motion prediction from stochastic-dynamic source models[END_REF][START_REF] Schmedes | Correlation of earthquake source parameters inferred from dynamic rupture simulations[END_REF] 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 [START_REF] Petersen | Fault displacement hazard for strike-slip faults[END_REF] 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 [START_REF] Liu | Efficient modeling of Q for 3D numerical simulation of wave propagation[END_REF] 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 [START_REF] Boore | GMRotD and GMRotI: Orientation-Independent Measures of Ground Motion[END_REF] (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 [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF]. 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 [START_REF] Rodriguez-Marek | An empirical geotechnical seismic site response procedure[END_REF], [START_REF] Spudich | Directivity in NGA earthquake ground motions: Analysis using isochrone theory[END_REF], 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 [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF]. 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 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 [START_REF] Rowshandel | Directivity correction for the next generation attenuation (NGA) relations[END_REF] 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 [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF] 3D synthetic calculations to represent basin amplification in younger shallower basins. Thus, it is interesting to compare the [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF] 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 [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF] response ratios of sites 1-3 to site 4 using the Z1.5 values from the 3D simulation model in the [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF] 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; [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF] 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 [START_REF] O'connell | Influence of 2D Soil Nonlinearity on Basin and Site Responses[END_REF] 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 [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF] 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 Choi et al. (2005), [START_REF] Day | Model for basin effects on long-period response spectra in southern California[END_REF], 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 [START_REF] Jennings | Engineering seismology, In: "Earthquakes: Observation, Theory, and Interpretation[END_REF] 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 (Chiou et al., 2008;[START_REF] Kaklamanos | Model validations and comparisons of the next generation attenuation of ground motions (NGA-West) project[END_REF]. 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% [START_REF] Stephenson | Blind shear-wave velocity comparison of ReMi and MASW results with boreholes to 200 m in Santa Clara Valley: Implications for earthquake ground-motion assessment[END_REF][START_REF] Boore | Comparisons of shear-wave slowness in the Santa Clara Valley, California, using blind interpretations of data from invasive and noninvasive methods[END_REF]. 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 Choi et al. (2005) 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.

  the convolution of the time evolution of the slip-time functions, responses between the fault and the site (Figure3.1) as,

Fig. 3

 3 Fig. 3.1. Schematic diagram of finite-fault rupture ground motion calculations. Three discrete subfault elements in the summation are shown. Rings and arrows emanating from the hypocenter represent the time evolution of the rupture. The Green functions actually consist of eight components of ground motion and three components of site ground velocities. Large arrows denote fault slip orientation, which is shown as predominantly reverse slip with a small component of right-lateral strike slip. Hatched circles schematically represent regions of high stress drop.

Fig. 3

 3 Fig. 3.2. Schematics of line source orientations for strike-slip (a) and thrust faults (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).

Fig. 4

 4 Fig. 4.1. Hyperbolic model of the stress-strain space for a soil under cyclic loading. Initial loading curve has a hyperbolic form, and the loading and unloading phases of the hysteresis path are formed following Masing's criterion.

Figure

  Figure 4.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[START_REF] Ishihara | Soil Behavior in Earthquake Geotechnics[END_REF][START_REF] Kramer | Geothechnical Earthquake Engineering[END_REF][START_REF] Beresnev | Nonlinear site response -a reality?[END_REF]

Fig. 4

 4 Fig. 4.2. Borehole transfer functions computed at KiK-net station TTRH02 in Japan. The orange shaded area represents the 95% confident limits of the transfer function using weakmotion events (PGA < 10cm/s 2 ). The solid line is the transfer function computed using the October 2000 Tottori mainshock data.

Fig. 4

 4 Fig. 4.3. Surface and borehole records of the 1995 Kobe earthquake at Port Island (left), and the 1993 Kushiro-Oki earthquake at Kushiro Port (right). The middle panel shows the shear wave velocity distribution at both sites.

Fig. 4

 4 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 andIshihara, 1985).

  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.

Fig. 4

 4 Fig. 4.6. The top panel shows the computed strain time history at the middle of the borehole. Middle panels show the computed stress by trial-and-error using the multispring model in order to find the best dilatancy parameters. Bottom panels indicate the computed stress time history from acceleration records (after Bonilla et al., 2005).

  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 4.9d), acceleration seismograms (Figures 4.9b and 4.9eb), and velocity seismograms (Figures 4.9cand 4.9f). Since station IWTH25 is located in the deformed hangingwall of a reverse fault in rugged topography, it is clear that even these 2D nonlinear calculations are a crude approximation to the field conditions and complex incident wavefield associated with the finite fault rupture. However, the 2D nonlinear calculations summarized in Figs. 4.8 and 4.9 for station IWTH25 clearly show that the 2D P-SV nonlinear approach of[START_REF] Bonilla | 1D and 2D linear and non linear site response in the Grenoble area[END_REF] provide a sound basis to evaluate first-order nonlinear horizontal and vertical nonlinear responses, even for cases of extremely large incident accelerations and velocities.

Fig. 4

 4 Fig. 4.8. Observed and simulated surface ITWH25 horizontal response spectra (a,d), and acceleration (b,e), and velocity (c,f) time histories for the north (a-c) and east (d-f) components.

  Fig. 4.9. Observed and simulated surface ITWH25 vertical response spectra (a,d), and acceleration (b,e), and velocity (c,f) time histories using the north-vertical (a-c) and eastvertical (d-f) components.

Fig. 5

 5 Fig. 5.1. Boore and Atkinson (2008) amplification functions (a) and original Vs30=300 m/s (black) and "Vs30=915 m/s remapped observed" responses spectra (red) (b) for M=7.0, distance of 2 km, and PGA=0.45 g.

Fig. 5

 5 Fig. 5.2. Campbell and Bozorgnia (2008) amplification functions (a) and original Vs30=300 m/s (black) and "Vs30=915 m/s remapped observed" responses spectra (red) (b) for M=7.0, distance of 2 km, and PGA=0.45 g.

  Fig. 5.3. Boore and Atkinson (2008) (a) and Boore et al. (1997) (b) response spectra normalized by peak ground acceleration for Vs30=900 m/s. The geometric mean spectral accelerations from the three observed Vs30 > 900 m/s ground motions in Table6is compared to the mean[START_REF] Boore | Ground-Motion Prediction Equations for the Average Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods Between 0.01s and 10.0 s[END_REF] andBoore et al. (1997) estimates in (c).

  Fig. 5.4. North Anatolia Fault segments and sites for 3D ground motion modeling.

Fig. 5

 5 Fig. 5.5. Shear-wave (Vs) cross sections through the 3D velocity model along profiles shown in map view in Fig. 5.4.

Fig. 5

 5 Fig. 5.6. Near-median bilateral three-segment rupture synthetic velocity seismograms for the five sites shown in Figs. 5.4 and 5.5a.

  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).

Fig. 5

 5 Fig. 5.7. Site 4 3D synthetic and NGA GMPE mean response spectra (a) and the 3D/NGA GMPE ratio (b).

Fig. 5

 5 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.

Fig. 5

 5 Fig. 5.9. Mean 3D site-specific simulation 3D amplification and Day et al. (2008) 3D amplifications for sites 1-3 relative to reference site 4.

  Table 1 list factors influencing source amplitudes,

	source phase,	  S ij  . Table 2 lists factors influencing 

Table 2

 2

	. Seismic Source Phase Factors (    

ij ) Oglesby et al.

Table 3

 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

  velocity materials near the surface amplify ground motions for frequencies > /(4*h), where h is the thickness of near-surface low velocity materials. Coupled interface modes can amplify and extend durations of ground motions.

	Nonlinear soil responses, G kij (equivalent linear)   u N , ,
	G	  u t
		G	u t

kij N , (fully nonlinear) Depending on the dynamic soil properties and pore pressure responses, nonlinear soil responses can decrease intermediate-and high-frequency amplitudes, amplify low-and high-frequency amplitudes, and extend or reduce duration of large amplitudes. The equivlanet linear approximation is   G u kij N , . The fully nonlinear form,   kij N , , can incorporate any time-dependent behavior such as pore-pressure responses. Frequency independent attenuation,

Table 4 .

 4 Seismic Wave Propagation Phase Factors (

  and reproducible using 1D nonlinear site response modeling[START_REF] O'connell | Assessing Ground Shaking[END_REF]. 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[START_REF] O'connell | Assessing Ground Shaking[END_REF]. In particular, application of the nonlinear approach of shear-modulus reduction advocated and tested byBersenev et al. (2002) to predict nonlinear vertical responses, failed to predict peak vertical accelerations in excess of 2g[START_REF] O'connell | Assessing Ground Shaking[END_REF]). 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 ofBonilla et al. (

  .Aoi et al. (2008) propose a conceptual model for this asymmetry. Their model uses a loose soil with nearly zero confining pressure near the surface. The soil particles separate under large downward acceleration, and in this quasi free-fall state, the downward accelerations at the surface only modestly exceed gravity. Conversely, large upward accelerations compact the soil and produce much larger upward accelerations. Aoi et al. (2008) report three cases of these anomalous large vertical acceleration amplifications in a search of 200,000 strong motion recordings.Hada et al. (2009) successfully reproduce the strong vertical asymmetric accelerations at IWTH25 with a simple 1D discrete element model, a model that is not a rigorous model of wave propagation.Yamada et al. (2009a) interpret the large upward spikes in acceleration as slapdown phases, which are also typically observed in near-field recordings of nuclear explosion tests. Our focus here is not the asymmetry of the IWTH25 vertical accelerations recorded at the surface, but showing that the simple total stress plane-

strain model of soil nonlinearity in

[START_REF] Bonilla | 1D and 2D linear and non linear site response in the Grenoble area[END_REF]

reproduces both the first-order peak horizontal and vertical velocities and accelerations and acceleration response spectra at station IWTH25 using the borehole motions at 260 m depth as inputs.

Yamada et al. (2009b)

conducted geophysical investigations at the site and found lower velocities in the top several meters than reported in

Aoi et al. (2008)

. Trail-and-error modeling was used to obtain the final refined velocity model consistent with the results of

Yamada et al. (2009b)

; a lowest-velocity first layer of about 2 m thickness and shear-wave velocity on the order of 200 m/s was required to produce the maximum horizontal spectral responses observed near 10 Hz.

Table 5 .

 5 NGA Near-Fault Strike-Slip Ground Motions

	(name)	Date (day,mon,yr)	M	Station	Vs30 (m/s)	JB Fault Distance (km)
	Parkfield	28 Jun. 1966 6.1	Cholame 2WA	185	3.5
	Imperial Valley	15 Oct. 1979 6.5 El Centro Array #7	212	3.1
	Superstition Hills	24 Nov. 1987 6.6	Parachute	349	1.0
	Erzincan	13 Mar. 1992 6.9	95 Erzincan	275	2.0
	Landers	28 Jun. 1993 7.3	Lucerne	665	1.1
	Kobe, Japan	16 Jan. 1995 6.9	KJMA	312	0.6
	Kocaeli, Turkey	17 Aug. 1999 7.4	Yarimca	297	2.6
	Kocaeli, Turkey	17 Aug. 1999 7.4	Sakarya	297	3.1
	Duzce, Turkey	12 Nov. 1999 7.1	Duzce	276	8.1
	Geometric Mean		6.9		299	2.1

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. The authors benefited from helpful discussions with David Boore, Joe Andrews, Paul Spudich, Art Frankel, Dave Perkins, Chris Wood, Ralph Archuleta, David Oglesby, Steve Day, Bruce Bolt, Rob Graves, Roger Denlinger, Bill Foxall, Larry Hutchings, Ned Field, Hiro Kanamori, Dave Wald, and Walt Silva. Shawn Larson provided the e3d software, and Paul Spudich provided isochrone software. Supported by U.S. Bureau of Reclamation Dam Safety Research projects SPVGM and SEIGM and USGS award no. 08HQGR0068. The National Information Center