Probabilistic Estimates of Ground Motion in the Los Angeles Basin from Scenario Earthquakes on the San Andreas Fault

Abstract

Kinematic source inversions of past earthquakes in the magnitude range of 6–8 are used to simulate 60 scenario earthquakes on the San Andreas fault. The unilateral rupture scenario earthquakes are hypothetically located at 6 locations spread out uniformly along the southern section of the fault, each associated with two hypocenters and rupture directions. Probabilities of occurrence over the next 30 years are assigned to each of these earthquakes by mapping the probabilities of 10,445 plausible earthquakes postulated for this section of the fault by the Uniform California Earthquake Rupture Forecast. Three-component broadband ground motion histories are computed at 636 sites in the greater Los Angeles metropolitan area by superposing short-period (0.2–2.0 s) empirical Green’s function synthetics on top of long-period (>2.0 s) spectral element synthetics. The earthquake probabilities and the computed ground motions are combined to develop probabilistic estimates of ground shaking in the region from San Andreas fault earthquakes over the next 30 years. The results could be useful in city planning, emergency management, and building code enhancement.

1. Introduction

The interactions of the north American and the Pacific tectonic plates across much of California have created a network of major and minor active faults in the proximity of major cities such as Los Angeles and San Francisco that are capable of generating earthquakes as large as $M_w$ 8.3 or so. Major north-south trending faults in the vicinity of Los Angeles include the San Andreas fault, the San Jacinto fault, the Elsinore fault, and the Newport Inglewood fault. East-west trending faults include the Santa Monica-Hollywood-Raymond fault, the Sierra Madre fault, and the Puente Hills blind-thrust fault. Other yet-to-be discovered blind-thrust faults may be present as well. The proximity of these faults to the Los Angeles metropolitan area, the existence of large number of tall steel structures in the region, and the unexpected brittle failures in several steel buildings during the 1994 Northridge earthquake have prompted several investigations of the performance of these types of buildings (mainly of the steel moment frame variety) under hypothetical earthquake scenarios (e.g., [1,2,3,4,5,6,7]). These studies have focused on one or two events and their effects on structures. While such studies are useful to gain insights into potential outcomes, they are rarely suitable for rational decision-making as such one-off hypothetical events cannot capture the totality of the seismic hazard at the site. Here, we attempt to address this problem for earthquakes on the San Andreas fault.

We select kinematic finite source inversions of past earthquakes on geometrically similar faults and map these onto multiple locations on the southern San Andreas fault. We allow these earthquakes to propagate in two alternate directions, north-to-south and south-to-north. Here, we will refer to these simulated earthquakes as “scenario earthquakes”. We simulate a total of 60 scenario earthquakes spread over a magnitude range of 6–8. For each scenario earthquake, we compute 3-component ground motion histories at 636 “analysis” or “target” sites on a 3.5 km grid in southern California using the spectral element method for the low frequencies (<0.5 Hz) and empirical Green’s functions for the high frequencies ($0.5$–$5.0$ Hz).

Rupture forecasts such as the Uniform California Earthquake Rupture Forecast (UCERF, [8,9]) combine data from several sources, including local earthquake catalogs, magnitude frequency distributions (usually derived using global earthquake catalogs), paleoseismic observations, and GPS measurements of tectonic movement, and rules about rupture propagation between faults to predict probabilities of all plausible earthquakes on known faults during a specified time interval. UCERF bases these probabilities on four modeling components: (i) a (fault) model of the physical geometry of known California faults; (ii) a deformation model of slip rates and related factors for each fault section; (iii) an earthquake rate model of the region; and (iv) a probability model. It hypothesizes hundreds of thousands of ruptures (referred to as “forecast earthquakes” in this article) on specific seismogenic locations of faults and provides yearly occurrence rates that are most consistent with observations. These rates are transformed to probabilities of occurrence assuming an underlying probability distribution such as the Poisson distribution. Methods to prune this exhaustive set of ruptures down to those that contribute significantly to the hazard at a given site have been under development [10,11]. Here, we describe a rational method to redistribute the UCERF3 forecast earthquake probabilities to the pre-determined (selected) set of scenario earthquakes. Combining the probability of occurrence of the scenario earthquakes with the predicted ground motions, we develop probabilistic estimates of ground shaking over the next 30 years in the Los Angeles basin from San Andreas fault earthquakes.

