ML050390454

Appendix a, Rupture Process of the 2003 San Simeon Earthquake

ML050390454
Person / Time
Site:	Diablo Canyon
Issue date:	11/30/2004
From:	Graves R, Ichinose G, Somerville P URS Corp
To:	Office of Nuclear Reactor Regulation, Pacific Gas & Electric Co
References
DCL-05-002
Download: ML050390454 (61)

Text

Enclosure 1 PG&E Letter DCL-05-002 APPENDIX A Rupture Process of the 2003 San Simeon Earthquake Report prepared for Pacific Gas and Electric Company by Gene Ichinose, Paul Somerville, and Robert Graves URS Corporation

Enclosure I PG&E Letter DCL-05-002 RUPTURE PROCESS OF THE 2003 SAN SIMEON EARTHQUAKE Report prepared for Pacific Gas and Electric Company by Gene Ichinose, Paul Somerville and Robert Graves URS Corporation, Pasadena, CA November 30, 2004

1. INTRODUCTION The objectives of the work described in this report are to develop a rupture model of the 2003 San Simeon earthquake and assess rupture directivity effects on ground motions from the earthquake. The first step was to derive reliable focal mechanisms and source depths of the larger aftershocks from regional long period seismograms. These mechanisms, origin times, and focal depths were then used to calibrate the regional seismic velocity structure, and Green's functions for the velocity structure were calculated for use in deriving a finite fault rupture model of the earthquake. Teleseismic data were used to derive a multiple point source model of the earthquake rupture process. The teleseismic data were then combined with regional broadband, strong motion and geodetic data and used to develop a finite fault rupture model of the earthquake. This rupture model was used in a broadband ground motion simulation procedure to generate ground motion maps, which were used to assess rupture directivity effects in the recorded strong ground motions of the earthquake.

2. SOURCE PARAMETERS OF LARGE AFTERSHOCKS 2.1. Introduction To accurately resolve the distribution of slip of the 2004/12/23 19:15UT (Mw 6.5) San Simeon earthquake, calibrated Green's functions are required which can accurately represent the ground motions for a point source. We accomplish this by inverting broadband records from (3.5<M,<5.5) earthquakes (approximately point sources in the far-field) for the 1-D layered crustal structure.

Before this requirement is accomplished, we first estimate the source parameters including moment tensors of the largest aftershocks in order to compute accurate synthetic seismograms. We use the moment tensor inversion method on long-period regional-waves to estimate moment tensor solutions of 19 aftershocks with M < 4 selected from the Pacific Gas and Electric (PG&E), Northern CA Seismic Net (NCSN) and US Geological Survey's National Earthquake Information Center (NEIC) catalogs.

Although these catalogs are probably complete at the magnitude 4 level, there may be many missing earthquakes from the catalogs early in the aftershock sequence. In addition, we could not analyze any of the aftershocks within the coda of the M, 6.5 mainshock that lasted for about 2 hours2.314815e-5 days <br />5.555556e-4 hours <br />3.306878e-6 weeks <br />7.61e-7 months <br /> duration.

2.2. Data and Regional-Wave Moment Tensor Inversion Methodology Broadband seismograms from the U.S. National Seismic Network (US), Berkeley Digital Seismic Network (BK) and TRINET.(CI) stations (Figure 2.1) are instrument corrected to displacement and filtered between 100 and 5 sec depending on the epicenter distance and event magnitude. These networks operate stations configured with Streckeisen STS-1, STS-2 or Guralp CMG-40 very broadband sensors with 0.01 or 0.033 Hz low frequency corners and can recover static offsets depending on sufficient signal to noise ratio. The data loggers typically record at both high and low gain at different sampling rates ranging from 1 to 40 samples per seconds.

We estimated the source parameters of 19 earthquakes with moment magnitudes between 3.8 and 4.7 (Table 2.1; examples shown in Figures 2.2 and 2.3) using the regional-wave moment tensor inversion methodology similar to that employed by Ritsema and Lay (1995) and Pasyanos et al. (1996). The A-1 PG&E Letter DCL-05-002 method linearly inverts the long-period filtered full waveforms recorded at regional distances for a point source seismic moment tensor approximation.

The seismic moment tensor (Mij) is a second rank symmetric tensor that provides a general representation of the seismic source. The components of the moment tensor represent a set of force couple vectors (dipoles) oriented in the three spatial dimensions aligned along a mutual axis to estimate the torque through the axis perpendicular to the plane containing the force and dipole. The moment tensor is commonly decomposed into three components: double couple (DC), compensated linear vector dipole (CLVD), and isotropic (ISO) components which sum to 100% of the estimated moment release.

Through decomposition of the moment tensor eigenvalues, the DC component consists of two vector dipoles of equal magnitude and opposite sign that resolve shear on faults oriented 450 to the principle eigenvalues of the tensor. The CLVD component consists of a major vector dipole with twice the strength and opposite sign of two orthogonal, minor vector dipoles and can describe the dilatation or compression of a dislocation constrained to have no net volume change. The isotropic component has three orthogonal vector dipoles of equal magnitude that resolves volumetric changes.

40 z38 36 g

KPKD 34 NO.

124 122 120 118 116 Longitude (OW)

Figure 2.1. Map of BDSN (BK) and TRINET (CI) stations mostly used in the moment tensor inversion of the 2003 San Simeon aftershocks.

The moment tensor inversion is performed by applying Green's functions computed assuming a one-dimensional layered earth velocity model. Long-period waves are insensitive to the complexities from finite source contributions and path propagation effects given that the dominant period of the waves is larger than the dimensions of the source and earth heterogeneities. The reflectivity and frequency and wave-number integration technique of Zeng and Anderson (1995) is used to compute the Green's functions.

We use the Western US (WUS) velocity model (C. Ammon, unpublished manuscript) constructed from phase velocity measurements of surface waves across the tectonically active regions of the western U.S. [e.g., Ritsema and Lay, 1995).

We use the local network locations which were either from PG&E, NCSN, or NEIC to compute a suite of Green's functions for 1 km depth increments. Origin times are shifted by increments of 1 sec in A-2 C2q PG&E Letter DCL-05-002 comparing synthetics with observed data. We then iteratively solve for the centroid depth and origin time using a grid search for minimizing the L2-norm objective function while maximizing the percent double-couple component. We solve for the deviatoric moment tensor assuming no net volume change but also solve for a full moment tensor including DC, CLVD, and ISO components in cases where the DC component is low. Over a wide range of frequencies and station distances, the inversion becomes heavily weighted toward higher frequencies and larger amplitudes from near-field stations. We do not apply distance weighting, although there is a natural weighting with distance because the farther stations usually include more data points in the inversion. These inversions are also weighted more heavily by the surface wave contribution than body wave contribution due to the difference in amplitudes.

2.3. Moment Tensor Inversion Results At each step in the grid search for centroid depth and origin time, we solve for the deviatoric moment tensor. The waveform misfit is measured by 3 objective functions, the L2-norm of the raw amplitude error, a percent variance reduction, and cross-correlation coefficient averaged over the 3 components.

The error and percent variance reduction are modulated by the percent double couple. Figure 2.2 illustrates this process for 3 selected earthquakes by showing the deviatoric moment tensor solution and misfit at each trial depth and origin time. The moment tensors are then scaled to the L2-norm raw amplitude error and modulated by the percent double couple. The dark shaded moment tensor is the solution with the best waveform fit with the highest double couple component (Table 2.2). Based on this plot, these earthquakes may have origin time errors of - :1 sec and depth resolutions of between 42 or 44 km. Because of cycle skips, there are also other solutions but the best fit typically shows an improved variance reduction and highest PDC.

Origin Time Relative To Origin Time Relative To Origin Time Relative To 2003/12/23 18:17 UTC 2003/12/25 11:50 UTC 2004/03/17 23:53 UTC 14 -

e 14- (0)0 0

9 D

- 14 @.

6 0

12 -

Q e

