Three-Dimensional Ground-Motion Simulations of 

Earthquakes for the Hanford Area, Washington 

By Arthur Frankel, Paul Thorne, and Alan Rohay

Open-File Report 2013–1289

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

U.S. Geological Survey, Reston, Virginia: 2014 

For more information on the USGS—the Federal source for science about the Earth,

its natural and living resources, natural hazards, and the environment—visit 

http://www.usgs.gov or call 1–888–ASK–USGS

For an overview of USGS information products, including maps, imagery, and publications,

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

Suggested citation:

Frankel, Arthur, Thorne, Paul, and Rohay, Alan, 2014, Three-dimensional ground-motion simulations of earthquakes

for the Hanford area, Washington: U.S. Geological Survey Open-File Report 2013-1289, 48 p., 

http://dx.doi.org/10.3133/ofr20131289. 

ISSN 2331-1258 (online)

Any use of trade, product, or firm names is for descriptive purposes only and does not imply 

endorsement by the U.S. Government.

Although this report is in the public domain, permission must be secured from the individual 

copyright owners to reproduce any copyrighted material contained within this report.
iii 

Contents

Introduction.................................................................................................................................................................... 1

Three-Dimensional Velocity Models .............................................................................................................................. 1

Finite-Difference Program and Grid Spacing ................................................................................................................. 2

Seismic Stations ............................................................................................................................................................ 3

Comparison of Synthetics with Data.............................................................................................................................. 4

Prosser Earthquake ................................................................................................................................................... 4

Wooded Island Swarm Events ................................................................................................................................... 6

Walla Walla Earthquake............................................................................................................................................. 7

Simulations for Magnitude 6.5–6.8 Earthquakes ........................................................................................................... 7

One-Dimensional Models and Reflectivity Method ....................................................................................................... 8

Rattlesnake Hills Fault, West, Moment Magnitude 6.8............................................................................................... 9

Rattlesnake Hills Fault, East, Moment Magnitude 6.7...............................................................................................11

Gable Mountain Fault, Moment Magnitude 6.6 .........................................................................................................11

Cascadia...................................................................................................................................................................11

Conclusions..................................................................................................................................................................12

Acknowledgments ........................................................................................................................................................13

References Cited..........................................................................................................................................................13

Figures

Figure 1. Map of depth (meters) to S-wave velocity (Vs) of 1,500 meters per second for the Hanford area, 

Washington.................................................................................................................................................................. 15

Figure 2. Cut-out view of three-dimensional (3D) velocity model that used depth to basement derived from gravityderived model (GDM) .................................................................................................................................................. 16

Figure 3. 3D views of the elevation of the top of basement (m) for the gravity-derived model (GDM; top) and the 

refraction-derived model (RDM; bottom....................................................................................................................... 17

Figure 4. Satellite image of Hanford area from GoogleEarthTM showing locations of strong-motion accelerometers 

(circles), broadband stations (squares), and earthquakes studied (stars) ................................................................... 18

Figure 5. Observed (black) and synthetic (red) velocity waveforms at station H4A for the Prosser earthquake ......... 19

Figure 6. One-dimensional (1D) velocity profiles at station H2E used in the modeling of the recorded earthquakes. 

The gravity-derived model (GDM) uses the depth to top of basement derived from gravity modeling......................... 20

Figure 7. Observed (black) and synthetic (red) velocity waveforms at station H4A for the Prosser earthquake ......... 21

Figure 8. Observed (black) and synthetic (red) velocity waveforms at station LTH for the Prosser earthquake ......... 22

Figure 9. Observed (black) and synthetic (red) velocity waveforms at station H2E for the Prosser earthquake ......... 23

Figure 10. Observed (black) and synthetic (red) velocity waveforms at station H3A for the Prosser earthquake. ...... 24

Figure 11. Observed (black) and synthetic (red) velocity waveforms at station HAWA for the Prosser earthquake. .. 25

Figure 12. Observed and synthetics velocity waveforms (NS component) at station HAWA for the Prosser

earthquake.................................................................................................................................................................. 26

Figure 13. Bias and standard deviation of the spectral accelerations from the observed recordings at 0.5 and 1.0 s 

with respect to the spectral accelerations from the synthetics for the Prosser earthquakes (stations H3A, H4A, H2E, 

and LTH)...................................................................................................................................................................... 27

Figure 14. Bias and standard deviation of observed versus synthetic spectral accelerations for the Prosser 

earthquake................................................................................................................................................................... 28

Figure 15. Bias and standard deviation of observed versus synthetic spectral accelerations for the Prosser 

earthquake using three-dimensional (3D) models with different depths to basement.................................................. 29
iv

Figure 16. Cumulative energy at station H4A for Prosser earthquake plotted with respect to time after origin

time for observed seismograms (black and red lines) compared to synthetics (green and blue lines) from

three-dimensional (3D) (top) and one-dimensional (1D) (bottom) velocity models ...................................................... 30

Figure 17. Bias and standard deviation of energy derived from observed seismograms compared with those

from the three-dimensional (3D) and one-dimensional (1D) synthetics, for the Prosser earthquake 

(stations H3A, H4A, H2E, and LTH). ........................................................................................................................... 31

Figure 18. Observed (black) and synthetic (red) velocity waveforms for the September 2011 Wooded Island 

earthquake................................................................................................................................................................... 32

Figure 19. Bias and standard deviation of the observed and synthetics spectral accelerations for the 

September 2011 Wooded Island earthquake. ............................................................................................................. 33

Figure 20. Observed (black) and synthetic (red) velocity waveforms for the October 2011 

Wooded Island earthquake.......................................................................................................................................... 34

Figure 21. Bias and standard deviation of the observed and synthetic spectral accelerations for the 

October 2011 Wooded Island earthquake ................................................................................................................... 35

Figure 22. Observed (black) and synthetic (red) velocity waveforms at station LTH for the 

Walla Walla earthquake............................................................................................................................................... 36

Figure 23. Map showing surficial S-wave velocity (Vs)(in meters per second) of three-dimensional model 

and surface projections (red rectangles) of rupture planes used in the simulations..................................................... 37

Figure 24. Final slip and rupture initiation times for one of the simulations for the Rattlesnake Hills West source. .... 38

Figure 25. S-wave velocity (Vs) profiles used in the one-dimensional (1D) simulations for the 

Rattlesnake Hills West and Gable Mountain sources (bottom panel shows shallow portion). ..................................... 39

Figure 26. Velocity seismograms for one of the simulations of a magnitude 6.8 earthquake on the 

Rattlesnake Hills West source ..................................................................................................................................... 40

Figure 27. Natural logarithms of ratios of spectral accelerations between the three-dimensional (3D) and onedimensional (1D) models for the magnitude 6.8 Rattlesnake Hills West source for the sites used for the probabilistic 

seismic hazard assessment......................................................................................................................................... 41

Figure 28. Velocity seismograms for one of the simulations of a M6.7 earthquake on the Rattlesnake Hills East 

source .......................................................................................................................................................................... 42

Figure 29. Natural logarithms of ratios of spectral accelerations between the three-dimensional and one-dimensional 

models for the magnitude 6.7 Rattlesnake Hills East source for the sites used for the probabilistic seismic hazard 

assessment. ................................................................................................................................................................ 43

Figure 30. Velocity seismograms for one of the simulations of a magnitude 6.6 earthquake on the 

Gable Mountain fault.................................................................................................................................................... 44

Figure 31. Natural logarithms of ratios of spectral accelerations between the three-dimensional and one-dimensional 

models for the magnitude 6.6 Gable Mountain source for the sites used for the probabilistic seismic hazard 

assessment. ................................................................................................................................................................ 45

Figure 32. Velocity seismograms for simulation for Cascadia plane-wave source Synthetics for three-dimensional 

model are in black, synthetics for one-dimensional models are in red, and vertical scale is in arbitrary units. ............ 46

Figure 33. Natural logarithms of ratios of spectral accelerations from the three-dimensional model relative to those 

from the one-dimensional model at each site for the Cascadia plane wave source..................................................... 47

Table

Table 1. Source parameters of earthquakes simulated in this study. .......................................................................... 48
v 

Conversion Factors and Datums

SI to Inch/Pound

Multiply By To obtain

Length

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi)