While the study is conceptually similar to the Southern California Earthquake Center’s CyberShake initiative, the following are critical differences in the methodology: (i) CyberShake utilizes pseudo-dynamic rupture models whereas our source models are kinematic source models of past earthquakes on geometrically similar faults with similar source mechanisms. (ii) The ground motions simulated as part of the last CyberShake release (V15.4) were limited in frequency to 1 Hz (Web Ref: CyberShake (https://scec.usc.edu/scecpedia/Comparison_of_CyberShake_Studies)) whereas the ground motions in this study contain frequencies up to 5 Hz. (iii) While CyberShake source locations and extents follow the UCERF rupture definitions in order to directly use the forecast probabilities therein requiring a far greater number of events (and intensive computational resources) to be simulated, this study describes a method to redistribute the UCERF earthquake probabilities over a smaller set of earthquake sources, significantly lowering the computational demands and making it possible to test and compare alternate approaches. We should also note that the scope of CyberShake is much broader with all known faults being considered, whereas our study is focused on the San Andreas fault.

2. Earthquake Source Models and Associated Probabilities

2.1. Source Models

The San Andreas fault is an almost vertically dipping (dip angle of 90${}^{\circ }$ approximately) right-lateral strike slip fault (rake angle of 180${}^{\circ }$ approximately) with an average seismogenic depth of 20 km approximately. The six earthquakes listed in

Table 1. List of past earthquakes with fault geometry and rupture mechanisms closely matching earthquakes on the San Andreas fault whose kinematic finite source inversions are used in this study. The salient source parameters are listed as well.

	Name	Date	Location	${\mathit{M}}_{\mathit{w}}$	Length (km)	Depth (km)	Dip (${}^{\circ }$)	Rake (${}^{\circ }$)	Reference
1					Denali	2002	AK, USA	7.89		290.0	20.0	90.0	180.0	[13]
2	Izmit	1999	Turkey	7.59	155.0	18.0	90.0	180.0	[14]
3	Landers	1992	CA, USA	7.28	78.0	15.0	89.0	180.0	[15]
4	Kobe	1995	Japan	6.92	60.0	20.0	85.0	180.0	[16]
5	Imperial Valley	1979	CA, USA	6.58	42.0	10.4	90.0	180.0	[17]
6	Parkfield	2004	CA, USA	6.00	40.0	14.5	83.0	180.9	[18]

, with magnitudes $M_w$ 6.0 (2004 Parkfield), 6.58 (1979 Imperial Valley), 6.92 (1995 Kobe), 7.28 (1992 Landers), 7.59 (1999 Izmit), and 7.89 (2002 Denali), occurred on faults that are geometrically similar to the San Andreas fault. It should be noted that except for the Kobe earthquake which is an inland crustal earthquake, all others occurred on plate boundaries. Kinematic finite source inversions of these events, developed by various research groups, are archived in the finite source rupture model database [12]. As evident from Table 1, their source mechanisms are similar to that of a typical earthquake on the San Andreas fault. Seismic source spectrum is closely related to fault geometry and rupture mechanism. Conforming the scenario earthquake source model geometry and rupture mechanism to those observed on the target fault may help produce source spectra that are more representative of earthquakes on the target fault(s). So, we select the source models of these earthquakes as seeds for the scenario earthquakes on the San Andreas fault.

The selected source models are used for both low- and high-frequency contents of the time-histories. The source modelers have made sure that the discretization in these models (Table 1) is sufficiently small to radiate seismic energy at 2 s and longer periods. However, we still have to ensure that these frequencies are coherently propagated to the target sites through the wave-speed model used for the ground motion simulation. Based on the average wave-speed of 3 km/s in the SCEC CVM-H 11.9.0 southern California model and an average rupture speed of 2.5 km/s in all of the source models, we estimate that resampling the source to a finer resolution of 0.5 km would ensure that a 2 s wave (and longer period waves) is reliably propagated to the target sites. This calculation is based on a simplified 1-D source idealization, in which the linear source rupture is divided into segments. The time interval between the wave arrivals at a station from slip on adjacent source segments dictates the shortest wavelength that can be realized at the station. The shorter the time interval, the shorter the wavelength and shorter the wave period. The fact that our source models are 2-D idealizations will only render this estimate conservative as the added discretization of the source along the depth dimension translates to more frequent wave arrivals (from the newly added subfaults along the depth direction) at the station. This, in turn, translates to shorter wavelengths and shorter wave periods being realized at the station. So, we resample each of the source models to a 0.5 km resolution, allocating to the resampled sub-faults the same slips as the original sub-faults (from the original model), and applying a Gaussian filter to marginally smoothen the slip distribution, with 98% of the slip in the original sub-fault being preserved within the resampled sub-faults.

Each of the resampled source models is mapped onto the San Andreas fault at 5 separate locations distributed along the southern section of the fault starting at Parkfield in central California and ending at Bombay Beach in southern California. Two alternate rupture propagation directions are considered at each location, north-to-south and south to north rupture. Figure 1a–j illustrate these mappings for the $M_w$ 7.9 scenario earthquake. A total of 60 scenario earthquakes (six magnitudes × five rupture locations × two rupture propagation directions) are simulated here to cover the broad range of potential San Andreas fault earthquakes that could be damaging to tall buildings in the greater Los Angeles region.

2.2. Scenario Earthquake Probabilities

The Uniform California Earthquake Rupture Forecast (UCERF, [8,9]) by the Working Group on California Earthquake Probabilities (WGCEP), a joint effort of the U.S. Geological Survey (USGS), the California Geological Survey (CGS), and the Southern California Earthquake Center (SCEC), postulates a large set of plausible earthquakes on Californian faults and estimates their annual rates of occurrence on the basis of geologic, geodetic, seismic, and paleoseismic data. All the mapped out faults in California are discretized into segments of 2 km to 13 km length. A plausible event (forecast earthquake) is a hypothetical earthquake that ruptures two or more of these segments. Annual rates of occurrence of all forecast earthquakes are estimated from a grand inversion of diverse datasets of measured fault slip-rates, creep rates, historical earthquake timelines from paleoseismic investigations, etc. The model and data uncertainties are accounted for by the use of a logic tree. The weighted average of the forecast earthquake rates from all branches of this logic tree are converted to time-independent probabilities of occurrence over the target time horizon by assuming a Poisson distribution [Note: in this study we use the latest version of UCERF (Version 3) which provides only the long-term time-independent earthquake rates at the present time].

UCERF-3 postulates 10,445 physically plausible earthquakes (forecast earthquakes) on the southern portion of the San Andreas fault starting from Parkfield in central California to Bombay Beach in southern California (Figure 2). It should be noted that some of these have a near-zero probability of occurrence and their contribution to the hazard is negligible. To estimate scenario earthquake probabilities, all forecast earthquakes with magnitudes between 5.90 and 8.34 whose rupture extent occurs wholly or partially within the southern San Andreas fault are allocated to one of the following magnitude bins: [5.90–6.42], (6.42–6.80], (6.80–7.15] , (7.15–7.45], (7.47–7.78], and (7.78–8.34]. Note that the scenario earthquake magnitudes of 6.0, 6.58, 6.92, 7.28, and 7.59 correspond to centers of the first five magnitude bins based on seismic moment. The last bin is extended to a magnitude of 8.34 to include the probability of the largest plausible forecast earthquake. In the current risk framework, the forecast earthquakes in each magnitude bin will be represented by the ten scenario earthquakes (five rupture locations and two rupture directions) matches the bin’s central magnitude. Accordingly, the UCERF yearly rates of the forecast earthquakes in a given magnitude bin are redistributed among the ten scenario earthquakes representing that bin. This involves converting the forecast earthquake yearly rates to seismic moment rates (multiplying by the seismic moment corresponding to forecast earthquake magnitude), deaggregating the moment rates to the segments being ruptured, summing the moment rate contributions of all forecast earthquakes to each segment, assigning the total moment rate of each segment to the closest scenario earthquake, aggregating the moment rate contributions to each scenario earthquake, and converting the scenario earthquake moment rate to a yearly rate (dividing by the seismic moment corresponding to scenario earthquake magnitude). The scenario earthquake yearly rates are converted to 30-year occurrence probabilities using a Poisson distribution [$P(M_w/loc)=1-e^{-r\Delta T}$, where r is the yearly rate and $\Delta T=30$]. The steps involved are as follows:

(i)

We start with “magnitude binning” of forecast earthquakes. We define as many bins as the number of distinct magnitude scenario earthquakes. The lower and upper magnitude limits of a magnitude bin are the magnitudes derived from the averages of the seismic moments of the corresponding scenario earthquake and the scenario earthquakes tied to the previous and next bins, respectively. Corresponding to the six scenario earthquake magnitudes of 6.00, 6.56, 6.92, 7.28, 7.59, and 7.89, the following six magnitude bins are defined: [5.90–6.42], (6.42–6.80], (6.80–7.15] , (7.15–7.45], (7.47–7.78], and (7.78–8.34]. The seismic moments of the $M_w$ 6.00, 6.56, 6.92, 7.28, 7.59, and 7.89 scenario earthquakes correspond to the average of the seismic moments of the upper and lower magnitude limits of the first five bins, respectively. The upper limit of the last bin is assumed higher to include all forecast earthquakes with magnitude greater than 7.89. Each of the forecast earthquakes will be assigned to one of these magnitude bins. For instance, a forecast earthquake with magnitude, say, between 6.42 and 6.80 will be assigned to the magnitude bin tied to the scenario earthquake with magnitude $M_w$ 6.56. Its probability of occurrence will be redistributed among the ten $M_w$ 6.56 scenario earthquakes (five rupture locations and two rupture directions). The dashed black lines in Figure 3 demarcate the magnitude bins.

(ii)

The seismic moment of forecast earthquake i, $M_o^i$, is multiplied by the UCERF yearly occurrence rate $r_i$ to arrive at what may be termed as the seismic moment release rate with a unit of “seismic moment/year”. Seismic moment release rates are determined for all forecast earthquakes in this manner.

(iii)

Within each magnitude bin, the seismic moment release rate of a forecast earthquake is distributed among the UCERF segments being ruptured by that forecast earthquake in proportion to their areas. This is based on the fact that seismic moment release rate, given by $\mu A\dot{D}$ where $\dot{D}$ is the average slip rate on the fault, $\mu$ is the shear modulus, and A is the area of rupture, scales linearly with segment area (see Equation (4.8) and Appendix G of [19]). Thus the seismic moment release rate contribution of the $i^{th}$ forecast earthquake to the $j^{th}$ UCERF segment equals $r_i{M}_o^i\frac{A_j}{A_i}$, where $A_i$ is the area of forecast earthquake i and $A_j$ is the area of the UCERF segment j.

(iv)

Within each magnitude bin, the contributions to fault segment j of all N forecast earthquakes in that bin are summed: ${\displaystyle \sum _{i=1}^N}{r}_i{M}_o^i\frac{A_j}{A_i}$. This represents the yearly seismic moment buildup in segment j that is expected to be released periodically by earthquakes with magnitudes lying within the bin. It may be termed as the seismic moment release rate for segment j in earthquakes from that magnitude bin.

(v)

Within each magnitude bin, the cumulative seismic moment release rate of segment j, determined in the previous step, is assigned to the scenario earthquake tied to that bin and whose rupture location is closest to segment j. Then the seismic moment release rate of that scenario earthquake is given by ${\displaystyle \sum _{j=1}^M}{\displaystyle \sum _{i=1}^N}{r}_i{M}_o^i\frac{A_j}{A_i}$, where M is the number of UCERF segments occurring within the rupture extent of that scenario earthquake. It is possible that the rupture extents of two or more scenario earthquakes may extend over the same fault segment(s). The moment release rates on such segments are evenly distributed among the overlapping scenario earthquakes.

(vi)

The seismic moment release rate, obtained from the last step, for scenario earthquake k is divided by its seismic moment $M_o^{S{E}_k}$ to obtain its yearly occurrence rate $q_k={\displaystyle \sum _{j=1}^M\sum _{i=1}^N}{r}_i{M}_o^i\frac{A_j}{A_i}/M_o^{S{E}_k}$.

(vii)

The probability of occurrence of scenario earthquake k over a period of $\Delta T$ years is then given by the Poisson distribution as: $P\left(S{E}_k\right)=1-e^{-q_k\Delta T}$.

(viii)

Steps (iii)–(vii) are repeated for all magnitude bins and the scenario earthquakes associated with them.

Figure 4 graphically illustrates the results the various steps of the method applied to a hypothetical fault. The horizontal blue line represents the section of the target fault of interest. For the case study this line would represent the southern portion of the San Andreas fault. The dotted black lines show the segmentation (labeled “Seg.”) of the target fault, say from UCERF. The red unfilled rectangles represent forecast earthquakes (labeled “F.E.”). The heights of the rectangles are proportional to the seismic moment release rates of forecast earthquakes (computed in step ii). The yellow shaded regions represent the seismic moment release rate contributions of several forecast earthquakes to a given UCERF segment (computed in step iii). The red shaded region represents their summation (step iv). The magenta colored region represents the sum of all seismic moment release rates for a scenario earthquake determined by summing the moment release rates of all the UCERF segments within the extent of the scenario earthquake rupture (step v). Magenta lines at the very bottom of the figure show scenario earthquake (labeled “S.E.”) rupture extent.

Shown in

Table 2. UCERF3 time-independent 30-year occurrence probabilities for the 30 scenario earthquakes (six magnitudes and five rupture locations). Half of these probabilities are assigned to north-to-south propagating ruptures and the other half to south-to-north propagating ruptures.

${\mathit{M}}_{\mathit{w}}$ [Bin]	Location 1	Location 2	Location 3	Location 4	Location 5	Total Probability
	(Parkfield)				(Bombay Beach)	(All Locations)
6.00 [5.90–6.42]	0.6449	0.0459	0.1910	0.2485	0.0685	0.8081
6.58 (6.42–6.80]	0.0051	0.0100	0.0854	0.1280	0.0183	0.2288
6.92 (6.80–7.15]	0.0180	0.0171	0.0060	0.0764	0.0271	0.1380
7.28 (7.15–7.45]	0.0211	0.0182	0.0059	0.0153	0.0365	0.0935
7.59 (7.47–7.78]	0.0124	0.0121	0.0061	0.0082	0.0192	0.0568
7.89 (7.78–8.34]	0.0339	0.0281	0.0236	0.0225	0.0215	0.1231
Total Probability
$M_w$ [5.90–8.34]	0.6760	0.1249	0.2904	0.4221	0.1773	0.8553

are the 30-year occurrence probabilities of the 30 scenario earthquakes (six earthquake magnitudes and five rupture locations) determined using this approach. Half of these probabilities are assigned to north-to-south propagating ruptures and the other half to south-to-north propagating ruptures. It should be noted that these probabilities are time-independent and should be treated as probabilities of occurrence in any given 30-year period, not just the next 30 years. Furthermore, based on stress renewal models and elastic rebound theory, the occurrence of any of these events will significantly alter the probabilities in this table. This table is valid only up until the next big event occurs. It should not be used for estimating the joint probability distribution of the random variables, magnitude and location, that would cover all time and locations, and would sum up to unity. It applies to a 30-year period (does not cover all time), does not consider all values of the two random variables, and the total probability does not sum to unity. To ensure that the arithmetic of the outlined method is correct, we compare these values against probabilities derived from event return rates found in literature. The return period for a large magnitude 1857 Fort Tejon-like earthquake ($M_w$ 7.89) is around 150–300 years [20]. Taking the yearly rate (1/225) implied by an average return period of 225 years, its 30-year probability of occurrence may be computed as $1-e^{-(1/225)\times 30}=0.125$, which is close to the total probability value of 0.123 for the (7.78–8.34] magnitude bin in our case study. The return period for the magnitude 6 Parkfield earthquake is about 22 years [21]. The 30-year probability of occurrence (independent of the 2004 event) may be evaluated to be $1-e^{-(1/22)\times 30}=0.74$, which is reasonably close to the total probability value of 0.81 for the [5.90–6.42] magnitude bin in our case study.

3. Ground Motion Simulation

In addition to a mathematical description of the earthquake source, a detailed mapping of the earth’s density and elasticity structure is needed to characterize the seismic wave speeds in the region, allowing for the deterministic simulation of site-specific ground motions. The spatial resolution of this mapping dictates the limiting wavelength (and frequency) of the seismic waves that can be reliably propagated through a finite-element/finite-difference/spectral-element model of the earth; the higher the resolution, the shorter the limiting wavelength and the higher the limiting frequency. Two regional wave-speed models of southern California exist, both developed and maintained by the Southern California Earthquake Center (SCEC): (i) the SCEC Community Velocity Model (CVM, [22,23,24]), and (ii) the SCEC-CVM-Harvard or SCEC-CVMH [25,26,27,28,29,30]. Both models are capable of propagating seismic waves with frequencies at least up to 0.5 Hz and have been used in long-period ground motion simulations in the Los Angeles and surrounding basins (for e.g., [1,31,32,33,34,35,36,37,38], etc.). To synthesize the higher frequencies (above 0.5 Hz) in the ground motion, stochastic (e.g., [39,40]) and empirical (e.g., [41]) methods have been developed. Broadband ground motion is produced by combining these with the deterministic low-frequency ground motion from finite-element, finite-difference, or spectral-element simulations.

Here, we follow the [41] methodology to produce broadband ground motions with frequencies up to 5 Hz. High-frequency seismograms generated using a variant of the classical empirical Green’s function (EGF) approach of summing recorded seismograms from small historical earthquakes (with suitable time shifts) are combined with low-frequency seismograms produced using the open-source seismic wave propagation package SPECFEM3D (V2.0 SESAME, [28,35,42,43]) that implements the spectral-element method. SESAME uses Version 11.9 of the SCEC-CVMH seismic wave-speed model, accounting for 3-D variations of seismic wave speeds, densities, topography, bathymetry, and attenuation. The SCEC-CVMH model incorporates tens of thousands of direct velocity measurements that describe the Los Angeles basin and other structures in southern California [25,30]. It includes background crustal tomography down to a depth of 35 km [44,45] enhanced using 3-D adjoint waveform methods [27], the Moho surface [30], and upper mantle teleseismic and surface wave-speed models extending down to a depth of 300 km [26]. The wave-speed model-compatible spectral element mesh of the Southern California region was developed by [46], who adapted the unstructured mesher CUBIT [47] into GeoCUBIT for large-scale geological applications such as this.

The classical empirical Green’s function (EGF) approach involves the use of aftershock earthquake records as the Green’s functions sampling the travel paths from the source to those stations [48,49,50,51,52,53,54,55,56,57,58]. The rupture plane of an event is divided into (uniform or non-uniform) sub-faults. A pre-selected Green’s function (selected on the basis of the closest match to the subfault-to-target site path) is used to represent the seismic wave radiated from a given sub-fault. The EGFs are selected from a pool of thousands of low-magnitude historic events ($M_w$ 2.5–4.5) that have occurred in the vicinity of the San Andreas fault over the past few decades. Selection was based on signal quality and the scanning of the hundreds of thousands of records was automated. The Green’s functions from all sub-faults are time-shifted and summed to yield the ground shaking at a target site. The key challenge in this approach is that it is difficult to replicate Brune’s spectrum [59] in both the high- and low-frequency regimes simultaneously. Scaling based on seismic moments, where the total seismic moment of the EGFs matches that of the simulated event, will correctly reproduce the low-frequency content of the ground motion. On the other hand, scaling based on areas, where the total area of the EGFs matches that of the simulated event, will correctly reproduce the high-frequency content [51,52]. Ref. [41] was recently successful in developing a variant of the EGF summation that allows for the simulation of high-frequency ground motion (0.5–5.0 Hz) without the use of any artificial filters to achieve agreement with Brune’s spectrum. Because energy release from an EGF is typically much smaller than the moment release associated with slip on a single subfault, the EGF has to be summed several times to match the moment release on the subfault. Furthermore, the duration over which EGF energy release occurs must equal the duration of slip on the subfault, i.e., the EGF energy release must occur over the duration of the source-time function. To achieve this, the EGFs have to be time-shifted and summed. Ref. [41] adopts a non-uniform temporal distribution of EGFs over the source time function of each subfault. This temporal distribution of time shifts is sparse at the start and end of the rupture of that subfault and dense in between. He used low-magnitude ($M_w$ 2.5–4.5) earthquakes as EGFs and combined the high-frequency waveforms generated using this approach with low-frequency waveforms from the deterministic spectral element approach (lowpass-filtered using a second order Butterworth filter with corner at 0.5 Hz) to reproduce ground motions at large distances under the $M_w$ 6.0 Parkfield and the $M_w$ 7.1 Hector Mine earthquakes. We use this hybrid approach to simulate ground motions at the 636 greater Los Angeles sites from the 60 scenario earthquakes on the San Andreas fault.

Figure 5 shows the median values of three commonly used ground shaking intensities, peak horizontal displacement, peak horizontal velocity, and 5%-damped spectral acceleration $S_a$ at 1 s and 0.2 s periods, for the ten ruptures corresponding to each magnitude level of the scenario earthquakes (see the blue lines). The vertical bars show the one standard deviation spread of the data on either side of the median correspond to median. Also shown for comparison are the corresponding values determined using the Campbell-Bozorgnia (CB-08) Next Generation Attenuation (NGA) relation [60]. The soil properties for the 636 sites, as characterized by the $V_S^{30}$ values from [61], the basin depths from the SCEC-CVMH model [30], and the Joyner-Boore distance, defined as the shortest distance from a site to the surface projection of the rupture plane, are used as inputs for the NGA computation. The choice of CB-08 over NGA-West2 was driven by the fact that NGA-West2 excludes peak ground displacements (PGD) due to its sensitivity to frequency-filtering parameters and record processing [62]. Large earthquakes generate long-period ground motion that is characterized by large PGDs and it is important to compare our simulations against those predicted by the GMPEs, especially when considering the implications for long-period structures such as tall buildings.

There is good agreement between simulations and CB-08 in the peak velocity and displacement intensity measures for the lower magnitude earthquakes (up to 6.92). For the larger earthquakes, the simulations predict larger peak horizontal velocities (and much larger variances as well), whereas CB-08 predicts higher peak ground displacements (with comparable variances). CB-08 relies on observed near-field permanent displacements to constrain the PGD attenuation relation. The large permanent ground displacements (up to 9 m) observed during the magnitude 7.6 Chi-Chi earthquake of 1999, one of the few large magnitude earthquakes for which seismic, geologic, and geodetic near-source data is available, may have a strong influence on the PGD attenuation relation. On the other hand, CB-08 relies on seismic data alone for the PGV relation. Unfortunately, there is a sparsity of records from large magnitude earthquakes, especially in deep sedimentary basins such as the Los Angeles basin. This may, in part, explain the differences between the predictions by the simulations and the attenuation relations.

The median values of 1 s $S_a$ predicted by CB-08 are higher for the magnitude 6.00 and 6.58 earthquakes, about the same for the magnitude 6.92, 7.28, and 7.59 earthquakes, and significantly lower for the magnitude 7.89 earthquakes, when compared against the those predicted by the simulations. CB-08 predictions for 0.2 s $S_a$ are higher for the lower magnitude 6.00 and 6.58 earthquakes, but lower for the higher magnitude 6.92, 7.28, 7.59, and 7.89 earthquakes. We should note that the low-magnitude events ($M_w$ 2.5–4.5) used as EGFs are generally deficient in the higher frequencies (higher than about 5 Hz) due to attenuation. Thus, there is a natural tendency for the synthetic ground motions from the hybrid approach to be somewhat deficient in these higher frequencies. Additionally, the two-pass Butterworth filter used in filtering out the higher frequency ground motions has a corner at 0.2 s or 5 Hz. Ground motion intensities fall off smoothly with increasing frequency beyond this filter corner frequency and the values for 0.2 s $S_a$ are probably further under-estimated.

Figure 6a,c,e show maps of the median values of the geometric mean of the horizontal ground velocity under the ten ruptures of the magnitude 7.28, 7.56, and 7.89 scenario earthquakes, respectively. The maps cover the 636 analysis sites in southern California at which ground motions are computed. The corresponding maps, generated using the CB-08 attenuation relations with the site-specific soil and basin depth (Figure 7) information for the 636 analysis sites, are shown in Figure 6b,d,f. The strong influence of the basins is clearly seen. Ground motions are significantly amplified in each of the three basins, San Fernando, Los Angeles (LA), and San Gabriel (SG). The San Fernando valley’s proximity to the San Andreas fault (and perhaps seismic wave-speed structure) results in far more intense shaking there as compared to the LA and SG basins. The simulated ground motions are significantly more intense than the intensities predicted by CB-08, with this difference growing with earthquake magnitude.

Spectral accelerations at 1 s and 3 s periods from the scenario earthquake simulations are compared against those generated using the CB-08 NGA relations in Figure 8. Mean $S_a$ and the one standard deviation spread on either side of the mean are shown plotted as a function of source-to-site distance in Figure 8a,b. The fact that the peaks occur not at shorter distances, but at 35–65 km distances is due to the combined effect of basins (the closest distance to which is about 40 km from the fault) and the Joyner-Boore definition of distance that does not take into account the location of slip asperity on the fault or rupture directivity, being based upon fault proximity alone instead. It is interesting to note that the simulated ground motions carry comparable power at 1 s and 3 s periods. If anything, the peaks in the 3 s $S_a$ are higher than those in the 1 s $S_a$ plots. This is not the case with the NGA predictions with the 3 s period spectral accelerations being significantly diminished when compared to the 1 s period spectral accelerations.

Ref. [7] identified fourteen locations in the greater Los Angeles region where a significant number of tall buildings exist. These include Irvine, downtown Los Angeles, Anaheim, Long Beach, Hollywood, El Segundo, Santa Monica, Century City, Universal City, and Park La Brea in the Los Angeles basin, Encino, Canoga Park in the San Fernando basin, and Glendale and Pasadena in the San Gabriel basin (see Figure 7 for locations).

Table 3. Comparison of ground motion intensities from the ten (five locations and two rupture directions) simulated $M_w$ 7.89 scenario earthquakes against CB-08 NGA predictions at fourteen locations in southern California where a significant number of tall buildings exist.

Site Location	Latitude	Longitude	Simulated								CB-08				Soil Type
			PGV (m/s)		PGD (m)		${\mathit{S}}_{\mathit{a}}^{1\mathit{s}}$ (g)		${\mathit{S}}_{\mathit{a}}^{3\mathit{s}}$ (g)		PGV (m/s)	PGD (m)	${\mathit{S}}_{\mathit{a}}^{1\mathit{s}}$ (g)	${\mathit{S}}_{\mathit{a}}^{3\mathit{s}}$ (g)	UBC	UBC
			$\mathit{Md}$	$\mathit{\sigma }$	$\mathit{Md}$	$\mathit{\sigma }$	$\mathit{Md}$	$\mathit{\sigma }$	$\mathit{Md}$	$\mathit{\sigma }$	$\mathit{Md}$	$\mathit{Md}$	$\mathit{Md}$	$\mathit{Md}$	94	97
Irvine	33.67	117.80	0.46	0.35	0.56	0.29	0.21	0.17	0.15	0.18	0.21	1.18	0.17	0.08	$S_3$	$S_d$
Encino	34.16	118.50	0.30	0.46	0.43	0.24	0.12	0.19	0.13	0.27	0.19	0.89	0.14	0.06	$S_2$	$S_c$
Downtown LA	34.05	118.25	0.79	0.75	0.79	0.60	0.26	0.15	0.38	0.43	0.28	1.70	0.23	0.10	$S_3$	$S_d$
Canoga Park	34.20	118.60	0.94	0.53	0.70	0.41	0.38	0.24	0.30	0.37	0.18	0.87	0.14	0.06	$S_2$	$S_c$
Pasadena	34.16	118.13	0.13	0.09	0.20	0.13	0.04	0.10	0.01	0.03	0.15	0.77	0.12	0.05	$S_3$	$S_d$
Anaheim	33.84	117.89	0.73	0.61	0.70	0.48	0.26	0.18	0.42	0.41	0.22	1.16	0.17	0.07	$S_2$	$S_c$
Long Beach	33.77	118.19	0.26	0.21	0.33	0.27	0.14	0.10	0.08	0.09	0.23	1.38	0.19	0.09	$S_3$	$S_d$
Glendale	34.17	118.25	0.26	0.33	0.40	0.30	0.15	0.09	0.08	0.15	0.20	0.93	0.15	0.06	$S_2$	$S_c$
Hollywood	34.10	119.33	0.31	0.41	0.49	0.27	0.16	0.12	0.19	0.29	0.18	0.85	0.14	0.06	$S_2$	$S_c$
El Segundo	33.92	118.41	0.63	0.39	0.60	0.29	0.18	0.13	0.22	0.19	0.20	1.09	0.16	0.07	$S_3$	$S_d$
Santa Monica	34.02	118.48	0.66	0.32	0.68	0.30	0.17	0.10	0.16	0.12	0.19	0.94	0.14	0.06	$S_2$	$S_c$
Century City	34.08	118.42	0.66	0.41	0.68	0.45	0.16	0.10	0.21	0.16	0.20	0.99	0.15	0.06	$S_2$	$S_c$
Universal City	34.14	118.35	0.27	0.13	0.38	0.19	0.06	0.05	0.06	0.06	0.14	0.54	0.10	0.04	$S_2$	$S_c$
Park La Brea	34.06	118.35	0.30	0.46	0.43	0.24	0.12	0.19	0.13	0.27	0.19	0.89	0.14	0.06	$S_2$	$S_c$

shows the median and standard deviation of the PGV, PGD, $S_a^{1s}$, and $S_a^{3s}$ at these fourteen locations from the ten (five locations and two rupture directions) simulated $M_w$ 7.89 scenario earthquakes on the San Andreas fault. Ground motion is particularly strong at downtown LA, Canoga Park, Anaheim, El Segundo, Santa Monica, and Century City. The corresponding tables for the $M_w$ 7.59 and 7.28 earthquakes can be found in [41].

Figure 9 illustrates the effect of source directivity on ground motions. The north-to-south rupture at location 1 (see Figure 1) directs a great amount of energy into the region of forward directivity, which is the San Fernando valley and the Los Angeles beyond. The south-to-north rupture, on the other hand, directs the energy away from the LA basin into the central valley to the north. The focusing effect is enhanced by the added proximity of the target region to the primary slip asperity in the source in the case of the north-to-south rupture scenario, while the opposite is true for the south-to-north rupture scenario. Note that in reversing the rupture direction, the slip distribution is reversed as well, such that an asperity on the south side of the north-to-south rupture is located on the north side of the south-to-north rupture. Peak horizontal velocity in the target region under the north-to-south rupture scenario is two to four times that under the south-to-north rupture scenario. For scenario earthquakes at rupture location 5, it is the south-to-north rupture that produces the stronger ground motions in the target region and the contrast is comparable to that in the location 1 scenario. The consideration of both N-to-S and S-to-N rupture directivities ensure that the results, when considered collectively, are not biased by unilateral rupture considerations. For example, when a source with strong rupture directivity toward the LA basin is flipped to propagate in the other direction, the latter ground motions are much weaker. When the ground motions from the two events are considered collectively, any bias associated with the former rupture is removed.

The simulated ShakeOut scenario earthquake, used in the Great California ShakeOut Exercise and Drill, is a $M_w$ 7.80 rupture, initiating at Bombay Beach and propagating northwest through the San Gorgonio pass, terminating 304 km away at Lake Hughes in the north. Using a source developed by [63] and the SCEC-CVM wave-speed model [22,23,24], ref. [64] simulated 3-component long-period ground motion waveforms in the greater Los Angeles region. The south-to-north propagating $M_w$ 7.89 scenario earthquake at location 5 [Figure 1j] closely resembles this earthquake in as far as location, rupture directivity, and magnitude (with scenario earthquake having a slightly higher moment magnitude) are concerned. The ShakeOut scenario has served as a benchmark for ground motion simulation methodologies [65] and in Figure 10a-b we compare the results of the simulations here against this established benchmark. The ground motions simulated in this study are more intense than those predicted for the ShakeOut scenario. But the overall pattern of basin amplification is quite similar. The differences may be attributed to the slightly lower magnitude of the ShakeOut earthquake (with 30% smaller energy release) as well as the differences in the source (e.g., peak slip of 16 m in the ShakeOut source versus 12 m in the Denali earthquake source used for the earthquake simulated here) and wave-speed (SCEC-CVM versus SCEC-CVMH) models. The predictions by the NGA relations are far lower [Figure 10c]. The large red blob in the ShakeOut motions, attributed to a wave-guide through Whittier-Narrows by [66], cannot be found in the NGA predictions. Rupture directivity and wave-guide focusing, that clearly may have a strong influence on ground motions, are not explicitly accounted for in the NGA relations. In our simulation, a larger feature encompassing the wave-guide-related feature of the ShakeOut earthquake can be seen.

4. Probabilistic Estimates of Ground Shaking

Having computed ground motions from the suite of 60 scenario earthquakes and having mapped the 30-year probability of occurrence of all plausible San Andreas fault earthquakes on to this set of earthquakes, we can now develop cumulative distribution functions (CDFs) of PGV and PGD. PGV and PGD predictions from the 60 earthquakes at the fourteen tall building locations in southern California are sorted in increasing order and their probabilities accumulated to create the CDFs. Shown in Figure 11 are the CDFs for the fourteen locations. Also marked on the figures are the 10% and 2% probabilities of exceedance levels for PGV and PGD. The PGV and PGD intercepts at all fourteen locations are listed in

Table 4. Probabilistic estimates of PGV and peak ground displacements (PGD) at 14 stations in southern California with probabilities of exceedance of 10% and 2% in 30 years from earthquakes on the San Andreas fault.

Location	PGV (m/s)		PGD (m)
	10%	2%	10%	2%
Irvine	0.16	0.70	0.16	0.89
Encino	0.09	0.46	0.11	0.71
Downtown LA	0.22	0.82	0.19	1.02
Canoga Park	0.16	1.63	0.19	1.22
Pasadena	0.04	0.20	0.09	0.44
Anaheim	0.40	1.38	0.34	1.23
Long Beach	0.08	0.40	0.12	0.63
Glendale	0.05	0.52	0.09	0.74
Hollywood	0.06	0.84	0.10	0.91
El Segundo	0.20	0.89	0.21	0.85
Santa Monica	0.16	0.83	0.20	0.95
Century City	0.18	0.94	0.23	1.05
Universal City	0.06	0.48	0.11	0.74
Park La Brea	0.09	0.46	0.11	0.71

. The following observations can be made: (i) The PGVs at the 10% probability of exceedance level range from 0.05–0.40 m/s, while at the 2% probability of exceedance level they range from 0.20–1.63 m/s. The corresponding PGD ranges are 0.09–0.34 m and 0.44–1.23 m; (ii) For most stations, the variation in the PGV and PGD values at the 2% level is far greater than those at the 10% level. This is reflective of the contrasting long period ground motions that may be produced by larger events influenced strongly by basin structure, rupture directivity, and details of the seismic source (slip distribution, slip velocity and rupture velocity); (iii) The high 2% level PGV of 1.63 m/s and the correspondingly high PGD of 1.22 m at Canoga Park is a result of its relatively proximal location to the fault as compared to the other stations as well as the basin effect associated with the sedimentary deposit in the San Fernando valley; (iv) The high values of 1.38 m/s and 1.23 m at Anaheim, on the other hand, are associated more with rupture directivity and basin effects with Anaheim being located quite a distance away from the fault; (v) Downtown Los Angeles, being located close to the edge of the basin, escapes the harshest effects of basin amplification with 2% level PGV and PGD estimates of 0.82 m/s and 1.02 m, respectively. However, a few cycles at this level of shaking can be quite damaging to the tall buildings located there [67].

5. Discussion and Limitations

Using a single scenario earthquake to represent all forecast earthquakes within a magnitude bin is error-prone by construction. For instance, in the southern San Andreas case study, a magnitude 7.28 scenario earthquake is used to represent seismic risk from all earthquakes with magnitudes 7.15–7.45. To eliminate bias in the results, we have selected the scenario earthquake magnitude to be at the bin center based on seismic moment. However, if an alternate method is adopted to choose the magnitude of the scenario earthquake to be simulated such that the scenario earthquake happens to have a magnitude closer to the upper (or lower) limit of the bin, the lower (or higher) occurrence probability for the scenario earthquake assigned by our seismic moment release rate-based method automatically compensates for the introduced bias. For example, let there be just two forecast earthquakes, one with magnitude 7.15 and a yearly rate of 0.0010, and the other with magnitude 7.45 and a yearly rate of 0.0005. If we use an $M_w$ 7.28 scenario earthquake to represent this bin in our rupture-to-rafters simulations-based case study, our method would result in a scenario earthquake probability of occurrence of 0.045 over the next 30 years. If, on the other hand, an $M_w$ 7.40 scenario earthquake is used to represent this bin, the probability of occurrence drops to 0.030. Obviously, our rupture-to-rafters simulations would predict higher ground motions, heavier building damage, and losses when the 7.40 magnitude scenario earthquake is used. Fortunately, the lower probability of occurrence estimated for the $M_w$ 7.40 earthquake would at least partially offset these increases, perhaps resulting in comparable 30-year losses. Likewise, if we use a $M_w$ 7.18 scenario earthquake to represent this bin, the 30-year probability of occurrence increases to 0.063. This time the lower ground motions and economic losses predicted by the simulations would be at least partially offset by the higher occurrence probability. It is important to note that the errors associated with differences in ground motions from forecast and scenario earthquakes decrease monotonically with increasing density of magnitude and location sampling of scenario earthquakes. This is because the scenario earthquakes sources would monotonically tend to forecast earthquake sources.

The deaggregation of forecast earthquake rates into seismic moment release rates of the segments comprising the rupture solves another commonly encountered problem in PSHA. No matter what domain is chosen for the earthquakes to be considered in the PSHA, at least one or more of the forecast earthquakes will straddle the domain boundary, i.e., only portions of these ruptures will lie within the domain. Where in traditional GMPE-based PSHA, this problem would be circumvented by choosing a region of interest large enough to encompass all of the seismic sources that may contribute significantly to the hazard, it is not practical with rupture-to-rafters simulations-based PSHA. Maintaining the same sampling in rupture location would require many more scenario earthquake simulations and increase the computing resource demand significantly. The question then arises as to what fraction of the probability of occurrence of these earthquakes should be assigned to our closest-occurring scenario earthquake? The deaggregation in our method breaks these ruptures down to the participating segments, and within a small margin of error these segments will lie wholly either inside or outside the domain, thus automatically resolving this problem.

The ground motion estimates in this study have been generated by a limited set of six strike-slip events. The results are likely to be highly influenced by source model selection. To eliminate bias in the results and make them more robust, at each magnitude level, multiple source models from different faults may be used to generate additional scenarios. The size and location of slip asperities, rupture velocity distribution and source-time functions play a critical role in the resulting ground motion in addition to rupture directivity and basin effects considered here. Considering plausible variations in these parameters would lend greater credibility to the ground motion estimates. Likewise, the distribution just 5 rupture locations for the smaller magnitude earthquakes is not dense enough to adequately capture the distribution of shaking from the smaller events. While the sources for the large magnitude events overlap, this is not the case with the small magnitude scenarios. So sampling the full extent of the fault more densely should be a goal for future efforts.

The uniform sampling of scenario earthquakes used in the case study here was chosen for convenience. Sometimes, a non-uniform distribution of earthquakes may be preferred. For instance, one may wish to sample more densely fault regions of high earthquake probability density and/or fault regions whose earthquakes may result in greater variability in basin effects, etc. The outlined method can be used in these cases without any modification.

One of the limitations of our approach is the somewhat incompatible mixing of the highly precise wave propagation simulations using a highly discretized source with the gross nature of the forecast earthquake probabilities in the UCERF model. Three points are to be noted in this context: First, earthquake probabilities, however crude they may be, are needed to make real-world decisions. Only greater data-gathering over time and through ubiquitous instrumentation can help reduce the epistemic uncertainties. Second, the resolution of the data that is used to generate the UCERF model is coarse and good enough to estimate probabilities of occurrence of forecast earthquakes in a smeared manner alone, not high enough to provide recurrence information on fine patches/segments of the source. While this is an inherent limitation of UCERF, it is still the only comprehensive source of earthquake probability information in California. Third, the outlined method does not result in any loss of the existing resolution in earthquake probabilities because segment sizes used in the mapping of the forecast earthquake probability space to the scenario earthquake probability space [in step (iii) of the method] are similar.

Finally, we should note that sometimes there may exist very low probability events that may extend beyond the domain of the rupture-to-rafters simulations. But because their contributions to structural collapse risk may be minuscule it may not be worthwhile to expand the domain of simulation to cover their full rupture extents. Yet their probability of occurrence may be included in the largest magnitude bin of the scenario earthquakes considered. A case in point is a wall-to-wall rupture on the San Andreas fault that extends from the (center of) SAF-Offshore section to the SAF-Coachella section have a recurrence time of 150,000 years. This rupture extends beyond the domain of our simulation to the north. But its probability of occurrence over, say, a period of 30 years is just 2 × 10${}^{-4}$. Even assuming that the probability of a structure collapsing under ground motions from this event is 1.0, this still means that the probability of collapse of the structure over a 30 year period increases by just 2 × 10${}^{-4}$. While the full rupture extent of this event is not considered in the study, its probability of occurrence is included as part of the highest magnitude bin of scenario earthquakes (7.78–8.34].

6. Conclusions

Quantitative probabilistic estimates of ground shaking may be produced by combining 3-D seismic wave propagation simulations of scenario earthquakes with earthquake probabilities from rupture forecasts developed using geodetic and paleoseismic data. The method involves choosing a set of scenario earthquakes spanning the fault of interest, mapping the probabilities of the typically vast set (thousands) of plausible earthquakes from the rupture forecast onto the set of scenario earthquakes, performing deterministic ground motion simulations of the scenario earthquakes, and developing CDFs of PGV and PGD over the targeted period of interest. In this paper, we have presented a case study of ground shaking estimates in the Los Angeles basin from earthquakes on the San Andreas fault over the next 30 years, illustrating the steps involved. Based on sixty simulated $M_w$ 6–8 scenario earthquakes on the San Andreas fault, we estimate the PGVs with 2% probability of exceedance over the next 30 years at the location of 14 tall building clusters to be 0.20–1.63 m/s. The corresponding PGDs are estimated at 0.44–1.23 m. Future work may involve greater number of rupture locations, multiple seismic source models at each magnitude level, bilateral ruptures, alternate methods of EGF selection, the inclusion of earthquakes on all regional faults (not just the San Andreas) and the consideration of time-dependent earthquake probabilities.