12- (D 0 (9 - 12 e

© s

o Q

t0 (8c o

cD~

~ 8

)

0

- 8D tt

( -8(

e 66 E)6

)

gg

• 0 O

~~ Ut3 0

0 Q D(

)9(

2

2.

C (D Z )()-

2. 3 1'0 12 0

2 46 8

Origin Time (sec)

Origin Time (sec)

Origin Time (sec)

Figure 2.2. The moment tensors and waveforni fit errors as a function of depth and origin time are shown for the 3 largest aftershocks selected from the 2003 San Simeon Earthquake. These were obtained by grid search for the best double-couple moment tensor with the lowest waveform misfit. Tjhe moment tensor sizes are scaled by a normalized L2-norm amplitude error and modulated by the percent double-couple (PDC). The best solutions are shown in dark shading of the extensional quadrants.

A-3 PG&E Letter DCL-05-002 TABLE 2.1. Hypocenter and Centroid Parameters of M > 4 aftershocks Date Origin Time AT Latitude Longitude Z'

Z2 M(")

MW(2' Agency Event ID (y/m/d)

(hr:mn:sec)

(sec)

(IN)

(°W)

(km)

(krm) 2003/12/22 20:41:43.3

-1.3 35.669 121.1250 5.00 8

4.40 ML 4.21 NEIC 1005 2003/12/22 21:31:36.2

-0.2 35.692 121.0923 10.3 8

PG&E 21:31:36.4

-0.4 35.711 121.0617 3.88 8

4.33 ML 4.26 NCSN 40148866 2003/12/22 23:52:36.2

+0.8 35.529 120.8773 0.24 6

PG&E 23:52:36.5

+0.5 35.535 120.8805 0.45 6

4.10Mw 4.14 NCSN 30226891 2003/12/23 02:06:55.0

+1.0 35.672 121.1235 9.58 7

PG&E 02:06:55.3

+0.7 35.686 121.1057 6.10 7

4.23 ML 4.16 NCSN 40149054 2003/12/23 02:56:49.2

-0.2 35.709 121.0617 5.24 6

3.94 Md 4.00 NCSN 30226898 2003/12/23 03:46:00.6

+0.4 35.680 121.1310 10.6 5

PG&E 03:46:00.8

+0.2 35.687 121.1268 9.49 5

4.30 Mw 4.29 NCSN 40149108 2003/12/23 05:30:19.4

-0.4 35.660 121.1032 6.38 6

4.21 Md 4.28 NCSN 30226902 2003/12/23 09:41:49.0 0.0 35.525 120.9300 8.00 2

4.10 mb 3.90 NEIC 2003/12/23 11:40:36.1

-0.1 35.673 121.1318 12.0 6

PG&E 11:40:36.5

-0.5 35.688 121.1073 7.39 6

4.18 ML 4.02 NCSN 40149295 2003/12/23 18:17:10.8

+0.2 35.638 121.0623 10.6 5

PG&E 18:17:11.1

-0.1 35.646 121.0548 5.80 5

4.60 MW 4.71 NCSN 40149402 2003/12/25 05:20:13.1

-0.1 35.653 121.1005 8.66 4

PG&E 05:20:13.4

-0.4 35.656 121.0912 2.78 4

4.20 ML 3.94 NCSN 40149777 2003/12/25 05:21:59.0 0.0 35.645 121.1067 8.63 4

PG&E 05:21:59.4

-0.4 35.652 121.0900 4.09 4

4.33 ML 4.01 NCSN 40149778 2003/12/25 09:45:54.8

+0.2 35.657 121.0987 1.02 4

PG&E 09:45:55.3

-0.3 35.664 121.0852 0.10 4

4.10 ML 3.93 NCSN 40149823 2003/12/25 11:50:01.3

+0.7 35.556 120.8435 9.43 10 PG&E 11:50:01.9

+0.1 35.542 120.8395 6.06 10 4.60ML 4.50 NCSN 40149841 2004/01/02 10:44:52.0 0.0 35.697 121.1540 7.98 9

4.00 Md 4.06 NCSN 30226924 2004/01/02 10:47:40.1

-1.1 35.702 121.1507 7.05 8

3.83 Md 4.22 NCSN 30226905 2004/03/17 23:53:06.9

-0.9 35.736 121.0772 8.38 4

PG&E 23:53:07.3

-1.3 35.730 121.0707 3.99 4

4.50Mw 4.57 NCSN 51139882 2004/07/13 23:45:45.3

-0.3 35.735 121.0648 8.78 5

PG&E 23:45:45.1

-0.1 35.725 121.0700 0.04 5

4.05 ML 4.05 NCSN 40159054 2004/07/27 03:38:16.7

-0.7 35.721 121.0565 8.81 4

PG&E 03:38:17.0

-1.0 35.717 121.0602 4.77 4

3.88 ML 3.78 NCSN 40159549 AT=Origin time shift from this study (see Table 2.2), Z(1) Hypocenter depth from Agency, Z(2' Moment tensor centroid depth from this study (see Table 2.2), M(1) Magnitude and type estimated from Agency, MW(2) Moment magnitude from this study (see Table 2.2)

A-4 PG&E Letter DCL-05-002 TABLE 2.2. Moment Tensor Solutions of Aftershocks Date Time(X)

Z M0 x 1022 Mw PDC VRED Nodal Plane 1 Nodal Plane 2 (y/m/d)

(h:m:s)

(km)

(dyne-cm)

(%)

(%)

(strike/dip/rake)

(strike/dip/rake) 2003/12/22 20:41:42 8

2.56 4.21 92 82.0 1520/480/1450 2670/650/480 2003/12/22 21:31:36 8

3.09 4.26 87 84.6 2890/440/720 1340/480/1070 2003/12/22 23:52:37 6

2.02 4.14 99 91.6 3040/310/1120 990/620/770 2003/12/23 02:06:56 7

2.15 4.16 83 93.6 1170/340/77 3130/570/990 2003/12/23 02:56:49 6

1.23 4.00 64 93.4 3390/430/1280 1110/570/600 2003/12/23 03:46:01 5

3.39 4.29 82 96.5 620/31 0/20 3300/890/1210 2003/12/23 05:30:19 6

3.31 4.28 94 93.0 2740/290/740 1110/620/980 2003/12/23 09:41:49 2

0.878 3.90 88 81.9 3580/350/1490 1140/730/590 2003/12/23 11:40:36 6

1.32 4.02 87 90.3 3060/260/1000 1160/640/850 2003/12/23 18:17:11 5

14.3 4.71 97 96.7 2770/420/820 1080/480/970 2003/12/25 05:20:13 4

1.02 3.94 74 88.5 296 /450/910 1140/450/890 2003/12/25 05:21:59 4

1.27 4.01 78 90.0 1250/410/990 2930/490/820 2003/12/25 09:45:55 4

0.979 3.93 86 87.2 3080/410/1000 1150/490/820 2003/12/25 11:50:02 10 7.01 4.50 96 90.0 2920/390/1000 980/520/820 2004/01/02 10:44:52 9

1.53 4.06 92 90.4 1080/410/750 3080/500/1030 2004/01/02 10:47:39 8

2.67 4.22 98 93.1 1190/370/910 2980/530/890 2004/03/17 23:53:06 4

8.85 4.57 92 93.6 3050/360/900 1250/540/900 2004/07/13 23:45:45 5

1.46 4.05 99 82.6 1260/390/1060 2860/530/770 2004/07/27 03:38:16 4

0.592 3.78 98 78.5 1280/460/1050 2870/460/76° Time(l)-Origin time in UTC corrected with centroid time shift, Z-Centroid Depth from this study (Table 2.1)

PDC-Percent Double Couple, VRED-Percent Variance Reduction A-5 PG&E Letter DCL-05-002 2003/12/23 18:17 UTC Mw 4.7 pkd 62 km 660 sao 129 km 3440 pacp 153 km 3520 phl 53 km 1200 sbc 181 km 1370 Z 22.6 J

Z 8.1 Z 6.1 Z 9.6 Z 2.2 R 9.6 R 6.8 R 6.9

=f

=

R 11.5 R 2.5 i-T48.8 T T12.9t T 7.9 T 16.4 T 4.5 smm102km1110 kcc243km400 0

50 0 50 100 seconds Z 3.7 Z 2.9 R 5.8 R 2.8 T 9.0 V\\

T 4.8 2003/12/25 11:50 UTC Mw 4.5 pkd 52 km 450 sao 146 km 3380 pacp 168 km 3460 phl 30 km 1200 smm 80 km 1080 Z 13.5 Z 2.0 =,F,=- Z 1.7 Z 16.3 Z 9.8 R 22.6 R 0.

9 R 1.5 R31.4 R R10.3=ZZ~

T 21.6 Az T 3.6

f T 3.1 T 12.6 T 3.4 A kcc 241 km 340 sbc 160 km 1400 isa 214 km 860 sncc 282 km 1540 Z 0.9 Z 0.6 Z 0.5 Z 0.6 c R 1.3 R 0.9

=

R 0 6

R 0.8 T 1.7T 2.4 T 1.6 T 2.0 0

50 100 seconds 2004/03/17 23:53 UTC Mw 4.6 sao 120 km 3440 pkd 60 km 750 pacp 143 km 3520 sncc 310 km 1520 Z 2.8 Z 17.5246 Z 2.7 Z0.5 R 1.0 R 7.2 R 1.7 -

R 0.6 T 3.7 T 30.2 T 3.8 T 2.2==JJfc kcc 237 km 41 0 smm 107 km 1150 isa 234 km 910 sbc 189 km 1390 Z2.2 Z3.0 Z 0.7 Z 0.6.

R 2.1 R 4.5 74 R 1.2 R 1.0 T 0.8 T 7.1 T 3.1 T

4

=A T 1.5 0

50 100 seconds Figure 2.3. Observed and predicted displacements for the 3 largest aftershocks of the 2003 San Simeon earthquake. The vertical-(Z), radial-(R), and transverse-(T) displacements are in microns. The station code, epicenter distances (kIn), and source to receiver azimuth are labeled. The moment tensor solutions are listed in Table 2.2 and shown in Figure 2.3.

A-6 PG&E Letter DCL-05-002 Mainshock Subevents 35.8 35.7 z

0 0

I-

-I 35.6 35.5 35.4 +--

121.3 I

I I

I 121.2 121.1 121 120.9 120.8 120.7 LONGITUDE (OW)

Figure 2.4. Map of 2003 San Simeon M > 3.8 aftershock moment tensors solutions. Their associated solutions and locations are listed in Tables 2.1 and 2.2. The 2003 San Simeon mainshock sub-event mechanisms #1 and #2 are determined using teleseismic body waves and depth phases. The surface projections of the northeast dipping planes are projected to the surface given the estimated dip and depth.

The depth cross sections A-A' and B-B' are shown in the Figure 2.5.

A-7 Cc?

PG&E Letter DCL-05-002 E 1 000-LIJ-0-

0-10_

A A'

B I'

L elevation exagerated_

m

-1000 CD 1_

- 0

.0 CD

-10g NE SW SW-NE 1i :-I

.I 0

6 io 5a 1(0 Distance (kin) 15 0

1*0 10 20 Distance (km) 30 Figure 2.5. Cross sectional views A-A' and B-B' from Figure 4.2. The 2003 San Simeon mainshock subevent mechanism #1 is shown in the cross sectional view A-A' and subevent mechanism #2 is shown in B-B'. These subevent mechanisms and depths are determined from teleseismic body waves and depth phases but we assume an epicentral location from the NCSN catalog for only subevent #1.

2.4. Analysis of Results The vertical depths of the 19 aftershocks determined from the inversion of long-period regional-waves are on average about 1+/-3 km deeper than those estimated using only phase arrival information from local or regional seismic networks. The dips of the focal mechanisms from the aftershock moment tensors do not steepen with depth. Since the epicenters assumed from the agency earthquake catalogs may have epicentral errors of as much as +2 km, the dipping planes of the seismicity and focal mechanisms may align along the mainshock rupture and back-thrust planes. An alternative explanation is that the faulting is complex and occurring on sub-parallel and conjugate planes (Figures 2.4 and 2.5).

A-8

Enclosure I PG&E Letter DCL-05-002

3. CRUSTAL VELOCITY MODEL CALIBRATION USING BROADBAND SEISMIC WAVES 3.1. Introduction We used aftershocks of the 2003 San Simeon earthquake to calibrate one-dimensional (ID) layered Earth velocity models for the computation of synthetic Green's functions required for the finite-fault slip inversion. We used the waveform inversion method, similar to that described by Randall (1994) and used by Xu and Wiens (1997), to calibrate the ID velocity structure of the crust and upper mantle along source-receiver paths. The premise is that, when the source parameters of earthquakes are well known, including location, origin time, depth, and mechanism, this information can be used to refine the velocity and density structure of the Earth using complete (no windowing of phases) regional broadband records.

We have been successful with this technique for calibrating crustal structure at strong motion sites used for kinematic slip inversions (e.g., Ichinose et al., 2003). The inversion method is also similar to that used to invert receiver functions (e.g., Ammon et al., 1990) except that we use the complete observed regional waveforms instead of deconvolved teleseismic P-waves.

We selected 6 aftershocks of the 2003 San Simeon earthquake to use as calibration events for the estimation of crustal structure of the Central California coast region by inverting complete local to regional (Range 100 krm) seismograms recorded by BDSN and TRINET broadband seismograph networks. These earthquakes, with magnitudes lower than M 5, can be considered as point sources for frequencies below 1Hz and their long period (0.01 to 0.1 Hz) seismograms are first used to estimate their centroid moment tensor, origin time and depth. In the inversion for crustal structure, we first filter the 3-component displacement seismograms for frequencies somewhere between 0.01 and 0.5 Hz. We then invert them for the S-wvave velocities for a laterally homogeneous and isotropic layered Earth model along each source and receiver path. Partial derivatives are computed using a fast reflectivity frequency-wavenumber summation technique, and the conjugate gradient method is used to perform the iterative waveform inversion (Nolet et al, 1986). We perform the inversion for each aftershock at each station independently and then perform a simultaneous inversion of at least 4 events per station. To select a final model we examine the variance reduction between the data and synthetic waveforms computed from models estimated from four independent inversions, simultaneous inversion, and the mean and median models from the independent inversions.

Finally, the resulting model at SAO is used to compute the synthetics at another nearby station located at PACP to show that they can generate realistic synthetic Green's functions up to frequencies of 0.5 Hz. We also show that the model is more appropriate than others used in previous slip model inversions of the 2003 San Simeon mainshock.

3.2. Iterative Waveform Inversion Theory We apply the technique described by Nolet et al. (1986), Shaw (1988), and Randall (1994), which was used by Xu and Wiens (1997) for inverting regional waveforms to solve for the changes in Earth velocities and densities by minimizing the misfit between observed and synthetic differential seismograms. All synthetic seismograms are computed from source parameters including the assumption of the earthquake location, source depth, seismic moment, and focal mechanism. We start with the observed seismograms from a seismic station, do= [ di, d2,..., d. ]T with a total of n time samples. The ID layered Earth model is parameterized as a stack of k=l,...L isotropic layers over a half-space with P-wave velocities ak = [ al, a2,..., aL], S-wave velocities /3k = [ 81

• 42, --,fi L], and densities pk = [ p I p2 -- -,PLI-We define the model vector as, A-9 PG&E Letter DCL-05-002 Mi= [a: P:p I" whereT notes transpose, and there are a total of j--m Earth layer parameters. The differential "residual" seismogram is, r, = di - Si( moJ),

which is the difference between di and synthetic. seismograms Si, computed using an initial model mO.

The partial derivatives make up an mxn matrix, AGV = Si(bmj) - S,(moj) where Em1 are the Earth model perturbations. The ith rows of AGU are constructed using the difference between the synthetic seismogram computed using the perturbed Earth model and the synthetic seismograms computed using the initial Earth model with respect to just thejh initial or perturbed Earth model parameter. From experience, perturbations < 10% of the initial model work best for the fastest convergence. The residual seismograms r can be modeled using a partial derivative matrix AG that maps model perturbations Sm to features in r. This is represented as a system of linear equations of the form r,

= AGU Smj. We solve for Sm using the conjugate gradient method (Paige and Saunders, 1982a; Paige and Saunders, 1982b). The relationship between the Earth model and the seismogram is nonlinear and a solution is often obtained using an iterative sequence of linearizations. We iterate the inversion by updating the initial Earth model using the model perturbations and then solving for new model perturbations. Several iterations are required and an objective function is evaluated until convergence is identified (e.g., Ichinose et al., 2003).

The differential seismograms are computed using a fast f-k summation reflectivity method (Zeng and Anderson, 1995), producing an N/2 fold improvement in speed and accuracy where N is the number of layers compared with a typical finite difference approach. The inversion is performed using a single event and receiver for path specific model calibration, but multiple stations can also be used if the region encompassing the stations is -ID. We make a simplification by solving for only the S-wave velocities since they dominate the regional seismograms. We assume a Poisson's solid to get the P-wave velocities and assume densities from the relation p = 0.32a + 0.77 (Berteussen, 1977).

3.3. Data We gathered broadband seismograms from the Berkeley Digital Seismic Network (University of California Berkeley) and TRINET (Caltech/US Geological Survey). These stations include CIPHL, CISMM, BKSAO, and BK PKD (network designation code/station code). The sites typically operate Streckeisen STS-2, STS-1, or Guralp CMG-40 type velocity seismometers that have a nominal long-period response down to 33 or 100 sec. The stations operate as many as 6 different channels with 2 different gains (low and high gain) and three different sampling rates (1, 20 and 100 samples/sec). We found that the paths to stations adjacent to the San Andreas fault, including SAO, PKD, PHL and SMM, were simpler than paths that cross the Great Valley to KCC and ISA, and that cross the Transverse Ranges to SBC and SNCC, so our study was concentrated on the former stations.

We selected six aftershocks of the 2003 (Mw 6.5) San Simeon earthquake (Table 3.1) that have magnitudes less than 4.7. We notice that the seismograms recorded at the two nearest stations to the source, PKD and PHL, have finite source effects even at the magnitudes of the two largest aftershocks (Mw 4.5 and 4.7). We therefore did not use the largest aftershock in the calibration of PKD and PHL as some of the crustal structure information in the waveforms is filtered out by the effects of finite rupture, A-10

- I ~Enclosure 1

PG&E Letter DCL-05-002 and complex source effects may introduce phase changes which may be mapped into the velocity model.

The source parameters and mechanisms of the earthquakes listed in Table 3.1 were estimated in Section 2 using the moment tensor inversion method from long period regional waves recorded by these same sites but also including other regional broadband stations. Figures 3.1 through 3.4 show maps of the stations and the four aftershocks examined for each site; event numbers are provided in Table 3.1.

TABLE 3.1. Aftershock Source Parameters Date Time"'l Latitude Longitude Z

M"x 1022 Mw Nodal Plane 1 Nodal Plane 2 (y/m/d)

(h:m:s)

(0N)

(OW) km dyne-cm strike/dip/rake strike/dip/rake 1 2003/12/23 03:46:01 35.680 121.1310 5

3.39 4.29 620/310/20 3300/890/1210 2

2003/12/23 05:30:19 35.660 121.1032 6

3.31 4.28 2740/290/740 1110/620/980 3

2003/12/23 18:17:11 35.638 121.0623 5

14.3 4.71 2770/420/820 1080/480/970 4

2003/12/25 11:50:02 35.556 120.8435 10 7.01 4.50 2920/390/1000 980/520/820 5

2004/01/02 10:47:39 35.702 121.1507 8

2.67 4.22 1190/370/910 2980/530/890 6

2004/03/17 23:53:06 35.736 121.0772 4

8.85 4.57 3050/360/900 1250/540/900 Time(l)-Origin time in UTC corrected with centroid time shift, Z-Centroid Depth 38 36 34 0

L0

>1 10 20 ITERATION

-124

-122

-120

-118 Figure 3.1. Location map and mechanisms of four 2004 San Simeon aftershocks (Table 3.1) recorded at station CIPHL (Park Hill). We invert the long period filtered records from these stations to first estimate the moment tensors and then invert these records to estimate the seismic velocity structure. Event E4 lies along a path between the other 3 aftershocks and station PHL. The percent variance reduction %VRED between observed and predicted 3 component waveforms is greatly improved within 10 iterations.

A-1l C3bZ-PG&E Letter DCL-05-002 38 36 34 0w Ofl

>1-0 10 20 ITERATION

-124

-122

-120

-118 Figure 3.2. Location map and mechanisms of four 2004 San Simeon aftershocks (Table 3.1) recorded at station BKPKD (Parkfield). The percent variance reduction %VRED between observed and predicted 3 component waveforms is greatly improved within 8 iterations.

38 36 34 0w Wl

>1N

-124 I I.

I I

I I

I I I

-122

-120 10 X

ITERATION

-118 Figure 3.3. Location map and mechanisms of four 2004 San Simeon aftershocks (Table 3.1) recorded at station BKSAO (San Andreas Observatory).

The percent variance reduction (%VRED) between observed and predicted 3 component waveforms is greatly improved within 8 iterations.

A-12 PG&E Letter DCL-05-002 38-0-

36 hi&4 L

O;

-10 IS'M1-E2 34MN 47 Mi SMVVE3 34 SMvMv-E4

_SME6 200 km

\\ '<

\\. \\ t --

200-1 2 -

.l U

I*~

0 10 20

-124

-122

-120

-118 ITERATION Figure 3.4. Location map and mechanisms of four 2004 San Simeon aftershocks (Table 3.1) recorded at station CISMM (Simmler). The percent variance reduction (%VRED) between observed and predicted 3 component waveforms is greatly improved within 4 iterations for events 2 and 4 while event 3 diverged from a local minimum after 6 iterations.

3.4. Procedure We instrument correct the seismograms to displacement and band-pass filter the waveforms for two pass bands, 0.01 to 0.2 Hz (100 to 5 sec) and 0.01 to 0.5 Hz (100 to 2 sec) using a two pass, 3 pole Butterworth filter. The seismograms are then aligned given the origin time determined from the moment tensor inversion and rotated into vertical, radial and transverse components.

We first inverted the 3-component displacement seismograms for each of the four stations separately to generate four different S-wave velocities along each source and receiver path. From these four models we compute an average and median models. We then invert all of the 3-component waveforms using the same four events simultaneously recorded at the same site. We compute a percent variance reduction

(%VRED) by (d (t) -

s(t)) 2

% VRED

=

2-2100 d (t) 2 where d(t) is the observed time series, s(t) is the synthetic time series. The model which generates synthetic seismograms with the largest variance reduction (largest positive values in Figure 3.5) with respect to the recordings is selected as the general final model for that site. When using data from various events or stations in a simultaneous inversion, the issue of weighting becomes important. We did not apply any weighting or normalization and therefore the data from the largest earthquake usually dominated the solution as the inversion tries to fit the waves with the largest amplitudes. The magnitude differences within each bin of aftershocks only varied by a factor of 2 to 6, so we obtained velocity models that are very similar to those obtained in single station inversions.

We only invert for the S-wave velocities because the amplitudes from shallow earthquakes are generally dominated by S-and surface-waves in this frequency range at regional distances. We also A-13 PG&E Letter DCL-05-002 assumed that the earth is approximately a Poisson solid to fit the P-waves, which usually have lower amplitudes. We assume a Poisson's solid to get the P-wave velocities and assume densities from the relation p = 0.32a + 0.77 (Berteussen, 1977). We also assume a generic intrinsic attenuation profile which usually is not a significant factor for frequencies below 1 Hz relative to the frequency independent geometrical spreading that is taken into account in the computation of the synthetic seismograms.

BK SAOf=0.01-0.21-z BK SAOf=0.01-0.51-100 1 '

a a i I

I I

I I

aPHLf=.01-o.21F CIPILf=0.01-0.51-0I

-lb W

e4 n

)

W.

o-4 en a)

M)

X Cn M)

M co Cn M

NO1

• 'O1 O r,4, N

N NN N NN N

N N N

-At-EVBM Mca

-e-eAS'J4m&

-^ -

m EAN MOM SBMTAH=J 1WNVEL

-A-EvNT1 MODEL EVBCMMDEL E\\BT4 MVODEL

-e -

OD KEMMODEL XEJANADU CIM+/-TANEOUS INWMCD Figure 3.5. Variance reduction obtained using different velocity models estimated using one of four events and a simultaneous inversion using four events within two frequency bands. (Left) The highest variance reduction at station SAO for all of the data (Events 2, 3, 4, and 6) was obtained using the mean model while the reduction for the model obtained from the simultaneous inversion of all of the data was only a few percent lower.

We selected the mean model as the final model at station SAO since it performed the best up to frequencies of about 0.5 Hz. (Right) The highest variance reduction at station PHL was obtained using the velocity model obtained from inverting all of the data simultaneously. The simultaneous and mean models were the only ones that produced reasonable variance reductions for all of the events. We selected the simultaneous inversion model as the final model at PHL.

The velocity model is composed of 14 layers usually with 2-4 km thickness with the exception of the top layers having thickness of 1 km. Ideally all of the layers should have about the same thickness to avoid under or over weighting the importance of each layer's contribution to the velocity perturbations, A-14 PG&E Letter DCL-05-002 but the changes in thickness are to account for the realistic decrease in heterogeneity of the crust with increasing depth with the exception of the area near the Moho discontinuity. The increasing thickness as a function of depth may also be necessary in order to not over emphasize the lack of resolution in parts of the model which do not have penetration of body waves but are still constrained by surface waves. In previous tests of Japanese and central US earthquakes, we found that models of only 5 layers resulted in velocity profiles that had an averaging effect when compared velocity profiles resulting from inversions using 14 to 21 layers.

Synthetic seismograms are first computed using an initial model. In most cases we used a modified CC model where we changed the crustal thickness where needed. The PG&E relocated hypocenters were used to constrain the epicenter location (usually less than 2 km from the ANSS catalog location). The other source parameters including mechanisms were taken from Table 3.1. Location errors can introduce a 1-6% change in overall velocity perturbation.

0 l I

.0 zmax(P)

Zmax(P)

-10 hi

-10

-20

-20 a-30 BKPKD

-30 C-_PHL O

__CC mode[l Ic M ol

-40 Mean

- Mean 5Simul Inv

._Simul v1

-50

.- 50 1

2 345 6

7 89 1

23 4

567 8

9 0

0I

-10 Zmax(S)

-10 Zmax(S)

E

~--20 zmaP)

-20 Zma(P)

-30 BK_SA

-30 CISMM O

-CC Mod

-CC model Mean Mean Simul Inv Simul Inv

-50

-50:

1 2

3 4

5 6

7 8

9 1

2 3

4 5 5

6 7

89 P-and S-Wave Velocity (km/s)

P-and S-Wave Velocity (km/s)

Figure 3.6. P-and S-wave velocity profiles from models which give the highest and second highest variance reduction shown in Figure 3.5. P-and S-wave velocity profiles using six aftershocks of the 2003 San Simeon Earthquake and four broadband seismic stations. We performed an iterative waveform inversion using four aftershocks individually and compute the mean velocity models all four inversions at each site (labeled "Mean"). We also invert all four aftershock waveforms at each site simultaneously (labeled "Siml Inv").

Zmax(P) and Zmax(S) denote the bottoming depths for P-and S-waves respectively. For comparison we also show the Central Coast (CC) model (McLaren and Savage, 2001).

Generally the wave velocities are higher in the upper crust (<15 km) except for the top 1 km. There are also crustal thickness variations between SAO to the north and SMM to the south. The models at PHL and PKD are not resolved below Zmax because of the small epicentral distance.

3.5. Results In Figures 3.7 through 3.10, we show the 3-component displacements compared to synthetic seismograms computed using the mean or simultaneous inversion models in Figure 3.5 and selected in Figure 3.6. The observed and synthetic displacements were bandpass filtered between 0.01 or 0.02 sec and 0.5 Hz using a 3 pole, 2 pass, Butterworth filter. The final models are listed in Tables 3.2-3.5.

A-15 C (o PG&E Letter DCL-05-002 EVENT #1 vertical Mw 4.3 ClPHL 57km/1 190 radial transverse 0.005 I

I I I I

]

I b

Iccnn E I sinrul inv edata Event #4 vertical Mw4.5 CIPHL 30km/1200 radial transverse 0.02 -

0.00 I

cc nodel sirrul inv

\\ data Event #2 vertical Mw4.3 CIPHL 63km/1210 radial transverse 0.008 -

o 0.000 -j 0.006 -

0 0.000 -j Event #5 vertical cc rrdel sirrul inv data cc rnodel

-sirrul inv data I

l I lI l

0 20 40 60 seconds Figure 3.7. Comparisons between the 3-component predicted and observed displacement seismograms at station PHL using the final velocity model for this site (Figure 3.6). The synthetics are computed for four different aftershocks recorded at the station (Table 3. 1). The dashed lines mark the approximate P-and S-wave arrival times. The synthetics computed using the CC model (McLaren and Savage, 2001), similar to the initial model, are shown for comparison. Arrows point to amplitude and phase arrivals that were significantly improved using the new velocity model relative to the CC model.

A-16 C37 PG&E Letter DCL-05-002 Event #4 Mw4.5 BKPKD D=30km f =1200 vertical radial transverse 0.02 -

E 0.00 -

Event #2 Mw 4.3 BKPKD D=5 vertical radial 0.009-I I I

I E

I l

I I

0.000- I Event #1 Mw vertical 0.006 -

U I

Event #5 MV vertical A, 0.009 II]

7km f =1 190 tranerse I

cc model

_~~

hsMUl i nv data 51 km f =1210 transverse X

\\

S cc model I

VIl si mul i nv

- 4 data 33km f =1210 transverse I I I I cc model

-sisimul inv

-KV nIAA/

data I

I I

I I

I I

0 20 40 60 seconds Figure 3.8. Comparisons between the 3-component predicted and observed displacement seismograms at station PKD using the final velocity model for this site (Figure 3.6). The synthetics are computed for four different aftershocks recorded at the station (Table 3.1). The dashed lines mark the approximate P-and S-wave arrival times. The synthetics computed using the CC model (McLaren and Savage, 2001), similar to the initial model, are shown for comparison. Arrows point major amplitude and phase misfits that were significantly improved using the new velocity model relative to the CC model.

A-17 PG&E Letter DCL-05-002 Event #4 Mw4.5 BK_SAO D=146 km f =3380 vertical radial transverse 0.006 -

0.000 -

cc model n'ean data Event #2 Mw 4.3 BK_SAO D=127 kmf=3460 vertical radial transverse

-cc model 0.015 -

0.02 -

0.00 -_

-cc rmdel Event #3 Mw4.7 BKSAO D=1 29km f =3440 vertical I

radial transverse

-cc model I

I I

I I

I I

I I

0 20 40 60 80 seconds Figure 3.9. Comparisons between the 3-component predicted and observed displacement seismograms at station PKD using the final velocity model for this site (Figure 3.6). The synthetics are computed for four different aftershocks recorded at the station (Table 3.1). The dashed lines mark the approximate P-and S-wave arrival times. The synthetics computed using the CC model (McLaren and Savage, 2001), similar to the initial model, are shown for comparison. Significant amplitude and phase misfits between direct S-waves, surface-waves, Moho reflected S, and any possible S-to-P converted phases were improved using the new velocity model relative to the CC model.

A-18 PG&E Letter DCL-05-002 Event #4 M, 4.5 C ISMM 80km/1 080 vertical radial transverse

^

fA u.u I -

n nn -

cc mode simul inv data v.WW Event #2 Mw 4.3 CI SMM 107km/11 0 vertical radial transverse 0.01 I

0K0 I'm [vV cc model o

Aftsi mul i nv 0.00]

data II FV dt Event #3 Mw4.7ClISMM 102km/1110 radial tran 0.02 -

0 0.00 -

0.02 -

0.00 -

cc mode sirrul inv data Event#6 Mw4.6ClISMM 107km/1150 verti cal radial transverse Pg U

I AT 1

Ie 0

20 40 60 80 seconds Figure 3.10. Comparisons between the 3-component predicted and observed displacement seismograms at station SMM using the final velocity model for this site (Figure 3.6). The synthetics are computed for four different aftershocks recorded at the station (Table 3.1). The dashed lines mark the approximate P-and S-wave arrival times. The synthetics computed using the CC model (McLaren and Savage, 2001),

similar to the initial model, are shown for comparison.

Significant improvements were only made to fitting data recorded by events #3 and #4.

A-19 c V-0 PG&E Letter DCL-05-002 3.6. Discussion The resulting S-wave velocity structure is best resolved to a "bottoming" depth before the waves turn or reflect and refract. We estimate a turning depth of P-and S-wave from the TauP software package (Crotwell et al., 1999) using the final velocity models. Derivatives of differential seismograms are usually zero for layers below this depth. Higher mode surface waves can have some sensitivity below this depth but the fundamental mode sensitivity eigenfunction also decays with depth but depending on epicenter distance and hypocenter depth, surface waves may not have a chance to develop. The resulting waveform fits are improved significantly relative to the initial model, especially for the P-and S-body waves, and Rayleigh, and Love surface-waves relative to the initial model shown by the variance reduction curves in Figures 3.1-3.4.

The improvement of the resulting seismic velocity models is demonstrated in the case of SAO by showing that the waveform fits obtained at station PACP near SAO are superior using the velocity model for SAO than using other regional crustal velocity models.

All of the P-and S-wave velocity model profiles (Figure 3.6) indicate a linear gradient in the shallow crust between approximately 0 and 6 km depth with a variation between 3 and 6 km starting from a P-wave velocity of about 3.5-4 krn/s to 5.8-6.3 km/s. This velocity range is very similar to that found by Walter and Mooney (1982) for the Gabilan Range, indicating sediments overlying granitic basement. The sharp gradient at 3-6 km depth causes a triplication that generates complex P-and S-wave arrivals for earthquakes that occur at a depth above the sharp gradient. The multiple high frequency P-and S-wave arrivals (due to triplications) are best observed in velocity seismograms but can also be observed in the displacement seismograms shown in Figures 3.7 through 3.10, especially at station SAO.

This triplication can cause wave propagation effects that could be misinterpreted as a source effect in a kinematic type slip inversion using local strong motion data.

The crustal P-and S-wave velocities (Figure 3.6) between 6 and 26 km, below the sharp velocity gradient, is fairly smooth (between 6 and 6.5 km/s), consistent with the Walter and Mooney (1982) average of 6.1 to 6.5 km/s between 2 and 24 km depth. All the stations see this increase in velocity between 8 and 15 km depth relative to the CC model, which approaches 6.1 km/s at 15 km depth. This probably represents granitic rocks (e.g., Quartz Monzonite) found in the Salinian block and has significant implications for the tectonic evolution of the Central Coast Ranges, e.g., Walter and Mooney (1982), Howie et al. (1993) and Hauksson et al. (2004). The small low velocity zones (LVZs) in the upper to mid crust may be artifacts resulting from the fitting of ID models to actual 3D structure.

Alternatively, the LVZs may indicate that Franciscan type rocks are embedded in the host rock around 10 km depth.

Stewart (1968) originally observed amplitude difference ("shadow effects") along seismic refraction profiles which were confirmed by Walter and Mooney (1982) in their reinterpretation of these data. The Moho depth and lower crustal gradients are not well constrained, so we assume a Moho depth between 26 to 30 km. This does not greatly affect the regional long period waveforms at the source to receiver distances of our data, but is important for short period arrival time fitting of phase date in earthquake location programs.

To validate the SAO mean model (Table 3.2 and Figure 3.6), we used this model for a blind prediction of the ground displacements recorded at station PACP (Figure 3.11) about 20 km northeast of the SAO station. We find that this model can predict the observed seismograms very well in body wave amplitude and phase up to epicentral distances of 150 km in all four aftershock records filtered between 0.01 and 0.5 Hz. There is a phase shift in the surface waves probably caused by the crustal transition between the Gabilan to Diablo Ranges along the San Andreas Fault zone.

Figure 3.12 shows synthetic seismograms computed at station SAO using the final velocity model compared to synthetic seismograms computed using various models (Figure 3.13) including the CC and Southern California model (SC) up to frequencies of about 0.3 Hz. These models do not perform as well in fitting the amplitude and phases of body and surface waves, and can cause artifacts in slip model inversions.

In particular, there are major differences along paths to PKD, a station heavily used in previous slip model inversions by Ji et al. (2004).

A-20 PG&E Letter DCL-05-002 Table 3.2. BKSAO Mean velocity model Thicknes ZTOP a

Q.

p p P s (km)

(km)

(kni/s)

(km/s)

(g/cm3) 1 1

0 3.53 50 2.04 25 1.90 2

2 1

4.55 50 2.63 25 2.23 3

4 3

5.86 200 3.38 100 2.65 4

4 7

5.62 200 3.24 100 2.57 5

4 11 6.39 200 3.69 100 2.81 6

4 15 6.37 200 3.68 100 2.81 7

4 19 6.38 200 3.68 100 2.81 8

4 23 6.46 200 3.73 100 2.84 9

4 27 6.76 200 3.90 100 2.93 10 4

31 7.79 500 4.50 250 3.26 11 4

35 7.77 500 4.49 250 3.26 12 4

39 7.76 500 4.48 250 3.25 13 4,

43 7.76 500 4.48 250 3.25 14 x0 47 7.74 500 4.47 250 3.25 Table 3.3. CISMM Simultaneous inversion velocity model Thickness ZTOP a

Q.

1 Qp p

(km)

(km)

(km/s)

(kmns)

(g/cm3) 1 1

0 3.86 50 2.23 25 2.01 2

1 1

4.23 50 2.44 25 2.12 3

2 2

4.71 200 2.72 100 2.28 4

3 4

5.95 200 3.44 100 2.68 5

3 7

5.84 200 3.37 100 2.64 6

4 10 6.37 200 3.68 100 2.81 7

4 14 6.26 200 3.61 100 2.77 8

4 18 6.00 200 3.46 100 2.69 9

4 22 5.99 200 3.46 100 2.69 10 4

26 7.84 500 4.53 250 3.28 11 4

30 7.76 500 4.48 250 3.25 12 4

34 7.76 500 4.48 250 3.25 13 4

38 7.90 500 4.56 250 3.30 14 X

42 7.99 500 4.61 250 3.33 A-21 PG&E Letter DCL-05-002 Table 3.4. BK PKD Simultaneous inversion velocity model Thickness ZTOP a

Q.

D QA p

(km)

(km)

(kmls)

(kimis)

(g/cm) 1 1

0 3.38 50 1.95 25 1.85 2

1 1

4.26 50 2.46 25 2.13 3

1 2

4.39 200 2.54 100 2.18 4

1 3

4.70 200 2.71 100 2.27 5

2 4

5.69 200 3.28 100 2.59 6

2 6

6.33 200 3.66 100 2.80 7

2 8

6.21 200 3.58 100 2.76 8

5 10 6.38 200 3.68 100 2.81 9

5 15 6.21 200 3.58 100 2.76 10 5

20 6.41 200 3.70 100 2.82 11 5

25 6.60 200 3.81 100 2.88 12 10 30 7.29 500 4.21 250 3.10 13 10 40 7.78 500 4.49 250 3.26 14 Xo 50 7.73 500 4.46 250 3.24 Table 3.5. CIPHL Simultaneous inversion velocity model Thickness ZTOP ax (km/s)

Q.

Q p

(km)

(km)

(kimis)

(g/cm3) 1 1

0 4.61 50 2.66 25 2.24 2

1 1

4.37 50 2.52 25 2.17 3

1 2

4.34 200 2.51 100 2.16 4

1 3

4.40 200 2.54 100 2.18 5

2 4

4.95 200 2.86 100 2.36 6

2 6

6.14 200 3.54 100 2.73 7

2 8

6.00 200 3.46 100 2.69 8

5 10 6.48 200 3.74 100 2.84 9

5 15 5.77 200 3.33 100 2.62 10 5

20 6.53 200 3.77 100 2.86 11 5

25 7.18 200 4.15 100 3.07 12 10 30 7.44 500 4.30 250 3.15 13 10 40 7.77 500 4.48 250 3.26 14 xo 50 7.67 500 4.43 250 3.22 A-22 PG&E Letter DCL-05-002 EVENT-4 2003/12/25 11:50 Mw4.5 BKPACP 168 km/346 0 0.006-Z I.

IA A A A EVENT-6 2004/03/17 23:53 Mw 4.6 BKPACP 143 km/3520 0.015 P

S I

I I

I I

I I

I 0

40 80sec EVENT-2 2003/12/23 05:30 Mw4.3 BKPACP 151 km/354 0 0.006-oI, -

\\AA (cm)

I__

~wAA 4A N

EVENT-3 2003/12/23 18:17 Mw 4.7 BKPACP 153 km/352 0 0.02(c (cm)

-§28~h I

I*

data P S I

I I

I I

I I--

n 0

40 80sec Figure 3.11. Forward calculated synthetics at station BKPACP using the model in Table 3.1 for four San Simeon aftershocks. The synthetics and data are band pass filtered between 100 and 2 sec period (0.01 to 0.5 Hz). Predicted P-and S-wave arrival times from the models are marked. The records were not used in the velocity model inversion. This represents a blind prediction at a station near BKSAO.

A-23 PG&E Letter DCL-05-002 2003/12/23 18:17 UT San Simeon Aftershock (Mw 4.6)

SAO 128km/3440 Pg PPSg Rayleigh It l1 synthetic-SC Model synthetic-CST Model synthetic-CC Model data synthetic-SAO Model synthetic-SC Model synthetic-CST Model synthetic-CC Model E

E 0.0 C u

'a 0.00-4 data synthetic-SAO Model synthetic-SC Model

-synthetic-CST Model

-synthetic-CC Model data synthetic-SAO Model I

I I

0 20 40 60 seconds Figure 3.12. Comparisons between waveform displacements recorded at BDSN station SAO (R=128km, AZ=344 0) and predicted synthetics computed using various models. The CC-model is derived from a ID joint hypocentral determination inversion by McLaren and Savage (2001). The CST-model was derived from a trial and error procedure by Poley and Eaton and is used by the Northern California Seismic Network for locating earthquakes in this region (Oppenheimer et al., 1993). The SC-model (Dreger and HeImberger, 1993) is the southern California velocity model.

The CST and SC models are normally used to compute ID Green's functions for slip model inversions and routine location of earthquakes but since they are constrained from only band limited P-waves, they might not be appropriate for modeling S-and surface waves.

The preferred SAO model obtained in this study (Table 3.2) is estimated from the iterative waveform inversion of complete waveforms filtered between 0.01-0.5 Hz. This model performs much better in predicting the amplitude and phase arrival times of many crustal phases.

A-24 PG&E Letter DCL-05-002 CC (M&S2001)

SC (D&H1993)

-I CST (Poley&Eaton)

-~

SAO (Table 3.2) 5-0a 0

2 3

4 5

6 7

8 P and S-velocity (km/sec)

Figure 3.13.

3.12).

1-D layered Earth velocity models used in the waveform comparisons (see previous Figure A-25 C43 PG&E Letter DCL-05-002

4. MAINSHOCK RUPTURE MODEL FROM TELESEISMIC DATA 4.1 Introduction After the 22 December 2003 San Simeon earthquake, several investigators (e.g. Hardebeck et al.,

2004) examined the rupture process using preliminary information that may not be accurate or well resolved. For example, mainshock hypocenter depth is usually poorly determined in routine earthquake locations, and precise relative locations cannot give information about absolute depth of faulting unless the master or reference event is well located. Inversion of waveform data for the slip distribution also assumes the mainshock hypocenter depth and may be biased by this assumption. The main purpose of this report section is to determine the depth and examine the rupture process of the San Simeon earthquake. In this section, we have only used teleseismic data, but a more comprehensive analysis using strong-motion, regional and teleseismic waveforms and geodetic data is presented in Section 5.

4.2 Method We use the method of Kikuchi and Kanamori (1991) that inverts teleseismic body waves to determine the mechanism and simple rupture pattern of large earthquakes. We call this method the multiple-point-source moment tensor inversion (MPSMTI or Kikuchi method). The rupture pattern is represented as a sequence of point sources (subevents) distributed in space but each with a finite duration in time. The subevent mechanisms and sizes are allowed to vary during the sequence. We use observations on rupture directivity and aftershock locations, described below, to limit the number of unknowns by constraining the area of rupture and resolve the fault plane ambiguity of the rupture.

4.3 Data We gathered global teleseismic data from the IRIS Data Center and Canadian National Seismic Network stations (Figure 4.1) which consisted mainly of 3 component broadband digital waveforms between epicentral distances of 300 and 1000. These were instrument corrected to velocity and also integrated to displacement within the band of 100 second period to 10 Hz. The data are decimated from 0.05 or 0.01 sec/sample to 0.2 sec/sample (Nyquist frequency of 2.5 Hz). The decimation does not greatly affect the data bandwidth because attenuation usually limits most of the teleseismic energy to frequencies below 2 Hz. We used travel-time calculations from the IASPEI-91 global velocity model to identify the arrival times of P-and S-waves and then window about 25 sec of P-waves. The P-wave polarities from the teleseismic dataset are displayed in Figure 4.2 along with some far-regional stations to show the change in polarity expected at the nodal plane.

4.4 Results We use observations on rupture directivity and aftershock locations to constrain the area of rupture and resolve the fault plane ambiguity of the rupture. One constraint is the relatively strong directivity effect seen in the teleseismic P-waves. In Figure 4.3, the time between peaks along azimuths to the southeast in the forward rupture directivity direction, at stations such as PAYG.Z (Galapagos Islands), is shorter than for azimuths to the northwest in the backward rupture directivity direction, at stations such as MAJO.Z (Matsushiro, Japan; Figure 4.1). These observations suggest two subevents separated in time by about 1-4 sec (average approx. 3.3 sec) depending on the source-receiver azimuth. We fit the time difference between the two peaks in the displacement spectra using a pulse width directivity equation based on the Doppler effect, 8r= L(1/V, - cos((p-4)/c) where time separation between peaks 8r is a function of L the fault length (or subevent separation), V, the rupture velocity, c the medium velocity, and ((p4) the difference between the strike of the rupture direction or fault and the azimuth of the station from the epicenter. We fit the time difference using a range of fault lengths (6-16 km) and rupture velocities (2-4 km/s) to estimate the best fitting fault strike or A-26 PG&E Letter DCL-05-002 rupture direction assuming an average P-wave crustal and mantle velocity of c = 9 km/s. We fit the data using both minimum deviation Li-norm and least squares (Figure 4.3). The error curves show the best rupture direction of 1350 using the Li -norm (solid line) and 1150 using the L2-norm (dashed line). Based on this analysis we will examine the effect of rupture direction on the fit of the waveforms for azimuths from 1000 to 1600. Note that the low rupture velocity of 2-2.5 km/s may be an artifact of using a simplification of a single medium velocity. Additional analysis using velocities based on ray bottoming depths or ray parameter values for individual stations may improve this estimate.

The three component strong motion records at Cambria Hwy #1 (Figure 5.5; about 13 km away) also indicate that there were two major releases of seismic energy (two subevents) visible in the acceleration, velocity and displacement records. In the displacement and velocity records, the time difference between these two pulses (S-waves) is about 3-3.5 seconds. Since the source directivity effect is seen consistently at all teleseismic stations, we can rule out the possibility that this phase was caused by some unknown near source Earth structure. We also use information from 1-day and 1-week aftershock locations to constrain the area of rupture and resolve the fault plane ambiguity of the rupture, which appears to be 20-30 km in length and 15-18 km in width given a fault dip of 60-50°.

LLN HIA YAK COLA PET4%

KDAK ALE KBS KEV KONO BJT u

42 ocESK ARU BRC OBN KlV par o

\\

~~~SN'&O

__TB 900 Figure 4.1. Teleseismic station locations (triangles) used in the multiple-point source moment tensor inversion. The dashed lines indicate the distance ranges of 30°, 60°, and 90° from the epicenter (star).

COA

,TNFigure 4.2.

Focal mechanism and first K[>4 \\22

/

/9 FCCmotion polarities of far-regional and teleseismic distance stations.

The focal PET ISCO mechanism with the shaded extensional YCM quadrant was estimated from the multiple-

P-CXF point source moment tensor inversion of JOT teleseismic body-waves, and represents the moment release of both subevents.

The IA \\

_dashed lines represent the moment tensor solution of the first subevent.

A-27 Cqht PG&E Letter DCL-05-002

_0.5 -.

.,.l 2

pl p2 L=7 km L=8km E 0.4 -L8k Vr = 2.0km/s k

Vr=m/s km/s AYG.Z ck c=9 km/s f=135I

=11 135° ITI X 0.3.

.,./gi A

I I

~pl Ip2 2

MAJO.Z

_ _ _ _ _ _..f=

306o A r 0

90 180 270 360 Azimuth from North (0)

Figure 4.3. The measured time difference between the first peak and second peak of teleseismic P-wave displacement seismograms versus source to receiver azimuth. On average we observe a separation of 3 sec for the two peaks that we identify as major subevents but for azimuths to the southeast there is a merging of the two pulses, which are nearly separated for paths to the northwest. We fit the time difference using the directivity equation for the pulse width based on the Doppler effect for a range of fault lengths and rupture velocities to estimate the best fitting fault strike or rupture direction. We fit the data using both minimum deviation Ll-norm (solid line) and least squares (dashed line). The error curve for the best rupture direction is at 1350 using the LI-norm and at 1150 using the L2-norm.

With some of the earlier observations (e.g., Figure 4.3), we made some initial assumptions in computing a grid of Green's functions (Earth wave propagation function from an impulse point source).

A grid of point source locations is computed based on the extent of the aftershock locations. We select a range of azimuths from 900 to 1600 (or fault plane strike=27 0 0 to 3400) for the line source strike with twelve point sources distributed in 2 km increments. We also try a series of hypocenter depths between 4 and 18 km in I or 2 km depth increments. We also use a range of depths between +/- 1 km increments 5 km above and below the hypocenter depth to compute the Green's functions for the subsequent subevents.

We assume only two subevents with triangle shape duration of 4 seconds (2-sec rise and 2-sec fall-off).

The second subevent can occur between 0 and 5 seconds after the initiation of rupture. We use the 1-D layered velocity model Central-California (CC) from McLaren and Savage (2001) listed in Table 4.1 for the source region and a generic iD model for the receiver sites (Table 4.2). Use of the calibrated models described in Section 3, which were developed after this teleseismic analysis, in place of the CC model would make no practical difference to our results because the models are similar and because teleseismic P-waves are not very sensitive to the crustal structure at the source. The Green's functions are computed using the Haskell propagator matrix method described by Bouchon (1976) and Haskell (1960; 1962). We account for attenuation using a t* of 1 sec for P-waves and t* of 4 sec for S-waves (Langston and Helmberger, 1975). The attenuation operator t* = t/Q, where t is travel time and Q is quality factor, integrated over the travel path.

A-28 PG&E Letter DCL-05-002 TABLE 4.1. Source Region Velocity Model VP V,

Density Thickness (km/s)

(km/s)

(g/cc)

(km) 4.00 2.31 2.05 4.0 5.70 3.29 2.59 6.0 6.00 3.46 2.69 4.0 6.20 3.58 2.75 10.0 8.00 4.62 3.33 0.0 TABLE 4.2. Receiver Region Velocity Model VP Vs Density Thickness (km/s)

(km/s)

(g/cc)

(km) 5.80 3.35 2.63 5.0 6.50 3.75 2.85 15.0 7.80 4.50 3.27 15.0 8.00 4.60 3.33 0.0 We performed a grid search over centroid depth, rupture velocity, rupture azimuth, and time shift. The rupture velocity constrains the time and location of the subevents. We estimate high rupture velocities of about 3.54 km/s (usually Vr, 0.8*V, for crustal earthquakes).

This estimate is not unreasonable considering the rupture needs to travel the 25-30 km length of the fault suggested by the 1 -day aftershock locations in about 10-12 sec. Other studies also indicate that there is a region of high slip 10-15 km from the epicenter, which is well located by the USGS with 1 station within 5 km and 10 stations within 100 km with an azimuthal gap < 90°. This indicates a rupture velocity of about 3-3.3 km/s, much higher than the preliminary estimate of 2.1 km/s suggested by D. Dreger in Hardebeck et al. (2004). Figure 4.4 (left) shows the centroid depth with the lowest error at 12 lan, 4-5 km deeper than the 7-8 km hypocentral depth estimated by the USGS using P-wave travel-times and double-difference methods. Figure 4.4 (right) shows the error space for rupture velocity and time shift at the centroid depth iteration at 12 km depth and a rupture azimuth of 1300. The contours show that rupture velocity is poorly resolved, which is expected using only teleseismic body-waves, but a range between 2.75 and 4 kml/s gives slight and insignificant improvements in waveform fitting., Since absolute timing is not used in the inversion of the teleseismic data (because the velocity model only includes the crust), a simple shift in time is applied to all the data so they are aligned with the synthetics. This alignment is performed before the inversion.

Figure 4.5 shows the resulting spatial distribution of the subevents, the mechanisms, and source time function for the centroid depth of 12 km (Table 4.3). The USGS epicenter location is assumed but not the hypocenter depth. The first subevent is located near the hypocenter indicating that a significant amount of seismic moment was released near this location, and is estimated to have an uncertainty in strike of about I0°. The second subevent occurred 3.2 sec later about 12 km to the southeast at 1100 azimuth from the hypocenter and at a depth of about 10 km. The strike of the second subevent is rotated clockwise by about 20-30°, with an uncertainty of about 100. This result is confirmed in double-difference relocations and geologic mapping indicating some structural complexity or a jog or bend in the fault. The moment release of the second subevent occurred between 3.2 and 7.2 seconds and was about 1.5 times larger than the first subevent. The surface projection of the focal planes from the point source locations to the surface suggests that rupture may have occurred along the Oceanic fault zone, although the surface projection lies west of the mapped fault zone. The surface projection of the plane of the second event has a more northerly strike than the mapped southern part of the Oceanic fault, but the difference in strike is within the uncertainty in strike of the second event.

A-29 PG&E Letter DCL-05-002 0.395 a:

o 0.42 -

a:

w 0

Z 0.40 -

N-J1 0,

S,

/

• ,*- 0

-0

/

4-Co 2-50 270 280 290 300 310 320 STRIKE (0) 330 340 I. I I

I III I.

Ij I

I I

'I III IDI 0 2 4 6 8 101214 tshift (sec) 3 0.50-a:

0 W 0.45 -

E 0

0.40 -

I I

I 3.314 3.8/13 3.8/8

-0 "4.0/1

,~~

"4 II 0.4 0.5 0.6 L2-norm Error 4

6 8

10 12 14 DEPTH (km) 16 18 Figure 4.4. (LEFT) Depth-error space. We performed a grid search over centroid depths between 4 and 18 km at 2 or 1 km increments and rupture direction azimuths from 110°-160°. The numbers labeled above the focal mechanisms are the best rupture velocity and time shift for that depth. The best centroid depth is 12 km but may only be resolved to within 10 and 14 km. The best fitting rupture azimuth is between 2900 and 3100 (or 1100 and 1300)

(RIGHT) At each depth step, we also perform a grid search for the rupture velocity (Vr) and time shift (ta,)f. The rupture velocity is used to determine the space and time distribution of point sources and the time shift allows for a shift between all the data and synthetics/Green's functions for the best alignment.

The error space for rupture velocity and thift at a centroid depth of 12 km is shown.

TABLE 4.3. Source Parameters Subevent Time R

(p (°)

Z MK Strike Strike Dip Dip Rake Rake (sec)

(km)

(km) x1025 P1 ()

P2 (0)

P1 P2 P1 (0) P2 (0) dyne (0)

(0) cm 1

0 0

110 12 1.64 138 298 35 57 106 79 2

3.2 12 n/a 10 2.48 116 320 46 47 72 108 TOTAL 3.97 124 312 40 50 84 95 PI-Nodal Plane #1, P2-Nodal Plane #2, Time-onset time of subevent rupture.

P2 is the preferred fault plane Z-centroid depth, (9-rupture direction, R-rupture length, A-30 PG&E Letter DCL-05-002 0

9-O a) 0 (n -A O#2 5

Q

.I. I 0

12 24 Distance (km) strike=2900

1. 1

1. 2 Total

.1 36 35.8 35.6 35.4

-121.4

-121.2

-121

-120.8

-120.6 Figure 4.5.

Map of the subevent locations and mechanisms. The first subevent was assumed based on the network locations using P-wave travel-times. The depth of the first subevent is at 12 km. The second subevent is 12 km from the first subevent at an azimuth of 1100 to the south-east and at a depth of 10 km.

The first and second subevents have seismic moments of 1.61 and 2.42x 1025 dyne-cm respectively. The second subevent occurred 3.2 seconds after the initiation of rupture. The heavy dashed lines represent the surface projections of the fault planes based on the location of the subevents and their individual focal mechanisms assuming rupture on the northeast dipping nodal planes. Also shown are some of the active fault zones. The change in strike of the second subevent is consistent with the bend of the Oceanic fault zone.

A-31 Con PG&E Letter DCL-05-002 P1V 1143 M2 I

nTE ua A a1 7A 3A Ji U2 127.3 J\\X JVV n.1 31 II5 MY S12 101 25l 4

m7m tA M

14 115 A

CP4 1tJ A

IP!

M1U P

4.s LIS LIS P1

.T l

1

.7a P

synheics for th ulilepit orc oenenos A-12 77133 13 3

Pr p

17 n VS synheic)frth mutpl pon orcemmntnos 23.3 12 2A1-3 S

C4t

Enclosure I PG&E Letter DCL-05-002 1IT W13 Jv IT p U IIMe 1W5X au rt rs k* n1 Ptt 77.1 SjA U

7A i aro A LN. U"

£25 k31 JaR.

1i 138 SITW A

-J I

3 V t To t

• A UU STU

'it 5t J j Y 5 Von 1341 Ino IOU VW

, 7 uhns e

gP a C~ t I-"

JTIR&

Jl= I 1179

• LI j 1

INA 3 #8M m9 au2 1i mu re

r. mi rr P 2 N

P51 Figure 4.6. Continued.

The Kikuchi and Kanamori method is limited by some basic assumptions made to simplify the rupture model. One is that a finite number (two in our case) of discrete subevents is forced to fit the source time function of the entire rupture.

This assumption causes underestimation of the total seismic moment because we are only fitting the two major episodes of energy release and missing the more continuous and detailed features of the slip. The lower seismic moment and moment magnitude (M, 6.3 instead of M, 6.5-6.6) is an artifact of this limitation. To get a better estimate of slip and moment we use a finite fault inversion described in Section 5.

4.5 Discussion and Conclusions The San Simeon earthquake initiated at a hypocentral depth in the range of 10 to 14 km, with a best estimate of 12 km depth. This initial event was followed 3.2 seconds later by a second event with a depth of 10 km located about 12 km to the southeast having a strike that is rotated 20 degrees clockwise with respect to the initial event. The location and orientation of the initial event is consistent with it having occurred on the Oceanic fault, and the location, orientation and change in strike of the second event are consistent with the change in strike of the Oceanic fault to the southeast of the epicenter. These results were obtained using only teleseismic waveform data, independently of geological and geodetic data.

A-33 PG&E Letter DCL-05-002

5. FINITE-FAULT KINEMATIC SLIP MODEL 5.1 Introduction In this Section we present a more comprehensive examination of the 2003 San Simeon earthquake rupture process using all of the most significant data including teleseismic P-waves, locally recorded strong ground motions, and horizontal coseismic displacements.

We use all available seismic and geological information about the mainshock faulting and aftershocks to build a fault model comprised of 147 subfaults having 2 by 2 km dimensions (Figure 5.1). The Green's functions which represent the impulse response of the Earth at each subfault location are calculated using the 1-D layered Earth velocity models that best predict the ground motions of earthquakes in this area up to frequencies of at least 0.2 Hz to about 0.5 Hz. The main objective of this portion of the study is to improve our understanding of the rupture process that will allow for a more detailed analysis in near-source ground motions and stress changes from this and possibly future earthquakes in this region. In Section 6, we compare our kinematic slip model with preliminary models by D. Dreger and C. Ji (see Hardebeck et al 2004, CISN report and Ji et al. 2004), and discuss the limitations of the various models.

5.2 Source Inversion Methodology and Multiple Time Window Inversion Technique Seismic and geodetic data are used to invert the representation theorem for the spatial and temporal evolution of slip and rake over the fault. The 3 component seismogram U(J is related to the space and time integration of slip distributed across the fault Z where, Efl UnUx,t= El

' 9j(e,r)

Gni~j(3E.e;t-r)] [it (I) is the nht component of observed displacement. The term Aijuj(S:,r) is the product of the rigidity, fault orientation unit vector, and slip function (fault slip) at point 4 and time I. The term Gfl (x,,,t -T)is the Green's function that describes the wave propagation from each point on the fault to the receiver. The vector x is the relative location of the source and receiver. The 4: and X are spatial and temporal variables of integration and i and j are orientation indices.

Our method for determining the slip function uj({,r) is explained in more detail by Thio et al (2004a,b) and is similar to the multiple time window method of Hartzell and Heaton (1983). The fault plane is described by a grid of points representing small subfaults for which Green's functions are computed at each grid point for each observation point. A propagating slip band is imposed from the hypocenter and propagates at a fixed number of time steps. This slip band is characterized by two parameters, maximum rupture velocity and rise time. Grid points that are contained within the slip band, at each time step, are combined and cast into a set of normal equations of the form Ax = b, where A contains the Green's functions from every grid point to every station, x is the solution vector containing the slip value at every grid point, and b is the vector containing all the data. The normal equation is solved simultaneously using a least squares solver with a positivity constraint to prohibit backslip (Lawson and Hanson, 1974).

Regularization and spatial smoothing are applied to prevent wildly oscillating solutions. This is achieved by imposing a smoothing condition in which the amplitude of slip at one point is forced to be equal to its neighbors. Lower smoothing will lead to models with varying complexity and large values of smoothing will lead to models with uniform slip.

Thio et al. (2004a,b) allows for variable rake, in which case every original rake is split into two conjugate vectors where the new rake angles are different from the original by 445' so that the same positivity constraint can be used.

The final rake vector is the vector summation of the two individual rake vectors. Allowing for variable rake thus leads to a doubling of the number of unknowns, so we need to regularize the two slip vectors to prevent strong, unphysical variation in rake. This is achieved in a manner analogous to the spatial smoothing, where instead of forcing the amplitude of a slip vector to be equal to its neighbors, we force it to be equal to its conjugate slip vector. In this case, a strong smoothing constraint will cause the two conjugate vectors to become equal, and thus yield the orientation of the A-34 PG&E Letter DCL-05-002 initial single rake vector.

We imposed a weighting scheme that uses overall weight per data type, and within each group, an automatic weighting based on the absolute amplitude and number of points (e.g., Thio et al., 2004a,b and Ichinose et al., 2003). This was made to prevent the contribution of one dataset from dominating the solution especially in the case of greatly varying amplitudes and sample lengths. Many of inversion weighting parameters including those for slip and rake smoothing were selected on the basis of trial and error.

Subevent

1. 1 z Subevent 36 35.8 z01-0 4S

-J 35.6 35.4 121.4 121.2 121 120.8 120.6 120.4 Longitude (OW)

Figure 5.1. Map of CSMIP, BDSN, and TRINET strong motion sites and GPS sites used in the inversion of the San Simeon Earthquake. We constructed a fault model using the focal mechanism solution from the multiple point source moment tensor solution (Figure 4.5) obtained using only teleseismic data (solid nodal planes). The epicenter of subevent #1 was adjusted so that its surface projection coincides with the Oceanic fault. The upper (southwest) edges of the two fault planes are at depths of 1 km.

5.3 Strong Motion Data, Global Positioning System Data and Velocity Models We gathered strong motion data from California Geological Survey (CGS) California Strong Motion Instrumentation Project (CSMIP) sites Cambria Hwy 1 bridge (CAMB) and Templeton Hospital Ist Floor A-35 PG&E Letter DCL-05-002 (TEMP), and the toe of the San Antonio Dam. Theses sites operated SMA type accelerometers recording in cmxs-? and integrated to velocity and displacement. We also included data from the Berkeley Digital Seismic Network including accelerometers at Parkfield (PKD) and San Andreas California (SAO). We also included data from TRINET operated by Caltech and USGS. These sites included Park Hill (PHL) and Simmler (SMM). These records were instrument corrected to cmxs-2 and integrated to velocity and displacement. We mainly focus on records from the four nearest stations CAMB, TEMP, PHL, SADM, and PKD (Figure 5.1) which have the most importance in the distribution and timing of slip along the fault. Global Positioning System sites were operated by the Southern California Integrated Geodetic Network (SCIGN). The current best estimates of coseismic displacements were provided by Michael Heflin (Jet Propulsion Laboratory) by taking the difference in average positions over 24 hours2.777778e-4 days <br />0.00667 hours <br />3.968254e-5 weeks <br />9.132e-6 months <br /> between December 21 and December 23. Other solutions obtained using different strategies varied by 1 to 8 mm.

We included only the three nearest and most significant sites RNCH, CRBT, and LOWS (Figure 5.1).

Teleseismic data and station descriptions were discussed in Section 4 and shown in Figure 4.1.

Green's functions are computed using the'Haskell propagator matrix method described by Bouchon (1976) and Haskell (1960; 1962). We account for attenuation using a t* of 1 sec for P-waves and 4 sec for S-waves (Langston and Helmberger, 1975). The teleseismic Green's functions were computed using a calibrated 1-D layered Earth velocity model for the source region and a generic model for the receiver sites. Green's functions for the local accelerometer and GPS sites were computed from layered Earth models determined using aftershocks in Section 3. We used the PHL model (Table 3.5) for seismograph sites and the SAO model (Table 3.2) for GPS sites. The Green's functions were computed using an f-k reflectivity summation method (Zeng and Anderson, 1995). The data and Green's functions were band pass filtered using a 3 pole 2 pass Butterworth filter between 0.01 and 0.7 Hz.

5.4. Finite Fault Model The fault geometry is based on the focal mechanisms obtained from inverting teleseismic data using the multiple point source moment tensor inversion (see Figures 5.1 and 4.5) that resolved two subevents.

The selection of the two nodal planes was based on the northeast dipping aftershock distribution listed by the PG&E, Hardebeck et al. (2004), Hauksson et al. (2004) and Advanced National Seismic System (ANSS) catalogs. The second subevent has a northeast dipping nodal plane (strike=3200, dip=470) with a strike that differed by 200 from the same northeast dipping nodal plane of the first subevent (strike=2980, dip=570). Since the second subevent occurred about 10-15 km from the hypocenter (assuming a rupture velocity betveen 3 and 3.5 km/s), we assumed that the first fault segment has a strike length of 10 km and a dip width of 18 km. The second fault segment has a length of 22 km and a width of 20 km. This change in strike is corroborated by the change in strike of the surface trace of the Oceanic fault and details in the aftershock locations. The subfault dimensions are 2 by 2 km. The upper (southwest) edges of the two fault planes are at depths of 1 km. The hypocenter is based partly on the PG&E relocation using a more appropriate local velocity model for this region. We modify this epicenter by shifting it horizontally by 1.5 km to the northeast to make the 11 km hypocenter depth consistent with the 570 dipping fault plane projected to the mapped surface trace of the Oceanic fault. The PG&E location is shifted to the northeast by I km from the location listed in the ANSS catalog (i.e., California Integrated Seismic Network hypocenter location). An additional shift of 1.5 km makes it still within the 2 km uncertainty assumed in local network locations especially when travel time information is used from farther stations.

We use a maximum rupture velocity of 3.5 km/s, a maximum rise time of 1.5 sec, and a minimum rise time of 0.5 sec. We use 18 time windows with 0.5 sec increments. We assume a triangle shaped source time function for each subfault with a rise of 0.5 sec and a fall-off of 0.5 sec. A rupture front and a healing front originate from the hypocenter and their subsequent locations are computed along the fault plane for each time step "time window" based on the maximum rise time and rupture velocity.

In this case, slip is only allowed within this 5.25 km wide slipping band limited by the top and bottom depth extent of the fault plane. Variable rake is allowed between 450 and 1350.

A-36 PG&E Letter DCL-05-002 22 December 2003 San Simeon Earthquake (Mw 6.58) 36 Z

0._

te-(U

-J 35.8 35.6 35.4 Longitude (OW)

Figure 5.2. The cumulative slip and rake distribution of the 2003 San Simeon earthquake estimated from the inversion of teleseismic P-waves, locally recorded strong motion velocity and displacement seismograms, and coseismic horizontal vector displacements measured from SCIGN GPS instruments.

The first week of aftershocks are plotted to show that they occur mostly along portions of the fault that did not rupture or ruptured with lower slip. The rake angles "slip vectors" shows the direction of the hanging wall relative to the foot wall.

5.5. Kinematic Slip Model Results Figure 5.2 shows the map view of the slip and rake distribution along the two fault segments of the Oceanic fault during the 2003 San Simeon earthquake from the inversion of vertical component teleseismic P-waves from 46 global stations, 3 component strong motion records from 5 local sites, and horizontal displacements from 4 sites recorded by GPS instruments. The slip is concentrated mostly near and below the hypocenter at 11 km depth along the first fault segment. The slip is less concentrated along the second fault segment but is highest along the edge or bend shared by the two segments. One scenario is that the slip of the first segment triggered the second subevent by decreasing the fault normal stresses along the second fault segment. The slip velocity function shown in Figure 5.3 indicates that slip was rather uniform except for a sharp increase in slip at around 4 and 6 sec. The overall shape of this function A-37 PG&E Letter DCL-05-002 can be represented with two triangles with duration of 4 sec each centered on 2 sec and 5 sec similar to that determined in the teleseismic analysis (Figure 4.5). One large slip region is coincident with the hypocenter and the other is located about 12 km southeast along strike from the hypocenter at about the same depth (about 2 km up-dip) consistent with previous results and studies.

4 0,728 0

=Q 0.000 0

10 seconds Figure 5.3. Slip velocity source time function.

Figure 5.4 shows the time progression of slip along the two fault segments at 0.5 sec increments in 17 time windows starting from the second time window (0.5 to 1.0 sec). We removed the first time window from the plot since it was very similar to the second time window.

The first subevent or asperity is consistent with the previous results in that it occurred near the hypocenter. The peak slip here is about 3.4

m. The second subevent ruptured on the second fault segment which began rupturing at about 3.5-4 sec into the rupture. The final seismic moment is 9.4x 1025 dyne-cm (Mw 6.59). This moment is consistent with our very long-period point source moment tensor solution using regional-waves (M,, 6.55). This indicates that the earthquake may have been slightly larger than initially estimated. The rupture area is 572 km2 and the average slip is 0.35 m. The rupture area scales with the relationship of Somerville et al.

[1999] and is within the scatter of global crustal earthquakes in tectonically active regions. The average slip is low compared to global averages, but the average asperity slip of 0.63 m scales with other global earthquakes. The peak slip is 3.4 m located near the hypocenter and the largest asperity has an average slip of 1.55 m. Rupture velocities ranged from 0 to 3.5 km/s and the average was 3 kln/s. The source parameters are summarized in Table 5.1.

Table 5.1. 2003 San Simeon source parameters from the kinematic slip model inversion Segment #1 Segment #2 Total Area 164 km2 408 km2 572 km2 Average Slip 0.33 m 0.33 m 0.35 m Peak, Slip 3.4 m 1.4 m 3.4 m Average Rake 750 1050 Seismic Moment (dyne cm) 2.52x1O2' 6.89x1025 9.41x10'5 Asperity(l) Area 232 km2 Average Asperity0l) Slip 0.63 m Rise time (minimum/maximum) 0.5/1.5 sec Rupture Velocity (minimum/maximum) 0-3.5 km/s

)Asperity Area is defined here as subfaults with slip greater than the overall average slip.

Figure 5.4. (next page). Map view of the time progression of slip and rake in 0.5 sec increments. The first time window (0.0 to 0.5 sec) is not shown and is similar to time window 2 (0.05-1.0 sec). Rake or slip vectors indicate hanging wall motion relative to the foot wall, and are scaled relative to the subfault slip value. The star marks the location of the hypocenter at about 11 km depth.

A-38 PG&E Letter DCL-05-002 A-39 PG&E Letter DCL-05-002 Figure 5.5 shows the fit to the strong motion data at 5 stations nearest to the earthquake. All recorded the earthquake on scale without clipping and the data and Green's functions were filtered between 0.01 and 0.7 Hz. The misfits at stations SADM and TEMP (Figure 5.5) would probably improve with the use of higher frequency Greens functions, as we cannot fully fit the TEMP higher frequency pulse and amplitude given the high frequency limits of the Green's function.

In addition, the misfits in displacements at CAMB, SADM, and TEMP, relative to PHL and PKD reflect the loss of lower frequencies from inadequacies of the instrumentation and during data processing of strong motion instruments. The observed and predicted telseismic P-wave displacements (Figure 5.6) correlate well in amplitude and phase. The addition of velocities (Figure 5.7) adds a higher frequency range and improves the resolution of slip distribution relative to using only the displacements.

The addition of higher frequencies may introduce artifacts as velocity model uncertainties can map into the slip model solution.

The results here show that both displacements and velocities fit equally well, so the artifacts are minimized especially with the use of a calibrated velocity model for the source region. The disadvantage of using only teleseismic data is that it cannot resolve variations in rupture velocity and rise time. This was shown in the modeling in Section 4. The addition of local strong motion data and horizontal coseismic displacements has improved the uniqueness of the solution and provided constraints on the higher rupture velocities resolved in the complete kinematic slip model solution.

PKD D = 64.0 km f = 730 PHL D=60.1 km f = 1230 CAMB D=12.7 km f =189° SADM D=22.3 km f =630 TEMP D =38.5 km f =1160 dispIaorent vodoaty

_ F ZAIV4

-4r di spl acement N-S --

W'L Veloci t N-S --

ed

-X

j
I edata At sy di SOlaceent E-W velocity E-W UP*

IEI

.; 3,

V1\\^

Ah' 20 -

50 seconds 2

0

_R 3 0 0

Figure 5.5. Observed and predicted displacement and velocity seismograms. Note that the scale for the records at TEMP are multiplied by 4x relative to the scale used for the other four stations.

A-40 ci
-

PG&E Letter DCL-05-002 p1 p2 AFI BRVK 10.2 4.3 ~YJ 68.8 IIf 91.0 23435 KBS 9.4 62.4 9

MAJO 6.8 77.1 306 ALEA 9.3~Ij 50.8 8

ANTO 1.6 101 0 20 'A/

ARU R

5.7 88.2 1

8 f KDIAK MDJ A~

8.1 flI 4

5.0 30.4 J~l 1 V 77.6 V' 1 326 316V

(

OBNA h 5.9 '1V 87.5 jvL 13

)TAV 22.3 53.1 122

)AYG 15.5 46.2 135 PET 6.8JI 56.9 Yld 315 PTC 24.9 61.0 189 RAR 19.0 67.4 219 RCBR 6.2 89.5 97 RPN 14.3 63.5 168 SJG 14.0 51.4 95 SNZO 3.1 96.6 223 SSE 3.1 91.5 311 I1 A

I f ESK 7.4 7532° TEIG LMA 14.2 32.71/f 109 TRQA 5.5 836303 /I F

PV fX/I\\

LwII 2.3 96.6 jv 339 2.1 J, 96.3 67.3 IVIY/

331 YA 6

BJT A

88.0 JVVV 320 HIA 6.7 79.7 324 INCN 5.3 83.9 312 JTS 12.7 41.4 119 KONO 7.9 76.6 23 BORG A~

5.1 4 iIV lJ 62.3 A 29\\

l /\\

3 1KURKA

~~3.7

.1 92i2 fij, MA2 6.9 324 PMG 2.5 96.8 263 -'V POHA 11.7 34.2 252 TATO data 2.4.AjgA.Idat 95.5 AV\\

synthetic 0

2

°25 seconds Figure 5.6.

Observed and predicted teleseisrmic P-wave displacements.

Amplitudes are in units of microns and labeled below the station code. The epicentral distance in degrees and source to receiver azimuth are also labeled below the station code. The station locations are shown in Figure 4.1.

A-41 C54

PG&E Letter DCL-05-002 AFRAA COLAA~A 7.2 3 -MJW 68.8 3 31I.73 234 340V~

ALE5 DWPF. Al 50~.8 1114I1~f 9 1 KIP k

.OTAV 11 1

,l 19.7 35.2 t 531 11 256 122 f SNZO 3.9 96.6 223 V

SSE 3.6 Iw 91.5 31 1 Yss 7.3 68.7 iUd\\

313 YAK jAl~1 data 5.1 i 1 J2 67.3 syntheic 331 1' 0

25 seconds KIV 1.2 -Ae-99.4 12 PAYG 14.5 kfJA 46.2 UlM.

135 ANTO 1.1 101.0 20 -e ENH KONO PET ;

990 76.6 56.9 J 317 'W~~23 315 TATO 2.3

• 1 95.5 306 ARU. ltliA HIA J~,KURK V..A.

PMG 5.4~i'~

4.3 1IP'&'2.3

-,V/sF, BBSR I

7.2 jUf 77 BFO i 85.0A'IM.

31 1 BJT 5.0 i 88.0 A'Ayrt'.b 320 79.7 f J tIL/ 'VL--

"U.o -V 324 t

348 263 INCN ASI MIA2 PTCN 6.5

7. 3 ULL'1.

83.9 t

59.1 61.0 312 324 189 JTSAIAMAJO RCBR aI 11.6 51A 5.3 6.5 41.4 77.1 89.5 flp.

119 306 97 KBS MDJ RPN 5.0 3.8 1 5.2 62.4 77.6 J~IfN,63.5 mu 9

316 168 ULN I 6.1 I~VV 86.3 330 WMQ 2.2 JIl 96.6 'VIr--

339 XAN 2.0 96.3 320 BRVK LhIAA.OBN

.1 91.0 4.6 'V4N/

6.7 91.0 r"N.

72.3 2i.,

87.5 35 11 13 SJG j 12.8 51.4 /

95 Figure 5.7.

Observed and predicted teleseismic P-wave velocities.

Amplitudes are in units of microns/sec and labeled below the station code. The epicentral distance in degrees and source to receiver azimuth are also labeled below the station code. The station locations are shown in Figure 4.1.

A-42 C 5S PG&E Letter DCL-05-002 5.6. Discussion Using Green's functions calculated from the calibrated regional ID seismic velocity models developed in Section 3, we proceeded to develop a finite fault rupture model of the earthquake (Figure 5.2) using strong motion recordings, regional and teleseismic body waves and surface waves, and regional geodetic data. This rupture model describes the location of the hypocenter, the orientation of the fault rupture planes, and the distribution of slip on the fault.

We developed a rupture geometry model of the earthquake that consists of two fault segments having different fault strikes, based on the teleseismic model of Figure 4.2. This change in the strike is indicated by three independent lines of evidence. First, the multiple point source model derived from teleseismic data indicated the presence of two subevents having different strike angles, as described in Figure 4.2.

Second, there is a change in strike of the aftershock zone. Hardebeck et al. (2004) describes this change in strike using two fault planes having strike angles of 295 and 330 degrees from double difference aftershock relocations, similar to the strike angles of 298 and 320 degrees that we obtained. Third, there is a change in the strike of the mapped surface faults in this region. Based on this evidence that the earthquake most likely occurred on the Oceanic fault system, we adjusted the epicenter location slightly so that the surface projection of the fault plane of the first subevent coincides with the mapped surface trace of the Oceanic fault.

The finite fault rupture model is shown in Figures 5.2 and 5.4. Rupture on the first fault segment, which contains the hypocenter, was mostly downdip of the hypocenter and occurred in the depth range of 12 to 15 km. This is consistent with the first point source derived from the teleseismic data alone.

Rupture of the second fault segment was largest on an asperity located about 12 km east-southeast of the hypocenter, at a depth of about 9 km. This is consistent with the second point source derived from the teleseismic data alone. Unlike the first fault segment, whose slip is concentrated near the hypocenter, the second fault segment had a fairly uniform amount of slip.

The rupture model has rupture velocities between 3 and 3.5 km/sec, which is approaching the local shear wave velocity of the upper crust. This high rupture velocity is physically plausible for along-strike rupture propagation of dip-slip faulting, because it constitutes a Type III crack (Scholz, 1990). For along-strike rupture propagation on strike-slip faults (Type II crack), the rupture velocity is bounded by the Rayleigh wave velocity, which is lower than the shear wave velocity.

The high rupture velocity contributed to strong rupture directivity effects, described further in the Section 6.

A-43 PG&E Letter DCL-05-002

6. COMPARISON WITH OTHER RUPTURE MODELS We compared the rupture models obtained in this study with rupture models derived by other investigators, including Hardebeck et al. (2004), Dreger (CISN, 2004) and Ji et al. (2004) (Table 6.1).

Differences among these models include the depth of the hypocenter, the depth range of faulting, and the presence of a change in fault strike southeast of the epicenter.

Hardebeck et al. (2004) contains two preliminary rupture models. Both rupture models assumed a hypocentral depth of 7.8 km, which is based on a crustal velocity model for northern California that may not be appropriate for the central California coast. A common feature of these models is a region of large slip (asperity) located about 15 km southeast of the epicenter. What differs among the models is the depth of the asperity. The rupture model derived by Ji using teleseismic data has an asperity at a depth of about 10 lan. This is consistent with the result that we obtain using teleseismic data; we find that the second subevent is located at a depth of about 10 km (Figure 4.5 and Table 4.3). In both of these studies, the depth is controlled by teleseismic depth phases.

When Ji et al. (2004) subsequently included geodetic and strong motion data with the teleseismic data, they obtained a depth of about 5 km for the asperity. Similarly, the rupture model derived by Dreger and described in Hardebeck et al. (2004) using regional broadband seismograms has an asperity depth of about 6 km. When we included the strong motion and regional data, the depth of the asperity on the second fault segment was 9 km, similar to the 10 km depth obtained from the teleseismic data alone. We attribute this difference in depth to our use of a different crustal velocity model, which was calibrated against regional recordings of the larger aftershocks as described in Section 3.

Ji et al. (2004) examined rupture models having two fault segments with different strikes, but concluded that a single fault plane model provided a better fit to the displacements at the GPS stations.

However, we found that the two fault segment model provided better waveform fits to regional broadband seismograms recorded east of the fault, while obtaining a reasonably good fit to the GPS displacements, compared with a single fault segment rupture model that we developed. Accordingly, we prefer the two fault segment rupture model.

The Templeton (TEMP) recording of the San Simeon earthquake has large strong motion amplitudes (horizontal PGV - 20 cm/s), which have been attributed by Hardebeck et al. (2004) to rupture directivity effects. Dreger (CISN, 2004) did not include Templeton in his rupture model inversion, but compared the recorded waveform with that calculated from his rupture model. He obtained a good fit in both amplitude and waveform, after scaling the synthetic seismograms up by a factor of 1.62 to account for the difference in near shear wave velocity between the synthetic seismograms (1.5 km/sec) and the Templeton site (0.4 km/sec).

Ji et al. (2004) also obtained a good fit to the waveform recorded at Templeton. Their model assumed a rupture velocity of 3 km/sec, and allowed the rupture time to vary in each fault element, as shown by the rupture time contours in Figure 3 of their paper. These contours indicate that the rupture velocity in the horizontal direction through the asperity was about 3.5 km/sec. This high rupture velocity, close to the shear wave velocity, produces large ground motion amplitudes, as was required by Ji et al. (2004) to match the recording at Templeton.

Our rupture model is based a maximum rupture velocity of 3.5 km/sec, similar to that of Ji et al.

(2004). We down weighted the contribution of the Templeton recording in our rupture model inversion, because this record contains high frequency components that are not modeled well by our choice of fault grid size (2 km x 2 km). The synthetic seismograms from our rupture inversion underpredict the large amplitude motions recorded at Templeton, as seen in Figure 5.5.

A-44

Enclosure I PG&E Letter DCL-05-002 Table 6.1. Summary of Fault Rupture Models of the San Simeon Earthquake Author Plane Strike Dip Length Width Top Bottom Seismic M,

(0 E

(°)

(km)

(km)

Depth Depth Moment of N)

(km)

(km)

(x lOe dyne.cm)

Dreger Ji URS 290 58 44 292 52 45 1

298 57 10 2

320 47 22 Total 22 0

21 0

18 1.0 20 1.0 17.4 15.5 16.0 15.6 5.7 6.2 6.5 6.5 9.4 6.6 A-45 PG&E Letter DCL-05-002

7. RUPTURE DIRECTIVITY EFFECTS IN STRONG GROUND MOTIONS We performed strong motion simulations using the crustal velocity model derived in Section 3 and the fault rupture model derived in Section 5 to simulate broadband (0-10 Hz) ground motions around the fault, using the procedure described by Graves and Pitarka (2004). This procedure takes account of site conditions in an approximate manner using a site conditions map of California based on geology and shear-wave velocity (Wills et al., 2000). The motions were simulated on a 2 km by 2 km grid of sites as shown in Figure 7.1. The white regions of the map have NEHRP soil categories of A through C, corresponding to shear wave velocities above 360 m/sec, the brown regions have C/D site conditions corresponding to shear wave velocities of 270 - 555 m/sec, and the yellow regions have D site conditions corresponding to shear wave velocities of 180 - 360 m/sec.

..X'.

X

.V..X....

36-lll a

X X

.i k

000

\\;,0 SADM........

\\'\\

35.8

=,,

~~~~.

, d Ad SA 35.6 35.4 35.2

-121.4

-121.2

-121

-120.8

-120.6

-120.4 Figure 7.1. The two segment fault model of the 2003 San Simeon Earthquake used in the simulation of ground motions on a 2 km by 2 km grid of sites. The solid lines represent the top edges of the fault plane at depths of 1 km. The hypocenter is marked with a star. The model is described in Section 5, Figures 5.1 and 5.2.

From the simulated time histories, we measured peak acceleration and peak velocity and these values are shown in map form in Figures 7.2 and 7.3, respectively. The largest simulated accelerations are near 1 g and occur just east of fault segment 1, over the northern portion of segment 2, and in the region between the faults where the site category is C/D (denoted by brown shading, which has stronger amplification than the surrounding region with site category C to A, denoted by white). A smoother peak acceleration map would be obtained if we were to average the ground motions over a suite of earthquake A-46 PG&E Letter DCL-05-002 scenarios, because each scenario is influenced by stochastic aspects of the simulation procedure. The largest simulated velocities are near 60 cm/s and occur on C/D site conditions in a zone south and east of fault segment 1, primarily over the central portion of fault segment 2.

The largest slip on the fault is at the hypocenter, which is located along the northern edge of fault segment 1. The regions of largest simulated acceleration and velocity occur southeast of the epicenter, suggesting that rupture directivity is influencing the distribution of simulated ground shaking.

Additionally, the peak velocities (Figure 7.3) show more southeastward extension than the peak accelerations (Figure 7.2), consistent with the observation that rupture directivity effects tend to be stronger at longer periods.

The strongest directivity effects in the simulated motions do not persist more than about 20-25 km southeast of the epicenter. Thus, our simulations underpredict the large amplitude motions observed at Templeton (Figure 7.4), which have been attributed to rupture directivity effects by Hardebeck et al.

(2004X Since most of the slip in our current rupture model is located on the deeper part of the fault, and within about 15 km of the hypocenter, it is difficult with this model to generate strong directivity effects extending to the distance of the Templeton station (about 40 km from the epicenter). We conclude that aspects of ground motion generation and propagation other than those included in our ground motion simulation procedure may be responsible for the large ground motions recorded at Templeton. These may include departures from 1D seismic velocity structure in the crust, and the presence of basin effects or unusual site amplification effects at Templeton.

36 - no 1

\\\\

35.8 35.6 -

35.4-35.2

-121.4

-121.2

-121

-120.8

-120.6

-120.4 Figure 7.2. Map of Simulated Peak Ground Acceleration (PGA).

A-47 PG&E Letter DCL-05-002 36 35.8 35.6 35.4 35.2

-121.4

-121.2

-121

-120.8

-120.6

-120.4 Figure 7.3. Map of Simulated Peak Ground Velocity (PGV).

TEMP EAST NORTH VERTICAL 28.3 55.5 16.3

..,1.

I

.. 1..

1l FIL A !Po 110MI-"L-00 yo.

_s 6o10sec recorded ground velocity simulated (actual hypocenter) simulated (reverse hypocenter)

Figure 7.4. Comparison of observed (top) and simulated ground velocities in units of (cm/sec) at station TEMP (Templeton). The 2003 San Simeon mainshock simulated ground motions were performed using the actual (center) and reversed (bottom) hypocenter locations.

A-48

5 PG&E Letter DCL-05-002 Rupture directivity effects depend strongly on where the fault rupture begins, i.e. on the location of the hypocenter. The San Simeon earthquake began to rupture at its northwest end, and the rupture propagated to the southeast, causing rupture directivity effects toward the southeast. To test the sensitivity of the ground motion maps to the location of the hypocenter, we placed the hypocenter at the southeast end of the fault rupture instead of at its actual location at the opposite (northwest) end of the fault. The resulting maps of peak acceleration and peak velocity using the reversed hypocenter are shown in Figures 7.5 and 7.6.

Rupture directivity effects result from the cumulative effect of rupture propagation toward a site, and so they depend on the distribution of slip on the fault as well as on the geometry of the spreading rupture front. Consequently, the ground motion maps for the actual and reversed hypocenters have differences that are greater than those of mirror images.

The peak acceleration map has a more concentrated region of high values at the northwest end of the second fault segment, and lower values in the southeast part of that fault segment, compared with the map that uses the actual hypocenter. As we expect, the peak velocity maps have much larger differences than the peak acceleration maps because rupture directivity effects are stronger for long period ground motions that for short period ground motions. The peak velocity map for the reversed hypocenter has its largest values near the northwestern end of the second fault segment, and large values extend further to the northwest off the end of that fault segment. The peak velocities in the southeast part of the second fault segment are much lower than for the actual hypocenter.

In Figure 7.4, we compare the ground motion velocity time histories recorded at Templeton with the simulated ground motions for the actual hypocenter and for the reversed hypocenter. We have already noted that our ground motion simulation for the actual hypocenter is unable to explain the large amplitudes of the ground motions recorded at Templeton, but the waveforms and durations of the recorded and simulated ground motions are quite similar. This is in contrast with the simulation using the reversed hypocenter, whose waveform lacks the large long period motions of both the recorded ground motions and those simulated using the actual hypocenter, and has longer duration. The amplitudes of the ground motions simulated using the reversed hypocenter are about half those simulated using the actual hypocenter. This figure clearly demonstrates the effect of different rupture directions on ground motion amplitudes, waveforms and durations, caused by differences between fonvard and backward rupture directivity effects.

From this comparison, we conclude that relatively strong ground motions occurred toward the southeast of the hypocenter of the San Simeon earthquake at sites located above the fault plane due to rupture directivity effects. The ground motions to the southeast would have been weaker, and those to the northwest would have been stronger, if the earthquake had begun to rupture at the southeast end of the fault and propagated to the northwest.

A49 PG&E Letter DCL-05-002 36 35.8 35.6 35.4 35.2

-121.4

-121.2

-121

-120.8

-120.6

-120.4 Figure 7.5. Map of simulated peak ground acceleration (PGA) using rupture velocity of 3.4 km/s and reverse hypocenter location compared to Figure 7.2.

36 -

0 60

\\

PGV (cm/s)

\\\\

35.8 35.6 35.4 10 km

{

35.2I

-121.4

-121.2

-121

-120.8

-120.6

-120.4 Figure 7.6. Map of simulated peak ground velocity (PGV) using a rupture velocity of 3.4 km/s and reverse hypocenter location compared to Figure 7.3.

A-SO C59 PG&E Letter DCL-05-002

8.

SUMMARY

The objectives of the work described in this report are to develop a rupture model of the 2003 San Simeon earthquake and to use it to assess the effects of rupture directivity on the ground motions from the earthquake. The first step was to derive reliable focal mechanisms and source depths of the larger aftershocks from regional long period seismograms. These mechanisms and focal depths were then used to calibrate the regional seismic velocity structure, and Green's functions for the velocity structure were calculated for use in deriving a finite fault rupture model of the earthquake. Teleseismic data were used to derive a multiple point source model of the earthquake rupture process. The teleseismic data were then combined with regional broadband, strong motion and geodetic data and used to develop a finite fault rupture model of the earthquake. This rupture model was used in a broadband ground motion simulation procedure to generate ground motion maps, which were used to assess rupture directivity effects in the recorded strong ground motions of the earthquake.

We acquired PG&E's waveform data for the mainshock and larger aftershocks recorded on its microearthquake and strong motion stations, especially waveforms that show both P and S waveforms.

We also acquired regional broadband seismic data, teleseismic data, and geodetic data. We used very long period regional seismic waves of the larger aftershocks to constrain their seismic moment tensors and focal depths. These focal depths and focal mechanisms are insensitive to the regional crustal structure, and so were useful for constraining the regional crustal structure using regional broadband seismograms.

We proceeded to use the source parameters of the large aftershocks to constrain the regional seismic velocity structure using waveform modeling. We obtained separate 1D seismic velocity models for each of six broadband recording stations. The velocity models that we obtained (Figure 3.6) are very similar to the Central California model used by McLaren and Savage (2001). These models were used in the inversion of the source rupture process of the mainshock.

We used teleseismic recordings of the earthquake to develop a representation of the source process in the form of multiple point sources (Figure 4.5). The teleseismic rupture model inversion was performed using the method of Kikuchi and Kanamori (1991). This method iteratively identifies subevents which are modeled as point sources but which represent extended regions on the fault having large energy release. An earthquake having a simple source on a single fault plane has a single subevent. We found that the San Simeon earthquake could be approximated by two subevents. These subevents are located approximately along strike of each other but have different strike angles, 298 and 320 degrees, indicating that a change in strike occurred as the rupture propagated toward the southeast. This change in strike is consistent with a change in strike of the mapped faults and of the aftershock zone. The first subevent is located near the hypocenter at a depth of 12 km, and the second subevent is located 12 km southeast of the hypocenter at a depth of 10 km. The seismic moment of the second subevent is 50% larger than that of the first subevent.

Using Green's functions calculated from these calibrated regional 1D seismic velocity models, we proceeded to develop a finite fault rupture model of the earthquake (Figure 5.2), using strong motion recordings, regional and teleseismic body waves and surface waves, and regional geodetic data. This rupture model describes the location of the hypocenter, the orientation of the fault rupture plane, and the distribution of slip on the fault. The rupture geometry model consists of two fault segments having different fault strikes. This change in the strike is indicated by three independent lines of evidence. First, the multiple point source model derived from teleseismic data indicated the presence of two subevents having different strike angles, as described above.

Second, there is a change in the strike of the aftershock zone. Third, there is a change in the strike of the mapped surface faults in this region.

Rupture on the first fault segment, which contains the hypocenter, was mostly downdip of the hypocenter and occurred in the depth range of 12 to 15 km. This is consistent with the first point source derived from the teleseismic data alone. Rupture of the second fault segment was largest on an asperity located about 12 km east-southeast of the hypocenter, at a depth of about 9 km. This is consistent with the second point source derived from the teleseismic data alone. Unlike the first fault segment, the A-S1 PG&E Letter DCL-05-002 second fault segment had a fairly uniform amount of slip.

The rupture model has a rupture velocity of 3.5 km/sec, which is roughly equal to the shear wave velocity of the upper crust.

The high rupture velocity, and the location of the hypocenter at the northwestern end of the fault rupture, contributed to relatively strong rupture directivity effects to the southeast, at sites located above the fault, as seen in the ground motion maps (Figures 7.2 and 7.3). The ground motions to the southeast would have been weaker, and those to the northwest would have been stronger, if the earthquake had begun to rupture at the southeast end of the fault and propagated to the northwest (Figures 7.5 and 7.6).

9. CONCLUSIONS

1. The centroid depth of the first subevent of the earthquake is 12 km +/- 2 kin, estimated from teleseismic body waves and surface waves.

2. The earthquake involved rupture on two fault segments having different fault strikes. The first large rupture event was located slightly downdip of the hypocenter. This initial event was followed 3.2 seconds later by a second event with a depth of 10 km located about 12 km to the southeast having a strike that is rotated 20 degrees clockwise with respect to the initial event.

3. This change in the strike is indicated by three independent lines of evidence. First, the multiple point source model derived from teleseismic data indicated the presence of two subevents having different strike angles. Second, there is a change in the strike of the aftershock zone. Third, there is a change in the strike of the mapped surface faults in this region.

4. The earthquake rupture propagated unilaterally to the southeast at a rupture velocity of about 3.5 km/sec, which is comparable to the shear wave velocity in the upper crust.

5. The high rupture velocity, and the location of the hypocenter at the northwestern end of the fault rupture, contributed to relatively strong rupture directivity effects to the southeast, at sites located above the fault. However, our earthquake source and wave propagation models are unable to explain the large ground motions recorded at Templeton. This may be due to the limited frequency resolution of our calculated Green's functions, or to the possibility that ground motion generation and propagation effects other than those included in our ground motion simulation procedure are responsible for the large ground motions recorded at Templeton.

Such unmodeled effects may include departures from 1D seismic velocity structure in the crust, and the presence of basin effects or unusual site amplification effects at Templeton.

A-52 PG&E Letter DCL-05-002

10. REFERENCES Ammon, C. J., G. E. Randall and G. Zandt (1990). On the non-uniqueness of receiver function inversions, J. Geophys. Res., 95, 15303-15318.

Bouchon, M. (1976). Teleseismic body wave radiation from a seismic source in a layered medium, Geophys. J. R. Astr. Soc., 47, 515-530.

Berteussen, K. (1977). Moho depths determination based on spectral-ratio analysis of NOSAR long-period P waves, Phys. Earth Planet. Inter., 15, 13-27.

CISN (2004). Performance of the CISN during the 2003 San Simeon earthquake (http://www.cisn.org/special/evt.03.12.22).

Crotwell, P. H., T. J. Owens, and J. Ritsema (1999). The TauP Toolkit: Flexible seismic travel-time and ray-path utilities, Seismol. Res. Lett., 70, 154-160.

Dreger, D. S., and D. V. Helmberger (1993). Determination of source parameters at regional distances with three-component sparse network data, J. Geophys. Res., 98, 8107-8125.

Graves R. G., and A. Pitarka (2004). Hybrid method for simulating broadband ground motions, World Conference Earthquake Engineering Proceedings.

Hardebeck, J. L., J. Boatwright, D. Dreger, R. Goel, V. Graizer, K. Hudnut, C. Ji, L. Jones, J. Langbein, J.

Lin, E. Roeloffs, R. Simpson, K. Stark, R. Stein, and J. C. Tinsley (2004). Preliminary report on the 22 December 2003, M 6.5 San Simeon, California Earthquake, Seismol. Res. Lett., 75, 155-172.

Hartzell, S. H., and T. H. Heaton (1983). Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California earthquake, Bull. Seismol. Soc.

Am., 73 1553-1583.

Haskell, N. A. (1960). Crustal relfection of plane SH-waves, J. Geophys. Res., 65, 4147-4150.

Haskell, N. A. (1962). Crustal reflection of plane P and SV-waves, J. Geophys. Res., 67,4751-4767.

Hauksson, E., D. Oppenhimer, and T. M. Brocher (2004). Imaging the source region of the 2003 San Simeon earthquake within the weak Franciscan subduction complex, central California, Geophys. Res.

Lett., 31, L20607, doi: 10.1029/2004GL021049.

Howie, J. M., K. C. Miller, and W. U. Savage (1993). Integrated crustal structure across the south central California Margin: Santa Lucia Escarpment to the San Andreas fault, J. Geophys. Res., 98(B5), 8173-8196.

Ichinose, G. A., H. K. Thio, and P. G. Somerville (2003). Rupture process of the 1944 Tonankai earthquake (Ms 8.1) from the inversion of teleseismic and regional seismograms, J. Geophys. Res., 108, B 10, 2497, doi: 10.1 029/2003JB002393.

Ji, C., K. M. Karson, Y. Tan, K. W. Hudnut, and K. Choi (2004). Slip history of the 2003 San Simeon earthquake constrained by combining 1-Hz GPS, strong motion, and teleseismic data, Geophys. Res.

Lett., 31, L17608, doi:10.1029/2004GL020448.

Kikuchi, M., and H. Kanamori (1991). Inversion of complex body waves-III, Bull. Seismol. Soc. Am.,

81, 2335-2350.

Langston, C. A., and D. V. Helmberger (1975). A procedure for modeling shallow dislocation sources, Geophys. J. R. Astr. Soc., 42, 117-130.

Lawson, C. L., and R. J. Hanson (1974). Solving Least Squares Problems, Englewood Cliffs, NJ:

Prentice-Hall.

Levander, A. R., and R. L. Kovach (1990). Shear velocity of the northern California lithosphere, J.

Geophys. Res., 95(B12), 19773-19784.

McLaren, M. K., and W. U. Savage (2001). Seismicity of South-Central Costal California: October 1987 through 1997, Bull. Seismol. Soc. Am., 91(6), 1629-1658.

A-53 PG&E Letter DCL-05-002 Nolet, G., J. van Trier, and R. Huisman (1986). A formalism for nonlinear inversion of seismic surface waves, Geophys. Res. Lett., 13(1), 26-29.

Oppenheimer, D., F. Klein, J. Eaton, and F. Lester (1993). The northern California Seismic Network Bulletin, January-December 1992, U.S. Geol. Surv., Open File Rept.93-578.

Paige, C. C., and M. A. Saunders (1982a). LSQR: an algorithm for sparse linear equations and sparse least squares: Assn. Comp. Mach Trans. Mathematical Software, 8,43-71.

Paige, C. C., and M. A. Saunders (1982b). Algorithm 583, LSQR: Sparse linear equations and least squares problems, Assn. Comp. Mach. Trans. Mathematical Software, 8, 195-209.

Pasyanos, M. E., D. S. Dreger, and B. Romanowicz, Toward real-time estimation of regional moment tensors, Bull. Seismol. Soc. Am., 86, 1225-1269, 1996.

Randall, G. E. (1994). Efficient calculation of complete differential seismograms for laterally homogeneous earth models, Geophys. J. Int., 118, 245-254.

Ritsema, J., and T. Lay (1995). Long-period regional wave moment tensor inversion for earthquakes in the western United States, J. Geophys. Res., 100, 9853-9864.

Scholz, C. H. (1990). The Mechanics of Earthquakes and Faulting, Cambridge University Press.

Shaw, P. R. (1998) Waveform inversion of time-corrected refraction data, Geophys. J. Int., 93, 75-90.

Somerville, P. G., K. Irikura, R. Graves, S. Sawada, D. Wald, N. Abrahamson, Y. Iwasaki, T. Kagawa, N. Smith, and A. Kowada (1999). Characterizing crustal earthquake slip models for the prediction of strong ground motion, Seismol Res. Lett., 70, 59-80.

Stewart, S. W. (1968). Preliminary comparison of seismic traveltimes and inferred crustal structure adjacent to the San Andreas fault in the Diablo and Gabilan Ranges of central California, in Proc. Conf.

Geological Problems of San Andreas fault system, W. R. Dickinson and A. Grantz, Editors, Stanford University Pub. Geol. Sci. 11, 218-230.

Thio, H. K., R. W. Graves, P. G. Somerville, T. Sato, and T. Ishii (2004a). A multiple time window rupture model for the 1999 Chi-Chi earthquake from a combined inversion of teleseismic, surface wave, strong motion, and GPS data, J. Geophys. Res., 109, B08309, doi:10.1029/2002JB002381.

Thio, H. K., R. W. Graves, P. G. Somerville, T. Sato, and T. Ishii (2004b). The source process of the 1999 Kocaeli (Turkey) Earthquake, J. Geophys. Res., (submitted to Geophys. J. Int.).

Walter, A. W., and W. Mooney (1982). Crustal structure of the Diablo and Gabilan Ranges, Central California: A reinterpretation of existing data, Bull. Seismol. Soc. Am., 72(5), 1567-1590.

Wills, C.J., M. Petersen, W.A. Bryant, M. Reichle, G.J. Saucedo, S. Tan, G. Taylor and J. Treiman (2000). A site-conditions map for California based on geology and shear-wave velocity. Bull. Seismol.

Soc. Am., 90, S187-S208.

Xu, Y., and D. A. Wiens (1997). Upper-mantle structure of the southwest Pacific from regional waveform inversion, J. Geophys. Res., 102, 27439-2745 1.

Zeng, Y., and J. G. Anderson (1995). A method for direct computation of the differential seismogram with respect to the velocity change in a layered elastic media, Bull. Seismol. Soc. Am., 85, 300-307.

A-54 PG&E Letter DCL-05-002 Engineering Update to Supplemental Report Dated March 29, 2004 San Simeon Earthquake of December 22, 2003 This enclosure provides a final update to the engineering related information pertaining to the San Simeon earthquake of December 22, 2003, reported by PG&E on March 29, 2004 (reference. "Supplement to Special Report 03-04: San Simeon Earthquake of December 22, 2003," PG&E Letter DCL-04-031). The following areas are addressed:

• The cumulative absolute velocity (CAV) value for free field earthquake ground motion records recorded at Diablo Canyon Power Plant (DCPP).

The excursion above the Design Earthquake response spectrum associated with the vertical component of the San Simeon earthquake record at the top of the Unit 1 containment structure.

1 of 6

-I PG&E Letter DCL-05-002 Engineering Update to Supplemental Report Dated March 29, 2004 San Simeon Earthquake of December 22, 2003 Updated Cumulative Absolute Acceleration Value Introduction The cumulative absolute velocity (CAV) value from the free field earthquake ground motion records serves as a check for exceedance of the Operating Basis Earthquake in accordance with Reference 1. Information on the CAV value associated with the DCPP ground motion records for the San Simeon Earthquake of December 22, 2003, was initially submitted to the NRC on March 29, 2004, (Reference 2).

Discussion The CAV value provided in Reference 2 was computed using the uncorrected time histories for the three components of free field ground motion, as measured by the seismic instrument at the Unit 1 containment basemat. Subsequently, baseline correction was applied to the time histories and PG&E recomputed the CAV value and found it to be slightly greater than that previously reported. As a result of the correction, the CAV increased from the 0.07 g-sec value reported in Reference 2 to 0.082 g-sec.

The baseline corrected CAV value is still significantly less than the threshold value of 0.16 g-sec., below which earthquake ground motions are not considered to have the potential for causing damage to nuclear power plants (Reverence 1). This updated information was previously reported to the NRC at a public meeting held at the NRC headquarters in Rockville, Maryland, on May 27, 2004.

References

1. "A Criterion for Determining Exceedance of the Operating Basis Earthquake,"

Electric Power Research Institute Report No. NP-5930, July 1988.

2. "Supplement to Special Report 03-04: San Simeon Earthquake of December 22, 2003," PG&E Letter DCL-04-031, dated March 29, 2004.

2 of 6

N PG&E Letter DCL-05-002 Engineering Update to Supplemental Report Dated March 29, 2004 San Simeon Earthquake of December 22, 2003 Engineering Study - Excursion above the Design Earthquake Response Spectrum Associated with the Vertical Component of the San Simeon Earthquake Record at the Top of the Unit 1 Containment Structure Introduction The engineering evaluation of the recorded response of DCPP during the San Simeon earthquake of December 22, 2003, demonstrated that dynamic amplification of the vertical component of the input seismic ground motion was present for the DCPP containment structure. The response spectrum associated with the vertical component of the San Simeon earthquake record at the top of the Unit 1 containment structure showed an excursion above the Design Earthquake (DE) response spectra for this location. This information was reported to the NRC in PG&E Letter DCL-04-031 dated March 29, 2004 (Reference 2). As a follow-up to PG&E Letter DCL-04-031, PG&E performed a detailed engineering study of the impact of vertical dynamic amplification on the seismic design of DCPP and reported the results of the study to the NRC at a public meeting held at the NRC office in Rockville, Maryland on May 27, 2004.

This report provides a brief summary of the information provided by PG&E during the public meeting with the NRC and provides closure to the issue of the earthquake response excursion reported in PG&E Letter DCL-04-031.

Background

The seismic design basis for DCPP associated with the DE and Double Design Earthquake (DDE) does not consider the effects of dynamic amplification of the vertical component of the input ground motion based on the industry practice at the time, which considered large reinforced concrete structures to be rigid in the vertical direction.

However, the top of the Unit I containment structure acted as a diaphragm to amplify the vertical motion during the San Simeon earthquake. This effect was recognized during the Hosgri Earthquake (HE) reevaluation of DCPP and was explicitly considered for the HE evaluation of structures, systems, and components (SSCs). Subsequently, this effect was also considered in the Long Term Seismic Program (LTSP) evaluation of DCPP.

The basis for the acceptability of the DCPP seismic design is that the HE reevaluation of SSCs, which resulted in significant strengthening of DCPP, is governing such that the nonconservative exclusion of vertical amplification for the DE and DDE does not adversely impact the seismic qualification of DCPP. In order to support this conclusion with quantitative analyses, PG&E undertook an engineering study to validate the adequacy of the DE/DDE seismic design of safety related SSCs, considering amplified vertical ground motion.

3 of 6 PG&E Letter DCL-05-002 Engineering Update to Supplemental Report Dated March 29, 2004 San Simeon Earthquake of December 22, 2003 Scone Selection for Engineering Study The engineering study addressed the seismic qualification of the following areas of the plant:

Major Structures Containment Structure

• Auxiliary Building Systems Piping Electrical Raceways HVAC Ducts Components

• Equipment Qualified by Testing

• Equipment Qualified by Analysis The scope of the study was based on a sample of SSCs and excluded many major components that are in plant locations where the amplification of vertical motion will not occur. These areas/SSCs include the following:

• SSCs in Vertically Rigid Areas: Reactor, Reactor Coolant Pumps, Steam Generators, Polar Crane Wall, Containment Cylinder, Engineered Safety Feature Pump Rooms

• SSCs on Grade: Emergency Diesel Generators, Auxiliary Saltwater Pumps and Buried Piping, Component Cooling Water Heat Exchangers, Emergency Water Tanks (Condensate, Refueling Water, Firewater)

Analysis and Evaluation of Structures The analysis and evaluation of structures included the development of finite element models to capture the vertical dynamic response, as well as generation of amplified vertical response spectra for the DE and DDE, for use in the seismic evaluation of safety-related systems and components. The resulting stresses of structural components were evaluated against the design basis allowable stresses for the applicable load combinations when subjected to amplified vertical input motion. The damping ratios applied in the structural analyses were consistent with those of Regulatory Guide 1.61 (Reference 1).

4 of 6 PG&E Letter DCL-05-002 Engineering Update to Supplemental Report Dated March 29, 2004 San Simeon Earthquake of December 22, 2003 The study included the containment structure exterior concrete and two floor slabs at elevation 140' in the auxiliary building. The results of the study demonstrated that these structures satisfied their design requirements when vertical amplification is considered.

Evaluation of Systems The study considered the major categories of systems at DCPP. Key conclusions reached are summarized as follows:

• Five selected piping systems in the Unit 1 portion of the auxiliary building and the Unit 1 containment structure were reevaluated using the DE and DDE vertically amplified response spectra as input. These evaluations demonstrated that pipe stress levels and support loads satisfied the applicable design criteria when vertical amplification is considered.

• The plant vent was reevaluated using the DE and DDE vertically amplified response spectra as input. However, since the current design basis considers the envelope of the DE, DDE, and HE, and the HE is still the governing case, the consideration of vertical amplification for the DE and DDE did not impact the qualification of the plant vent.

• Electrical raceway supports and HVAC duct supports were reevaluated using the DE and DDE vertically amplified response spectra as input. Since the current design basis does not require qualification for the DE and the seismic input is based on the envelope of the DDE and HE, the consideration of vertical amplification for the DE and DDE did not impact the qualification of these commodities.

Evaluation of Components The seismic qualification of equipment for DCPP is performed through shake table testing, analysis, or a combination of both methods. The approach used in this study was dependent on the qualification method.

Shake table testing is based on the provisions of IEEE 344. Qualification by shake table testing to the Safe Shutdown Earthquake (SSE) level (corresponding to the HE) and to the Operating Basis Earthquake (OBE) level (at least 50 percent of the HE) accounts for amplification of the vertical input motion, since amplification is considered in the development of the vertical HE response spectra.

• Qualification by analysis compared the existing evaluations for the HE, which included vertical amplification, with the vertical DE response spectra generated for this study, considering the ratio of allowable stresses for the HE vs. DE to determine the governing load combination. For a sample of equipment calculations, and for 5 of 6 PG&E Letter DCL-05-002 Engineering Update to Supplemental Report Dated March 29, 2004 San Simeon Earthquake of December 22, 2003 the response spectra generated for this study, it was found that the qualification for the HE governs over the DE, even when vertical amplification is considered.

Conclusions When amplification of the vertical component of the DE and DDE is considered, the SSCs covered in the study meet required acceptance criteria. The current design basis, which considers the vertical amplification of the HE assures that there is adequate margin in the seismic design of DCPP, even though vertical amplification is not considered for the smaller DE and DDE.

References

1. "Damping Values for Seismic Design of Nuclear Power Plants," Regulatory Guide 1.61, October 1973.

2. "Supplement to Special Report 03-04: San Simeon Earthquake of December 22, 2003," PG&E Letter DCL-04-031, dated March 29, 2004.

6 of 6