kilometer (km) 0.5400 mile, nautical (nmi) 

meter (m) 1.094 yard (yd) 

Flow rate

meter per second (m/s) 3.281 foot per second (ft/s) 

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as follows:

°F=(1.8×°C)+32

Datums

Vertical coordinate information is referenced to the North American Vertical Datum of 1988 (NAVD 88)

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83)

Altitude, as used in this report, refers to distance above the vertical datum.
vi

This page left intentionally blank 
1 

Three-Dimensional Ground-Motion Simulations of 

Earthquakes for the Hanford Area, Washington

By Arthur Frankel1, Paul Thorne2, and Alan Rohay2

Introduction

This report describes the results of ground-motion simulations of earthquakes using 

three-dimensional (3D) and one-dimensional (1D) crustal models conducted for the probabilistic 

seismic hazard assessment (PSHA) of the Hanford facility, Washington, under the Senior 

Seismic Hazard Analysis Committee (SSHAC) guidelines. The first portion of this report 

demonstrates that the 3D seismic velocity model for the area produces synthetic seismograms 

with characteristics (spectral response values, duration) that better match those of the observed 

recordings of local earthquakes, compared to a 1D model with horizontal layers. The second part 

of the report compares the response spectra of synthetics from 3D and 1D models for moment 

magnitude (M) 6.6–6.8 earthquakes on three nearby faults and for a dipping plane wave source 

meant to approximate regional S-waves from a Cascadia great earthquake. The 1D models are 

specific to each site used for the PSHA. The use of the 3D model produces spectral response 

accelerations at periods of 0.5–2.0 seconds as much as a factor of 4.5 greater than those from the 

1D models for the crustal fault sources. The spectral accelerations of the 3D synthetics for the 

Cascadia plane-wave source are as much as a factor of 9 greater than those from the 1D models. 

The differences between the spectral accelerations for the 3D and 1D models are most 

pronounced for sites with thicker supra-basalt sediments and for earthquakes on the Rattlesnake 

Hills fault and for the Cascadia plane-wave source. 

Three-Dimensional Velocity Models

Paul Thorne, Alan Rohay, and Stephen Reidel developed two 3D crustal models 

describing the P-wave velocity (Vp), S-wave velocity (Vs), and density for the Hanford area. 

The models were based on maps of thicknesses of the various geologic units in the area. Each 

geologic unit was assigned a Vp, Vs, and density. For the shallow units, these values were based 

on borehole velocity measurements. The parameters for the deeper layers were constrained from 

seismic refraction surveys and gravity data (for example, Burns and others, 2011). 

___________________

1

U.S. Geological Survey.

2

Pacific Northwest National Laboratory.
2 

The shallow portion of the Hanford region is characterized by a relatively thin layer 

(< 150 m) of Holocene to Miocene sediments on top of a sequence of Columbia River basalts

approximately 3–5 km thick. The depth to Vs of 1,500 m/s for the model, which is equivalent to 

a map of the thickness of these supra-basalt sediments, is shown in figure 1. There is a strong 

velocity contrast at the base of the surficial sediments and the top of the Columbia River basalts.

As we will demonstrate from the observations and simulations, these supra-basalt sediments can 

be a waveguide for high-frequency (1–2 hertz [Hz]) surface waves and can produce 

reverberations of nearly vertically propagating S-waves. The supra-basalt sediments include 

Holocene deposits of alluvium and loess, Pleistocene flood deposits (Hanford Formation), PlioPleistocene deposits, and the Ringold Formation of Miocene-Pliocene age (Geomatrix, 1996). 

The minimum Vs is 425 m/s. There are several different basaltic units, including the Saddle 

Mountains, Wanapum, and Grand Ronde basalts. Each of these units was assigned a Vp, Vs, and 

density. The model did not explicitly include the thin interbedded layers of sediment within the 

basalts. The velocities and density of the Saddle Mountains basalt were adjusted to account for 

the sediment interbeds. 

There is a large thickness of pre-Miocene sediments beneath the basalts, although 

substantial uncertainty exists about this thickness. Two models were produced. One was based 

on modeling of gravity data and had an inferred thickness of about 1–2 km for the sub-basalt 

sediments. This depth-to-basement model was provided by Stephen Reidel and is denoted as the 

gravity-derived model (GDM). A cut-out view of this model is shown in figure 2. The other 

model was based on seismic refraction data and the thicknesses of the sub-basalt sediments were 

about about 5–6 km (Glover, 1985). This model is called the refraction-derived model (RDM). 

The 3D configurations of the top of basement in the two models are shown in figure 3. There is a 

Vp and Vs gradient with depth within the sub-basalt sediments. The basement rock was assigned 

a Vp of 6,100 m/s and Vs of 3,500 m/s. The maximum depth of the model was 25 km. 

The finite-difference program used in this study requires a flat, free surface. The 3D 

models were developed so that the thicknesses of the supra-basalt sediments and the other 

geologic units at any location were preserved. In the 3D models, these thicknesses were applied 

starting at a free surface with a constant elevation. 

Qs values were determined from Vs using the formulas from Frankel and others (2007, 

2009), which were used in simulations for the Seattle basin and produced good fits to the 

observed spectral amplifications and waveforms. For Vs less than 1,000 m/s, Qs=0.16428*Vs -

14.2857, where Vs is in m/s. For Vs > 1,000 m/s, Qs=0.15*Vs. As will be shown later, these Qs 

values produce durations of synthetics similar to the data. Qp was taken to be 2*Qs. 

Finite-Difference Program and Grid Spacing

We used the 3D finite-difference computer program DISFD written by Pengcheng Liu of 

the Bureau of Reclamation (Liu and Archuleta, 2002). This program consists of a time-stepping 

algorithm to solve for the velocity vectors and stress tensors on a grid with complex material 

properties. The code uses fourth order finite-difference approximations for the partial derivatives 

in space and second order approximation for the partial derivative in time (Liu and Archuleta, 

2002). A key advantage of the code is that it allows for the grid spacing in the top portion of the 

model to be one-third of that in the bottom portion of the model. This saves substantial 

computing time, because the finer grid spacing needed for accuracy in the lower velocity layers 

near the surface does not have to be continued to deeper depths. The finite-difference code also 

allows for an increase of the vertical grid spacing with depth. The code uses the linear velocity-
3 

stress equations with frequency-independent Qp and Qs that can vary with position. Each source 

is specified by a seismic moment, rise time, start time, focal mechanism, and location in the grid. 

The code uses the message-passing interface (MPI) protocol so that multiple processors can be 

used. A typical run involved 96 processors on the Hydra Cluster in Menlo Park, California. 

For the top 500 m of the model, the horizontal and vertical grid spacing used in the finitedifference simulations was 25 m. The seismic velocities and densities in the top 500 m of the 

model were specified with a horizontal grid spacing of 100 m and vertical grid spacing of 25 m. 

The finite-difference code interpolated to a horizontal grid spacing of 25 m. At depths greater 

than 500 m, the vertical grid spacing increased to 75 m. At depths from 500 to 6,500 m, the 

model was developed with a horizontal grid spacing of 200 m and vertical grid spacing of 75 m.

At depths greater than 6,500 m, the model was specified with a vertical grid spacing of 225 m 

and horizontal spacing of 400 m. The finite-difference code interpolated to a horizontal grid 

spacing of 75 m at depths greater than 500 m. The vertical grid spacing used by the finitedifference code for depths greater than 8,000 m was 225 m.

The simulations are accurate up to 2 Hz. At the 425 m/s minimum Vs, 2 Hz corresponds 

to a wavelength of 213 m, which is 8.5 grid spacings (25 m for the top 500 m). For fourth order 

codes, a general guideline is that at least six grid points per wavelength are required for accuracy 

(Alford and others, 1974).

All synthetics and data in this study are low-pass filtered with a corner frequency of 2.0 

Hz. The observed seismograms from the accelerometers also are high-pass filtered at 0.5 Hz, 

because of poor signal-to-noise less than 0.5 Hz. This high-pass filtering did not significantly 

affect determinations of spectral accelerations at 0.5 and 1.0 s. Seismograms from the broadband 

sensors did not require high-pass filtering.

Seismic Stations

Stations and earthquakes modeled in this study are shown in figure 4. The key stations 

are the five strong-motion accelerometers currently operated by the Pacific Northwest Seismic 

Network (PNSN) of the University of Washington and originally installed by the Pacific 

Northwest National Laboratory (PNNL). We also used recordings from the broadband 

instruments at station HAWA (U.S. National Seismic Network), station LTH, and two stations 

from the Earthscope Transportable Array. The data from the broadband instruments were 

obtained from the IRIS Data Management Center. The seismograms from these instruments were 

converted to ground velocity using the poles, zeroes, and gains listed in their instrument files. 

The data search was assisted by the database compiled by Adrian Rodriguez-Marek for the 

Hanford PSHA (A. Rodriguez-Marek, written comm., 2013). 
4 

Comparison of Synthetics with Data

Prosser Earthquake

The PNSN earthquake database was systematically searched to identify earthquakes that 

were larger than duration magnitude (Md) 3.0 that occurred within the 3D model; that were 

recorded by the strong-motion stations on the Hanford facility with adequate signal to noise 

ratios; and that had focal mechanisms determined by the PNSN. Only three earthquakes matched 

those criteria (see table 1). A fourth earthquake somewhat outside the 3D model near Walla 

Walla, Washington, which was recorded by station LTH in the basin, also was identified, but 

pre-dated the strong-motion instruments. 

The Prosser earthquake of October 18, 2008, had a duration magnitude of 3.7 and a depth 

of about 20 km (PNSN). This earthquake is especially important because it was near the downdip 

portion of the Rattlesnake Hills fault. The deep depth of the earthquake also provides a good test 

on whether the southern edge of the supra-basalt sediments (northern edge of Rattlesnake Hills)

generates surface waves from incident, upward-propagating S-waves.

We used the focal mechanism for this earthquake derived from Robert Herrmann using 

inversion of regional surface waves. This focal mechanism differed substantially from that 

determined by P-wave first motions by the PNSN. The moment magnitude, Mw, determined by 

Herrmann was 3.27. We applied the relocated hypocenter listed in the Rodriguez-Marek 

database, which had a depth of 18 km. 

A comparison of the velocity waveforms at station H4A for the data and the 3D 

simulation using the RDM is shown in figure 5 (left). Because of the high frequencies involved, 

we were not expecting to match the data “wiggle for wiggle.” However, the synthetics for the 

north-south (NS) component match the observed S-wave arrival time, peak velocities, and long 

duration of substantial ground motions after the initial S-wave arrival. For the east-west (EW)

and vertical components, the 3D synthetics capture some of the arrivals after the S-wave, but not 

all of them.

Two 1D, flat-layered velocity models were constructed based on the velocities at station 

H2E near the waste-treatment plant (fig. 6). This model has a thickness of supra-basalt sediments 

of 100 m. The 1D models used the thickness of sub-basalt sediments from the GDM and RDM. 

For these 1D simulations for modeling the recorded earthquakes, we did not use a velocity 

gradient for the sub-basalt sediments, applying depth-averaged average velocities and density. 

This did not significantly affect the waveforms of the 1D synthetics. The 3D finite-difference 

code was run using these flat-layered models. 

The synthetics for the 1D model clearly do not contain the large arrivals after the S-wave 

that are observed in the data in all three components (fig. 5). This applies to the 1D models using 

GDM and (RDM values for the top of basement. The 1D models produce very small surface 

waves, relative to the S-waves, because the epicentral distance is only slightly greater than the 

depth of the hypocenter. The main difference in the synthetics from the GDM and RDM 1D 

models is the predicted S-wave arrival time. However, this discrepancy between the observed 

and predicted arrival times from the 1D model is just the result of the difference of thickness of 

the sub-basalt sediments between stations H4A and H2E where the 1D models were derived.
5 

The waveforms from the 3D models using the different depths to basement are very 

similar. The observed waveforms at station H4A and the synthetics from the 3D models are 

shown in figure 7. The synthetics from both 3D models reproduce the substantial arrivals for 10 s 

after the S-waves that are observed in the data. For subsequent figures in this report, we used the 

3D model with the RDM values for the depth to top of basement, unless otherwise indicated. 

The arrivals observed after the S-wave are likely to be surface waves created near the 

southern edge of the supra-basalt sediments from the conversion of incident S-waves. There are 

several pieces of evidence for this interpretation: (1) these arrivals show dispersion in the 

observed acceleration traces with the long-period energy arriving first in the wave train, (2) their 

arrival 5 s after the S-wave is consistent with their creation at the southern edge of the basin and 

the 500 m/s group velocity expected for the Vs of the supra-basalt sediments, and (3) the 1D 

models show that the direct surface waves have very small amplitudes for this distance and 

source depth. These high-frequency surface waves are largely trapped in the supra-basalt 

sediments, even though the sediments are less than 100 m thick. The observed recordings at 

station LTH exhibits similar large arrivals after the initial S-wave, which are matched by the 

synthetics from the 3D model but not by those from the 1D models (fig. 8).

The observed recordings at station H2E show ringing with long durations of tens of 

seconds after the initial S-wave (fig. 9). The synthetics from the 3D model (RDM) match the 

peak amplitudes and capture some of the long-duration character of the observations. In contrast, 

the synthetics from the 1D model are smaller than the observed recordings and have virtually no 

amplitude after the initial S-wave. On the EW component, the 1D synthetics show a very small 

S-wave arrival, in contrast with the larger S-wave observed in the data and the long duration of 

motion. The 3D and 1D synthetics are similar at station H3A (fig. 10), because of the relatively 

thin supra-basalt sediments at this location. Both sets of synthetics do not contain the significant 

arrival observed about 2 s after the S-wave.

Station HAWA is located on rock and the data displays a simple S-wave for this station

in contrast to the more complex, long-duration arrivals for stations located on the supra-basalt 

sediments. These simple waveforms at station HAWA are matched by the 3D and 1D models

(fig. 11). A phase on the radial component that arrives between the P- and S-waves is shown in 

the seismograms in figure 12. The polarization and arrival time indicate that this arrival is an S to 

P conversion (SP phase) at the top of the basement. Synthetics using the GDM predict an SP 

arrival later than the observed arrival, indicating that the top of basement is deeper than in the 

GDM (fig. 12). The RDM produces a SP phase closer to the timing of the data, but the S-wave 

arrival is delayed relative to the data. This may indicate that the S-wave velocities in the RDM 

are too slow for the path from the hypocenter to station HAWA. 

The synthetics from the 3D model have spectral accelerations closer to the observed 

accelerations, compared to the 1D synthetics. The bias and standard deviation of the natural 

logarithms of the response spectral accelerations (5 percent critical damping; pseudo-spectral 

accelerations) at 0.5 and 1.0 s for the observations relative to the synthetics are shown in figure 

13. Here we considered the stations on the supra-basalt sediments: H2E, H3A, H4A, and LTH. 

The geometrical average of the spectral acceleration from the two horizontal components was 

used. The bias of the 3D synthetics is much lower than the bias from the 1D synthetics, 

indicating that the 3D synthetics better match the observed spectral accelerations than the 1D 
6 

synthetics. The 1D synthetics substantially under predict the observed spectral accelerations at 

0.5 and 1.0 s (positive bias indicates that the observed data have larger spectral accelerations 

than the synthetics). The standard deviations are somewhat smaller for the 3D synthetics relative 

to the 1D synthetics, illustrating that the 3D model produces a better fit to the observed 

variations of spectral accelerations between stations than the 1D model. 

The focal mechanism makes a significant difference in the fits to the response spectra. 

The bias and standard deviation between the 3D synthetics made with the Herrmann focal 

mechanism and those from the PNSN focal mechanism, both using the same 3D model are 

shown in figure 14. For 1 s period, the Hermann focal mechanism produces a better fit to the 

observed variations of spectral accelerations between the stations, as shown by the lower 

standard deviation. 

There is no significant difference in the fits to the data for the GDM and RDM of the 

depth to basement (3D models; fig. 15) for 0.4 and 1.0 s. This reflects the idea that it is the suprabasalt sediments that form the waveguide for the surface waves and the deeper structure is 

relatively inconsequential.

The energy was calculated from summing the square of the ground velocity as a function 

of time. The cumulative energy of the observed data and synthetics at station H4A is shown in 

figure 16. For the NS component, the energy of the observed data and the 3D synthetics are 

similar. The EW synthetics in the 3D model underestimate the observed energy on that 

component. The synthetics from the 1D model severely under predict the observed energy for 

both components. This is not unexpected given the short duration of the 1D synthetics compared 

to the data.

The 3D synthetics better predict the total energy of the observed records than the 1D 

synthetics. The bias and standard deviation of the total energy for observed and synthetic 

waveforms at stations H3A, H4A, H2E, and LTH are shown in figure 17. The bias is much 

lower for the 3D synthetics than the 1D synthetics, again demonstrating that the 3D model better 

fits the characteristic of the observed recordings than the 1D model and that the 1D model under 

predicts the observed energy.

Wooded Island Swarm Events

Two earthquakes in the Wooded Island area were simulated using the 3D and 1D models. 

These earthquakes were very shallow (depths less than 2 km). For both earthquakes, the 3D and 

1D simulations used an Mw of Md-0.4, based on the Mw to Md difference found for the Prosser 

earthquake. Herrmann did not determine focal mechanisms or moments for these earthquakes, so 

the PNSN first-motion focal mechanisms were applied. 

The observed and synthetic waveforms at stations H4A, H2E, and H2W for the 

September 2011 earthquake are shown in figure 18. At station H4A closer to the hypocenter, the 

3D and 1D synthetics underestimate the observed duration after the S-wave arrival. The 3D 

synthetics for stations H2E and H2W generally reproduce the long duration of the observed 

records, whereas the 1D synthetics do not (fig. 18). However, the 3D model overestimates the 

peak velocities at stations H2E and H2W. The 1D model also overestimates the peak velocities at 

these stations. The bias and standard deviation between the 0.5 and 1.0-s spectral accelerations 

from the synthetic and the observed data are similar between the 3D and 1D models (fig. 19). For 

this comparison, we used stations H3A, H4A, H2E, H2W, and H1K. Because of the shallow 
7 

depths of these earthquakes, the waveforms may be affected by scattering from topography that 

is not contained in the 3D and 1D models. This scattering by topography may be especially 

important for stations near the hypocenter, which is located approximately beneath the bluffs 

along the Columbia River. 

The results are similar for the October 2011 earthquake. The 3D and 1D models 

underestimate the ringing observed after the S-wave at station H4A, close to the epicenter (fig. 

20). For stations H2E and H2W, the 3D synthetics capture the long duration of shaking of the 

observed recordings. The 1D synthetics have shorter durations of shaking than the observed 

recordings. The bias of the observed and synthetic accelerations at 0.5 and 1.0 s are similar for 

the 3D and 1D models (fig. 21). The standard deviation at 1.0 s is lower for the 1D model than 

the 3D model, indicating that the 1D model better fits the observed variations of spectral 

accelerations at 1.0 s for this earthquake.

Walla Walla Earthquake

The earthquake near Walla Walla, Washington, is somewhat outside of the area of the 3D 

model and was not recorded on the strong-motion instruments. However, the earthquake was 

recorded at station LTH. The 3D model was extended to the east and to the south by continuing 

the velocity model at the southern and eastern edges. Synthetics from the 3D model show that it 

contains substantial arrivals after the S-wave at station LTH (NS component), although these do 

not have as much duration as the arrivals in the observed data (fig. 22). The synthetics for the 1D 

model are much shorter duration than the observed data at station LTH.

Simulations for Magnitude 6.5–6.8 Earthquakes

Given the importance of 3D propagation effects in the observed data, it is critical to 

assess how the 3D model affects the spectral accelerations predicted for M6 earthquakes from 

nearby faults, compared to those from 1D models. Here, we considered three sources: M6.8 

earthquakes on the western part of the Rattlesnake Hills fault, M6.7 earthquakes on the eastern 

part of the Rattlesnake Hills fault, and M6.6 earthquakes on the Gable Mountain fault. The 

surface projections of the rupture zones are shown in figure 23. The rupture zones for the 

Rattlesnake Hills fault dipped to the south at 55 degrees. Ruptures on the Gable Mountain fault 

dipped to the south at 70 degrees. In addition, a dipping plane source was studied to approximate 

the regional S-waves from an M8–9 earthquake on the Cascadia subduction zone. Synthetics 

were calculated at the five sites used for the PSHA (fig. 23): A (200 East area and adjacent to the 

Waste Treatment Plant), B (200 West area), C (Columbia Generating Station), D (100 BC area), 

and E (300 area).

The general methodology for generating the rupture history for the M6 earthquakes 

follows that developed by Frankel (2009) for the deterministic portion of the calculation. First, 

slip distributions were developed, characterized by random fields with a wavenumber (k) 

spectrum proportional to k

-2 and a Gaussian probability distribution. The along-strike and downdip correlation distances for each magnitude were taken from the formulas of Mai and Beroza 

(2002). A 500 m sub-fault dimension was used, with the rupture plane extending from depths of 

3 to 15 km. One of the slip distributions used for the Rattlesnake Hills West rupture is shown in 
8 

figure 24. This slip for each sub-fault was then converted into seismic moment, using the depthdependent Vs and density along the rupture plane from the 3D velocity model. The moment 

magnitudes of the earthquakes were determined from the fault lengths of the three rupture 

planes, using the surface rupture length-magnitude formula in Wells and Coppersmith (1994) as 

a guide. 

The rise times of slip for each sub-fault were determined by specifying an average slip 

velocity following the procedure of Frankel (2009), which compared the spectral acceleration 

from synthetics with those predicted by the Next Generation of Attenuation (NGA) relations. 

Frankel (2009) determined that an average slip velocity of 2.7 m/s reproduced the 1.0 and 3.0 s 

spectral accelerations of the NGA relations for earthquakes with M5.5, 6.5, and 7.5. The 

dynamic stress drop is proportional to the slip velocity divided by the rupture velocity (Vr). 

Because the Vr is specified as a fraction of the local shear-wave velocity, the average slip 

velocity was multiplied by the sub-fault Vs divided by3,500 m/s, to approximately maintain a 

constant dynamic stress drop as the Vs varies with the location of the sub-fault. The slip velocity 

was randomized with a uniform distribution with ±20 percent variation. The rise time for each 

sub-fault was determined from the ratio of the final slip divided by the slip velocity. A Brune 

(1970) pulse was used for the source time function for each sub-fault.

One key set of parameters is the rupture initiation time for each sub-fault. The fraction of 

Vr/Vs for each sub-fault was varied between 0.62 and 0.98. The variation in Vr was specified as 

20 percent standard deviation around 0.8 of the local Vs. The perturbation in the rupture velocity 

was proportional to the final slip for that sub-fault, such that areas with higher slip have faster 

rupture velocity. This follows results of dynamic rupture simulations. The Vs varied with 

position along each rupture plane according to the location of the rupture in the 3D model. Twodimensional ray tracing along the rupture plane was used to calculate the initiation time of each

sub-fault. Separate ray-tracing runs were done to calculate the initiation time of each sub-fault.

The local rupture velocity for each ray-tracing run was calculated from the Vr/Vs fraction 

determined for that given sub-fault. In this manner, areas with high slip will rupture sooner than 

adjacent areas with low slip. We also added a small random component to the rupture time of 

each sub-fault, similar to that done in Frankel (2009). The rupture initiation times for one of the 

simulations of the Rattlesnake Hills West rupture is shown in figure 24. There is quite a bit of 

randomness to the initiation times, similar to that reported in dynamic rupture simulations with

complex pre-stress and strength. The rake of each sub-fault is derived from a uniform random 

distribution with a variation of ±20 degrees around the average rake of 120 degrees. 

One-Dimensional Models and Reflectivity Method

Ten rupture scenarios were done for each fault source, resulting in 30 3D simulations. For 

each scenario, random draws were taken for the along-strike and down-dip locations of the 

hypocenter and for choosing among five slip distributions generated with different random 

number seeds. The down-dip location of the hypocenter was constrained to be within the lower 

one-half of the rupture. This has been the case for all M6.5 or larger oblique or thrust 

earthquakes in California since 1970 that have well-resolved hypocentral and aftershock depths. 

The 3D and 1D simulations used the same source files, specifying the seismic moment, 

start time, rise time, and focal mechanism for each sub-fault.
9 

Green’s function for the 1D model for each station were calculated using the frequencywavenumber integration code of Zhu and Rivera (2002). The 1D model for each of the five 

stations was developed by Paul Thorne. We decided to follow an approach suggested by the 

Technical Integration Team for the SSHAC study. Beneath the Grand Ronde basalts, we used the 

same velocity profile for all five stations. The 1D velocity model for depths greater than 3 km 

corresponded to the Vs profile at the center (measured along strike) of the rupture plane for each 

fault, as determined from the 3D model. This ensures that the Vp, Vs, and density values with 

depth for the sub-faults are approximately the same between the 3D and 1D models. For the 

Rattlesnake Hills West and Gable Mountain ruptures, it was found that the Vs profiles between 3 

and 15 km were equivalent in the 3D model. Therefore, the same 1D velocity models were used 

for these faults at depths greater than 3 km. The Vs profiles for the five sites used in the 1D 

simulations for the Rattlesnake Hills West and Gable Mountain sources are shown in figure 25. 

Separate sets of Green’s functions were calculated for the five sites for the Rattlesnake Hills East 

and Cascadia sources. The depth to basement for the Rattlesnake Hills East source was about 820 

m less than that for the Rattlesnake Hills West and Gable Mountain sources. The depth to 

basement for the Cascadia source was 340 m greater than the Rattlesnake Hills West and Gable 

Mountain sources. There also were some differences in the depth to top of pre-Miocene 

sediments between these profiles that were incorporated in the velocity profiles. All together, 

there were 15 sets of 1D Green’s functions calculated (5 sites and 3 different profiles for depths 

greater than 3 km). There were150 1D simulations completed: 3 rupture zones, 10 scenarios for 

each rupture zone, and 5 Vs profiles for each scenario. 

Because of the large number of 1D simulations requested by the Technical Integration

Team, the 1D simulations could not be completed using the 3D finite-difference code in the 

allotted time. We compared the waveforms and spectral accelerations between the synthetics 

from the reflectivity summation code and the 3D code using the same 1D velocity profile and 

source file. We determined that the spectral accelerations were within 20 percent at 0.5 and 0.75 

s and within 10 percent for periods greater than or equal to 1.0 s. The velocity waveforms at any 

given station were almost identical. The ratios of spectral accelerations shown below are for the 

pseudo spectral accelerations, using the geometrical average of the two horizontal components. 

Rattlesnake Hills Fault, West, Moment Magnitude 6.8

Synthetics from the 3D model at sites A (200 East area) and B (200 West area) for the 

Rattlesnake Hills West source show substantial later arrivals after the initial S-wave, arrivals that 

are not present in the synthetics from the 1D models with the same source (fig. 26). The 1D 

synthetics shown for each site were simulated with the 1D model for that site. Sites A and B 

have the largest thickness of supra-basalt sediments, 94 and 39 m, respectively (fig. 1). At site A, 

the peak velocity in the 3D synthetics occurs 10–20 s after the initial S-wave, in contrast to the 

short-duration 1D synthetics (fig. 26).These later arrivals in the 3D synthetics likely are surface 

waves produced by conversion of incident S-waves at the southern edge of the supra-basalt 

sediments. These surface waves also may be reflected and amplified by other features in the 

shallow portion of the 3D velocity model. For the other sites (C, D, and E), the 3D and 1D 

synthetics are similar, although there is short-period energy superimposed on the long-period 

waveforms in the 3D synthetics that is not present in the 1D synthetics (fig. 26). This indicates 
10

that the reverberations of S-waves in the supra-basalt sediments may be stronger in the 3D model 

than the 1D model. This could be caused by the non-horizontal interface in the 3D model at the 

base of the supra-basalt sediments. The thickness of supra-basalt sediments in the 1D models for 

sites C, D, and E are 14, 25, and 17 m, respectively.

Most sites show higher spectral accelerations at periods of 0.5, 0.75, and 1.0 s for the 3D 

model than the 1D models. The natural logarithm (ln) of the ratios of the spectral accelerations 

between the 3D and 1D synthetics are shown in figure 27. This figure displays the 3D/1D ratios 

at each site for the 10 simulations for that rupture zone, as well as the mean and standard 

deviation of the ln ratio at each period. The spectral accelerations for periods of 3–10 s generally 

are similar between the 3D and 1D synthetics for all the sites. There is substantial variation in the 

ratios at each site between the 10 source scenarios, although the general dependence on period is 

similar. The standard deviations of the ratios generally decrease with increasing period. At 

periods of less than or equal to 2 s, the standard deviations of the ln ratios for the 10 simulations 

can be as large as about 0.3 ln units. 

The ratios between the 3D and 1D models are largest for sites A and B with the thickest 

supra-basalt sediments. For site A, the ratios start to increase at 2 s period, reach a peak oln ratio 

of 1.5 (factor of 4.5) at 0.75 s, and have a ln ratio of 1.1 (factor of 3.0) at 0.5 s. At site B, the 

ratios start to increase at 1.5 s, reaching ln ratio of 1.1 (factor of 3.0) at 0.5 s. The 3D/1D ratios at 

site C start to increase at 1 s and have a maximum ln ratio of 0.75 (factor of 2.1) at 0.5 s. The 

ratios at site D start to increase at 2 s and peak at 0.75 s at a ln ratio of 1.0. For site D, the ln of 

the ratios are slightly negative for 3–5 s, indicating that the spectral accelerations for the 3D 

synthetics are slightly smaller than those for the 1D synthetics at those periods. The ln of the 

3D/1D ratios at site E are smaller, with no values greater than 0.5 (less than a factor of 1.6). 

The ratios generally increase with decreasing period, starting at some transition period. 

This transition period is greater for sites A and B with the thickest supra-basalt sediments than 

those for the other sites. 

We were initially concerned that the 3D/1D ratios for site C could be affected by the 

discretization of the 3D model. For site C, the thickness of the supra-basalt sediments is 14 m, 

less than the vertical grid spacing (25 m) used in the 3D model. The thickness of the sediments in 

the 3D model is set to one grid spacing for this site (fig. 1). For vertically propagating S-waves, 

the resonance frequency is about 5 Hz for a 25 m thick layer with Vs of 425 m/s. So the 

resonance frequency is lower in the 3D model than in the 1D models, where the supra-basalt 

layers are not discretized. We were concerned that this may cause the 3D synthetics to have 

higher spectral accelerations at 0.5 and 0.75 s than the 1D synthetics. However, we tested this 

idea by making 1D synthetics at site C using Green’s functions calculated for site D, which has a 

thickness of 25 m. We found no significant difference in the spectral accelerations at 0.5, 0.75, 

and 1.0 s at site C between the 1D models with thicknesses of 25 and 14 m, indicating that the

differences between the 3D and 1D values are not driven by the variation in sediment thickness 

between the 3D and 1D models. The observation that the 3D/1D ratios for site E, which has a 

sediment thickness of 17 m, are close to 1 also indicates that the discretization of the 3D model is 

not the cause of the higher ratios at site C.
11

Rattlesnake Hills Fault, East, Moment Magnitude 6.7

The synthetics from the 3D model for the Rattlesnake Hills East source also display 

prominent later arrivals after the S-waves at sites A and B, indicating basin surface waves that 

are not present in the 1D synthetics (fig. 28). The peak arrivals in the 3D synthetics for these 

sites occur well after the initial S-wave, indicating that these late arrivals are basin surface waves 

not present in the 1D models. The 3D and 1D synthetics at sites C, D, and E are similar, although 

there is more high-frequency energy apparent in the 3D synthetics. The synthetics at site E 

contain a strong S-wave pulse from rupture directivity. 

The ratios of the 3D/1D spectral accelerations show that the 3D model produces higher 

spectral accelerations at 0.5–1.0 s at all the sites, with the largest ratios at sites A and B (fig. 29).

For site A, the ratios start to increase at 2 s and peak at 0.75 s at a ln ratio of 1.3 (factor of 3.7). 

The ln ratios at site B are slightly negative at 3 and 4 s. The 3D/1D ratios start to increase at 2 s 

and reach a ln ratio of 1.0 at 0.75 and 0.5 s. The ratios at site C are positive for 0.5 and 0.75 s, 

and have a peak 1n ratio of about 0.75 (factor of 0.75). The 3D/1D ratios are small at site D, 

reaching ln ratio of 0.5 at 0.5 and 0.75 s. For site E, the ln ratio increases to about 0.9 (factor of 

2.5) at 0.5 s. 

Gable Mountain Fault, Moment Magnitude 6.6

This source is different from the Rattlesnake Hills sources in that most of the PSHA sites 

are directly or nearly above the rupture. S-waves traveling from the source to these sites are 

upward propagating and do not interact with the edge of the area with supra-basalt sediments, 

except perhaps by reflection. However, the 3D structure could potentially cause focusing of the 

upward-traveling S-waves.

1D and 3D synthetics for one of the Gable Mountain sources are shown in figure 30. For 

sites A,B,C, and D, the 3D and 1D synthetics for each site look similar and contain strong Swave pulses produced by forward rupture directivity. There appears to be more high-frequency 

energy in the 3D synthetics than in the 1Dsynthetics. The 3D synthetics at site E, the farthest 

from the source, have prominent later arrivals that are not present in the 1D synthetics. This is 

especially notable on the vertical component. The 3D model is apparently producing more basin 

surface waves at this relatively distant site than the 1D model.

For most of the sites, the 3D/1D ratios of spectral accelerations (fig. 31) are smaller than 

those determined for the Rattlesnake Hills sources, but still can be significant in some cases. The 

3D/1D ratio for the spectral accelerations was largest for site E, with values reaching a ln ratio of 

1.0 (factor of 2.8) at 0.5 s. The 3D/1D amplifications were small for site A and actually negative 

at 0.5 s. The ratios at sites B, C, and D reach a ln ratio of 0.6 (factor of 1.8) at 0.5 s. 

Cascadia

To approximate an incoming S-wave from an earthquake on the Cascadia subduction 

zone, we used a set of sources arranged as a dipping plane beneath the western portion of the 3D 

velocity model. The planar source consisted of sources at 500 m spacing with a north-south 

strike and a dip of 41 degrees, down to the east. This dip was selected because it represents the 

angle of incidence of a critical reflection from the Moho, for a Vs of 3,500 m/s in the crust and 
12

4,600 m/s in the top of the mantle. The sources at all of these grid points are initiated at the same 

time. A focal mechanism with a strike of 0 degrees, dip of 41 degrees, and rake of 90 degrees

was applied to each source. A rise time of 6 s was used for each source point, approximating that 

expected for a great Cascadia earthquake. This rise time is not important to the spectral 

acceleration ratios. The same source file was used for the 3D and 1D simulations. 

The 3D and 1D synthetics generally display a simple S-wave pulse (fig. 32) with shortperiod energy superimposed. The 3D synthetics have stronger, short-period reverberations after 

the initial arrival than the 1D synthetics. Site E, farthest to the east, shows a long-period arrival 

about 10 s after the initial S-wave on the vertical and east-west (radial) components of the 3D 

and 1D synthetics. 

There is substantial amplification of the spectral accelerations from the 3D model with 

respect to those from the 1D model. The ratios of the spectral accelerations for the 3D and 1D

simulations for the five PSHA sites are shown in figure 33. Again, the 1D simulations used 

Green’s functions calculated for the 1D model for each site. Amplification of the synthetics from 

the 3D simulation relative to the 1D simulations increased from periods of 1.0–0.5 s. The ratio 

for site A peaks at 0.75 s with ln ratio of 2.2 (factor of 9) and was ln ratio of 1.5 at 0.5 s and ln 

ratio of 1.3 at 1.0 s. For site C, the 3D/1D ratio was ln ratio of 0.5 (factor of 1.6) at 4 s and 

increased to ln ratio of 1.4 (factor of 4.1) at 0.75 and 0.5 s. For site D, the 3D/1D ln ratio was1.4 

(factor of 4.1) at 0.5 s. Site E shows the lowest amplification, reaching ln ratio of 0.5 (factor of 

1.6) at 0.75 and 0.5 s. 

It should be emphasized that the dipping-plane wave source for Cascadia is a pragmatic 

simplification. Obviously, the input regional S-wave will have a more complicated configuration 

because of the structure between the subduction zone and the western portion of the 3D model 

used here. This configuration also cannot assess the amplification of regional surface waves from 

the Cascadia source.

Conclusions

The three-dimensional (3D) model produces synthetics that are a better fit to the observed 

waveforms, durations, spectral accelerations, and energy for the Prosser earthquake than the 

synthetics from the one-dimensional (1D) model. The 3D model reproduces the observed basin 

surface waves from the Prosser earthquake, surface waves that were likely created along the 

southern edge of the supra-basalt sediments and whose large amplitudes cannot be explained 

with the 1D model. The synthetics from the 3D model for the Wooded Island and Walla Walla, 

Washington, earthquakes show durations of shaking more similar to the observed records, 

whereas the 1D synthetics do not display the later arrivals seen in the data. 

For the simulations of magnitude 6 earthquakes on the Rattlesnake Hills fault, the 3D 

model produced synthetics with substantially higher spectral accelerations at periods from 0.5 to 

1.0 seconds (s) than the 1D models. The durations of the 3D synthetics at sites A and B are much 

longer than those from the 1D models. The ratios between the 3D and 1D spectral accelerations 

reach a factor 4.5 and generally are larger for the sites with the thickest supra-basalt sediments 

(sites A and B). The 3D/1D ratios typically are more subdued for the Gable Mountain source, 

although they still can be significant. For the Cascadia plane-wave source, the amplification of 

0.2 and 1.0 s spectral accelerations of the 3D model relative to the 1D models is substantial (as 

much as a factor of 9 for one site) at most of the sites and significant amplification extended to 

periods as long as 4 s period for one site. 
13

This study has demonstrated that simulations in a 3D velocity model for the Hanford area 

are required to accurately predict the spectral acceleration, energy, and duration of strong ground 

motions from moderate and large earthquakes on the Rattlesnake Hills and Gable Mountain 

faults, as well as from great earthquakes on the Cascadia subduction zone. These 3D effects need 

to be considered to properly assess the probabilistic seismic hazard for structures and storage 

tanks at the Hanford facility and the Columbia Nuclear Generating Station. The deployment of 

dense arrays of surface seismometers and borehole arrays on the Hanford facility would help to 

better quantify and predict strong ground motions from future earthquakes.

Acknowledgments

The U.S. Department of Energy provided funding for this study, as part of the Senior 

Seismic Hazard Analysis Committee probabilistic seismic hazard assessment for the Hanford 

facility. We thank Robert Bryce for his support of this project, for his assistance with 

coordinating this effort, and for his comments on the report. Stephen Reidel helped produce the 

3D velocity model and his insights on the constraints for the model were illuminating. We thank 

Pengcheng Liu for providing his 3D finite-difference program. Larry Baker helped with 

optimizing the program. We thank Julian Bommer, Robert Youngs, Gabriel Toro, Ken 

Campbell, Adrian Rodriguez-Marek, Linda Al Atik, and Walter Silva for their comments and 

suggestions. Earthquake seismograms, locations, focal mechanisms, and magnitudes provided by 

the Pacific Northwest Seismic Network were essential to this study. We thank Craig Weaver for 

his support and advice. Paul Bodin and John Vidale reviewed the report and provided helpful 

comments. 

References Cited

Alford, R.M., Kelley, K.R., and Boore, D.M., 1974, Accuracy of finite-difference modeling of 

the acoustic wave equation: Geophysics, v. 39, p. 834–842. 

Brune, J.N., 1970, Tectonic stress and the spectra of seismic shear waves from earthquakes: 

Journal of Geophysical Research, v. 76, p. 4,997–5,002. 

Burns, E.R., Morgan, D.S., Peavler, R.S., and Kahle, S.C., 2011, Three-dimensional model of the 

geologic framework for the Columbia Plateau Regional Aquifer System, Idaho, Oregon, and 

Washington: U.S. Geological Survey Scientific Investigations Report 2010-5246, 44 p., 

http://pubs.usgs.gov/sir/2010/5246/. 

Frankel, A.D., 2009, A constant stress-drop model for producing broadband synthetic 

seismograms—comparison with the Next Generation Attenuation relations: Bulletin of the 

Seismological Society of America, v. 99, p. 664–680. 

Frankel, A.D., Stephenson, W.J., Carver, D.L., Williams, R.A., Odum, J.K., and Rhea, S., 2007, 

Seismic hazard maps for Seattle incorporating 3D sedimentary basin effects, nonlinear site 

response, and rupture directivity: U.S. Geological Survey Open-File Report, 2007–1175, 70 p. 
14

Frankel, A., Stephenson, W., and Carver, D., 2009, Sedimentary basin effects in Seattle, 

Washington: ground-motion observations and 3D simulation: Bulletin of the Seismological 

Society of America, v. 99, p. 1579-1611. 

Geomatrix, 1996, Probabilistic seismic hazard analysis—DOE Hanford site, Washington: San

Francisco, California, Geomatrix Consultants, Inc., WHC-SD-W236A-TI-002, rev. 1, 

Glover, D.W., 1985, Crustal structure of the Columbia Basin, Washington, from borehole and 

refraction data: University of Washington, Geophysics Program, MS Thesis. 

Liu, P-C, and Archuleta, R.J., 2002, The effect of a low-velocity surface layer on simulated 

ground motion: Seismological Research Letters, v. 73, p. 267. 

Mai, P.M., and Beroza, G.C., 2002, A spatial random field model to characterize complexity in 

earthquake slip: Journal of Geophysical Research, v. 107, p. 1–21. 

Wells, D.L., and Coppersmith, K.J., 1994, New empirical relationships among magnitude, 

rupture length, rupture width, and surface displacements: Bulletin of the Seismological Society 

of America, v. 84, p. 974–1,002. 

Zhu, L., and Rivera, L., 2002, Computation of dynamic and static displacement from a point 

source in multi-layered media: Geophysical Journal International, v. 148, p. 619–627. 
15

Figure 1. Map of depth (meters) to S-wave velocity (Vs) of 1,500 meters per second for the Hanford area, 

Washington. This equals the thickness of the supra-basalt sediments in the three-dimensional (3D) model. 

Contour interval is 25 meters, the same as the vertical grid spacing in the 3D model. Black circles are 

seismic stations (labeled). Red circles are locations of sites used in the probabilistic seismic hazard 

assessment. 
16

Figure 2. Cut-out view of three-dimensional (3D) velocity model that used depth to basement derived from 

gravity-derived model (GDM). Blue is basement, yellow designates sub-basalt sediments, gray represents 

basalts, and pink and purple denote supra-basalt sediments. 
17

Figure 3. 3D views of the elevation of the top of basement (m) for the gravity-derived model (GDM; top) 

and the refraction-derived model (RDM; bottom).
18

Figure 4. Satellite image of Hanford area from GoogleEarthTM showing locations of strong-motion 

accelerometers (circles), broadband stations (squares), and earthquakes studied (stars). Red star indicates 

initial location of Prosser earthquake from Pacific Northwest Seismic Network (PNSN); white star 

represents re-location from database of A. Rodriguez-Marek. Other earthquake epicenters from PNSN. The 

two Wooded Island earthquakes have very similar epicenters.
19

Figure 5. Observed (black) and synthetic (red) velocity waveforms at station H4A for the Prosser 

earthquake. Vertical scale is in meters per second. Left panel shows synthetics from 3D velocity model. 

Center and right panels show synthetics using 1D models in figure 5. The synthetics from the 3D model 

show substantial arrivals after the S-wave similar to the observed waveforms. These later arrivals are 

surface waves (labeled SW) created by the 3D structure. The synthetics from the 1D models do not 

reproduce the observed arrivals after the S-wave. GDM, gravity-derived model; RDM, refraction-derived 

model; NS, north-south; EW, east-west; UD, up-down; OT, origin time.
20

Figure 6. One-dimensional (1D) velocity profiles at station H2E used in the modeling of the recorded 

earthquakes. The gravity-derived model (GDM) uses the depth to top of basement derived from gravity 

modeling. The refraction-derived model (RDM) applies the depth to top of basement constrained by seismic 

refraction studies.
21

Figure 7. Observed (black) and synthetic (red) velocity waveforms at station H4A for the Prosser 

earthquake. The synthetics were made with two three-dimensional (3D) crustal models (refraction-derived 

model [RDM] and gravity-derived model [GDM]), with different depths to top of basement. 
22

Figure 8. Observed (black) and synthetic (red) velocity waveforms at station LTH for the Prosser 

earthquake. Note that the synthetics from the three-dimensional (3D) model capture much of the energy 

observed after the S-wave, whereas the one-dimensional gravity-derived model (1D GDM) synthetics do 

not.
23

Figure 9. Observed (black) and synthetic (red) velocity waveforms at station H2E for the Prosser 

earthquake. Left panel shows synthetics for the 3D model. Right panel shows synthetics for the 1D model.
24

Figure 10. Observed (black) and synthetic (red) velocity waveforms at station H3A for the Prosser 

earthquake. Left panel shows synthetics for the 3D model. Right panel shows synthetics for the 1D model.
25

Figure 11. Observed (black) and synthetic (red) velocity waveforms at station HAWA for the Prosser 

earthquake. Left panel shows synthetics for the 3D model. Right panel shows synthetics for the 1D model.. 
26

Figure 12. Observed and synthetics velocity waveforms (NS component) at station HAWA for the Prosser 

earthquake. Note arrival between P- and S-wave is inferred to be an S to P conversion (SP) at the top of 

basement, based on travel time calculations. 
27

Figure 13. Bias and standard deviation of the spectral accelerations from the observed recordings at 0.5 

and 1.0 s with respect to the spectral accelerations from the synthetics for the Prosser earthquakes 

(stations H3A, H4A, H2E, and LTH). Note that the bias is larger for the synthetics from the one-dimensional 

(1D) model than those from the three-dimensional (3D) model (refraction-derived model).
28

Figure 14. Bias and standard deviation of observed versus synthetic spectral accelerations for the Prosser 

earthquake. Different focal mechanisms and slightly different hypocenters were considered in threedimensional (3D) simulations for the same 3D velocity model (refraction-derived model). 
29

Figure 15. Bias and standard deviation of observed versus synthetic spectral accelerations for the Prosser 

earthquake using three-dimensional (3D) models with different depths to basement.
30

Figure 16. Cumulative energy at station H4A for Prosser earthquake plotted with respect to time after origin 

time for observed seismograms (black and red lines) compared to synthetics (green and blue lines) from

three-dimensional (3D) (top) and one-dimensional (1D) (bottom) velocity models. Results are plotted for the 

east-west and north-south components.
31

Figure 17. Bias and standard deviation of energy derived from observed seismograms compared with 

those from the three-dimensional (3D) and one-dimensional (1D) synthetics, for the Prosser earthquake 

(stations H3A, H4A, H2E, and LTH).
32

Figure 18. Observed (black) and synthetic (red) velocity waveforms for the September 2011 Wooded 

Island earthquake. Left panel shows synthetics for the 3D model. Right panel shows synthetics for the 1D 

model.
33

Figure 19. Bias and standard deviation of the observed and synthetics spectral accelerations for the 

September 2011 Wooded Island earthquake. 
34

Figure 20. Observed (black) and synthetic (red) velocity waveforms for the October 2011 Wooded Island 

earthquake. Left panel is for the 3D model. Right panel is for the 1D model. NS, north-south; EW, eastwest; UD, up-down; OT,origin time. 
35

Figure 21. Bias and standard deviation of the observed and synthetic spectral accelerations for the 

October 2011 Wooded Island earthquake.
36

Figure 22. Observed (black) and synthetic (red) velocity waveforms at station LTH for the Walla Walla 

earthquake. Left panel shows synthetics for the 3D model. Right panel shows synthetics for the 1D model.
37

Figure 23. Map showing surficial S-wave velocity (Vs)(in meters per second) of three-dimensional model 

and surface projections (red rectangles) of rupture planes used in the simulations. Circles denote sites 

used in the probabilistic seismic hazard assessment and for the receivers in the simulations. A (200 East 

area and adjacent to the Waste Treatment Plant), B (200 West area), C (Columbia Generating Station), D 

(100 BC area), and E (300 area). 
38

Figure 24. Final slip and rupture initiation times for one of the simulations for the Rattlesnake Hills West 

source.
39

Figure 25. S-wave velocity (Vs) profiles used in the one-dimensional (1D) simulations for the Rattlesnake 

Hills West and Gable Mountain sources (bottom panel shows shallow portion). The thickness of the basalts 

and the depth to the top of basement vary somewhat for the Vs profiles used for the Rattlesnake Hills East 

and Cascadia 1D simulations (see text).
40

Figure 26. Velocity seismograms for one of the simulations of a magnitude 6.8 earthquake on the 

Rattlesnake Hills West source. Synthetics for three-dimensional model are in black, synthetics for onedimensional models are in red, and vertical scale is in meters per second. 
41

Figure 27. Natural logarithms of ratios of spectral accelerations between the three-dimensional (3D) and 

one-dimensional (1D) models for the magnitude 6.8 Rattlesnake Hills West source for the sites used for the 

probabilistic seismic hazard assessment. Each colored line represents the ratios for one rupture simulation. 

Positive values indicate that the spectral accelerations from the 3D model exceed those from the 1D model 

for that site. Dots with error bars are the mean and standard deviation, respectively, of the natural 

logarithms of the ratios for each period.
42

Figure 28. Velocity seismograms for one of the simulations of a M6.7 earthquake on the Rattlesnake Hills 

East source. Synthetics for three-dimensional model are in black, synthetics for one-dimensional models 

are in red, vertical scale is in meter per second. 
43

Figure 29. Natural logarithms of ratios of spectral accelerations between the three-dimensional and onedimensional models for the magnitude 6.7 Rattlesnake Hills East source for the sites used for the 

probabilistic seismic hazard assessment. 
44

Figure 30. Velocity seismograms for one of the simulations of a magnitude 6.6 earthquake on the Gable 

Mountain fault. Synthetics for three-dimensional model are in black., synthetics for one-dimensional models 

are in red, and vertical scale is in meters per second. 
45

Figure 31. Natural logarithms of ratios of spectral accelerations between the three-dimensional and onedimensional models for the magnitude 6.6 Gable Mountain source for the sites used for the probabilistic 

seismic hazard assessment. 
46

Figure 32. Velocity seismograms for simulation for Cascadia plane-wave source Synthetics for threedimensional model are in black, synthetics for one-dimensional models are in red, and vertical scale is in 

arbitrary units.
47

Figure 33. Natural logarithms of ratios of spectral accelerations from the three-dimensional model relative 

to those from the one-dimensional model at each site for the Cascadia plane wave source.
48

Table 1. Source parameters of earthquakes simulated in this study. 

[Source: Pacific Northwest Seismic Network (PNSN), unless otherwise indicated. Abbreviation: km, kilometer]

Yr Mo Day Hour 

Min

Latitude

(°N)

Longitude 

(°W)

Depth 

(km)

MagnitudeStrike Dip Rake

Prosser 

(PNSN focal 

mechanism 

and location)

2008 05 18 15:20 46.168 119.550 20.1 Md 3.7 60 45 80

Prosser 

(Herrmann 

focal 

mechanism 

and database 

location)

2008 05 18 15:20 46.156 119.528 18.0 Mw 3.27 60 55 35

Wooded 

Island

2011 09 04 04:13 46.411 119.260 1.8 Md 3.7 340 45 130

Wooded 

Island

2011 10 15 06:11 46.408 119.262 1.4 Md 3.4 75 40 21

Walla Walla 2006 12 20 09:43 46.095 118.513 13.6 Md 3.4 120 55 110