Bulletin of the Seismological Society of America, Vol. 98, No. 2, pp. 778792, April 2008, doi: 10.1785/0120070119
Explosion-Source Energy Partitioning and Lg-Wave Excitation: Contributions of Free-Surface Scattering
by Yaofeng He, Xiao-Bi Xie, and Thorne Lay

Abstract

A 2D boundary-element numerical simulation approach and a local slowness analysis method for an embedded array are used to quantify effects of topographic scattering on near-source energy partitioning for simple underground explosion sources. Various parameters, including free-surface models with different root mean square (rms) random topographic fluctuations and correlation lengths, source depths from 0.25 to 3 km, and Q values from 50 to infinity are included in the numerical simulations, with energy responses of different phases being determined as functions of frequency. The results reveal that for a crustal model with a relatively high surface velocity, near-source free-surface scattering provides an important coupling mechanism that can impart additional explosion energy to the Lg wave. At relatively low frequencies, and for a moderately rugged free-surface, the Rg-to-Lg transfer is quite efficient, while at higher frequencies or for a very rugged free surface, the body wave to Lg transfer may dominate the process. The Rg excitation functions, source depth, and topographic correlation length all contribute frequency dependence to the Lg excitation function. The presence of a low Q value within the uppermost crust severely attenuates the high-frequency energy transferred to the Lg wave.

Introduction

Regional seismic phases (e.g., Pn, Sn, Pg, and Lg) play an important role in global monitoring of low-yield underground nuclear tests. Numerous empirical observations have shown that regional phases hold the keys to small-event magnitude and yield estimation and to discrimination between small explosions and earthquakes (e.g., Nuttli, 1986; Taylor et al., 1989; Kim et al., 1993; Walter et al., 1995; Fisk et al., 1996; Taylor, 1996; Kim et al., 1997; Hartse et al., 1997; Taylor and Hartse, 1997; Fan and Lay, 1998a,b,c; Patton, 2001; Bottone et al., 2002; Xie, 2002). However, building a sound physical basis for the application of the empirical relationships requires detailed knowledge of the regional phases (e.g., relative intensity and frequency dependence of different phases excited by earthquakes and explosions) and relationships linking the wave characteristics with a range of source and structure parameters. The complex excitation and energy-partitioning mechanisms yielding regional phases are difficult to empirically separate by data analysis. Thus, numerical modeling approaches are valuable for investigating excitation and propagation of regional seismic phases (e.g., McLaughlin and Jih, 1988; Xie and Lay, 1994; Bradley and Jones, 1998, 1999; Wu et al., 2000a,b; Bonner, et al., 2003; Stevens et al., 2003; Myers et al., 2005; Stevens et al., 2005; Xie, Ge, and Lay, 2005; Xie, Lay, and Wu, 2005). 778 Although there are continuing controversies over the relative importance of various energy-partitioning mechanisms affecting regional phases, most investigators agree that appreciable S-wave excitation for explosion sources occurs in the near-source region, reducing the performance of event discrimination approaches. Several possible near-source S-wave energy excitation mechanisms have been proposed to explain the generation of explosion Lg waves. Among these, near-source coupling between P, S, and Rg waves due to scattering at a rugged free surface may play an important role in Lg-wave excitation. This has been investigated by different authors from both observational and theoretical perspectives. From data analysis, Gupta et al. (1992, 2005) suggested that near-source scattering of explosion-generated Rg into S makes a significant contribution to low-frequency Lg signals. With Rg being strongly excited for very shallow explosions relative to deeper earthquakes, efficient Rg-to-Lg scattering may cause the low-frequency P=Lg ratio to fail to discriminate source type. Patton and Taylor (1995) analyzed the Lg spectral ratios from the Nevada Test Site (NTS) explosions and suggested that the Lg wave is generated by near-source scattering of Rg waves into body waves, which become trapped in the crust. Myers et al. (1999) investigated the 1997 depth of burial experiment at the former Soviet test site at Balapan, Kazakhstan. By comparing the

Explosion-Source Energy Partitioning and Lg-Wave Excitation regional and local recordings from explosions at different depths, they suggested that the data support Rg scattering as a major source of explosion S waves. Patton and Taylor (1995) and Gupta et al. (1997) introduced theoretical models of the explosion spall source to explain the observed similarity between Rg and Lg spectra. McLaughlin and Jih (1988) used finite-difference simulation to investigate topography influences on teleseismic P waveforms, and indicated possible Rg-to-P scattering due to the near-source topography. More recently, Bonner et al. (2003) and Wu et al. (2003) provided strong evidence in favor of the Rg-to-S scattering mechanism for the generation of the lowfrequency S and Lg for explosions. Stevens et al. (2005) modeled the Rg-to-Lg scattering by assuming that Lg is generated by a distribution of surface scatters. Xie, Ge, and Lay (2005) and Xie, Lay, and Wu (2005) investigated the contribution of shallow volumetric scattering to explosion-source energy partitioning and calculated the frequency-dependent Lg excitation functions. They found that the high-frequency Lg energy is mainly from P-pS-to-Lg and P-to-Lg scattering, while the low-frequency energy is mainly from Rg-to-Lg scattering. Myers et al. (2005), using numerical simulation, investigated the effect of surface topography on the P-to-S conversion. They concluded that near-source topography and geologic complexity in the upper crust strongly contribute to the generation of S waves. Given the existence of a rugged free surface, the actual formation and coupling between waves in the near-source region is expected to be rather complex, so it is important to study this phenomenon with the source excitation included in the model, not simply as a remote propagation effect. Scattering at the rough free surface can change the propagation direction of pP and pS waves, causing more of their energy to become trapped in the crustal wave guide to contribute to the Lg wave than would occur for a flat surface. A rugged free surface and/or shallow heterogeneity also provides coupling between surface waves and body waves. Both body-wave to surface-wave and surface-wave to body-wave scattering can occur. Multiple scattering, variable source depth, and attenuation in the shallow layers are all factors that may affect the frequency-dependent regional wave energy partitioning and these effects need to be quantified. In this article, we use a 2D P-SV boundary-element simulation (Ge et al., 2005) and an embedded array-slowness analysis method (Xie, Ge, and Lay, 2005) to investigate the effect of topographic scattering on explosion-source energy partitioning. Crustal models with random rough free surfaces are used in numerical simulations and frequency-dependent response functions for different source/model parameters (e. g., source depth, free-surface roughness, and intrinsic attenuation) are considered relative to the flat free-surface case. An isotropic explosion source is used throughout this article, and nonlinear effects, such as spall or cracking near the source, are not considered. The results show that surface scattering does cause coupling between the body and surface

Y. He, X.-B. Xie, and T. Lay

Distance km

Random 0 topography
Explosion source Crust model (Vp,Vs, )

Vertical receiver array

Figure 1.
array. Configuration of the source, model, and receiver
complex. Factors such as the source depth, local layered structure, attenuation, random volumetric velocity perturbations, and free-surface fluctuations all affect the partitioning. These effects often contribute to the partitioning in a coupled way, and the entire process is not necessarily linear or simply separable. In this study, we focus on the contribution from surface scattering. We use accurate numerical modeling to simulate the complex partitioning process and the slowness analysis to calculate the discrete response functions. The previously mentioned symbolic equations provide us with a basic formalism for understanding the process. The embedded-array method cannot characterize all partitioning coefficients, but is particularly well suited to characterizing relative changes in the trapped energy distribution for waves that will travel to large distances in the crustal wave guide.

Depth km

The Energy-Partitioning Formalism For convenience, we symbolically write the near-source energy-partitioning process for an explosion source as EK p; f S f RK p; f ;
Numerical Examples of Free-Surface Scattering
Space-Domain Representation of Surface Scattering Illustrated in Figure 2 are boundary-element-generated snapshots for wave fields in a model with free-surface scattering. The parameters of the three-layer velocity model are listed in Table 1. A random fluctuation with a correlation length of 0.5 km and an rms fluctuation of 0.15 km is used for the free surface and extends between distances of 30 and 50 km. The source is located at distance 20 km and at a depth of 0.5 km. Figure 2a,b presents horizontal and vertical displacements. In addition to familiar major phases (e.g., P, pS, and Rg) expected for a shallow explosion in a flat-earth model, scattered body and surface waves from the rough free surface are present in the wave field. The surface-to-body and body-to-body wave scattering is distributed through the entire medium following the direct waves, and the bodyto-surface and surface-to-surface wave scattering is concentrated at very shallow depths following direct waves as they graze the surface. To isolate the scattered phases, we subtract the wave field generated for a flat surface from the wave field for the model having a rough surface, yielding the results presented in Figure 2c,d. Most of the scattered body waves are shear waves. Because of the coupling between different wavenumbers, the scattered body waves have a very broad range of propagation directions. The horizontal component includes mostly shear waves propagating with steep dip angles that will tend to contribute to teleseismic S waves. The shear waves on the vertical component mostly have shallow angles and will contribute to crustal guided waves such as regional Lg. Although these space-domain images are instructive for understanding some aspects of the surface scattering, to fully explore the wave-propagation characteristics in the complex near-source environment, we conduct slowness analysis for these wave fields. We note that any differences in wave excitation due to model variations relative to the reference flat-surface case are captured in the differenced wave fields and may not be distinguishable from scattering

where EK p; f is the near-source energy partitioned to the type K wave (K can be P, S, Lg, Rg, or other wave types), p is the slowness, f is the frequency, and S f is the spectrum of an isotropic explosion source. RK p; f is the energy response function of the near-source structure for exciting type K wave, and can be expressed as RK p; f RK p; f F X

X RJ p; f T JK p; f F

p; f T KJ p; f :
On the right-hand side of this equation, RK p; f is the reF sponse of a flat, homogeneous layered earth model, partitioning the source energy into different phases. The transfer function T JK provides the J-to-K coupling, which modifies the original partitioning by moving energy from one phase to another. The second term on the right-hand side denotes energy being imparted into the K wave through coupling, and the third term denotes energy lost from the K wave to other phases. The combined effect gives the total partitioning of the energy radiated from an isotropic source into the K-wave energy distributed in slowness and frequency domains. This energy will develop into different regional phases, which propagate to remote distances. Having a complete description of the slowness distribution allows us to accurately predict energy imparted to the distant regional phases based on the near-source energy budget. Investigating these response and transfer functions provides a way to estimate the underlying process of energy partitioning. Because of the diverse mechanisms involved, the actual near-source energy partitioning can be highly
Explosion-Source Energy Partitioning and Lg-Wave Excitation

Figure 2.

Wave field snapshot at t 10:0 sec for a model with random free-surface fluctuation, where (a) and (b) are horizontal and vertical components of the wave field, and (c) and (d) are horizontal and vertical components of the scattered wave field obtained by subtracting the flat model wave field from the random surface wave field. The source is at 20-km horizontal position and 0.5-km depth.
effects, but statistical averaging over multiple-model realizations does tend to isolate the scattering effects. Slowness-Domain Presentation of Surface Scattering To quantify the scattered wave field, we transform the signal to the slowness domain. We compute wave fields for the two-layer velocity model listed in Table 2. The random surface topographic fluctuation is located above the source and extends in both directions for 20 km. The random topography has an exponential power spectrum, and its correlation distance, a, is 0.5 km. The rms fluctuations used are 0.15 km (maximum peak-to-trough 0.625 km) and 0.3 km (maximum peak-to-trough 1.281 km). Either a shallow source (depth 0:5 km) or a deep source (depth 3:0 km) is used in the simulation. The synthetic seismograms computed across a vertical array (see Fig. 1) are processed using the slowness-analysis method (for details, see Xie, Ge, and Lay, 2005). Illustrated in Figure 3 are energy distributions in a mixed horizontal slowness and depth domain, where different rows are for different combinations of source Table 1

Three-Layer Velocity Model
Bottom of Layer (km) V P (km=sec) V S (km=sec) (g=cm3 )
and free-surface parameters, and the successive frames in each row are for different time windows. The time progression allows ready identification of major phases. Phases such as Rg, pS, S , and trapped energy from different mechanisms are labeled in the Figure. Frequency filters can be applied to the synthetic data before conducting the slowness analysis. When this is viewed collectively, we know the energy distribution in combined domains of space, slowness, time, and frequency. Hereafter, we will call this SPTF domain characterization, with the P standing for slowness. In Figure 3, the solid vertical lines indicate the upper mantle S-wave slowness. Wave energy to the left of these lines has incidence angles steeper than the critical angle on the Moho, and the energy will leak to the upper mantle through multiple reflections. For wave energy to the right of these lines, total reflection will keep the energy in the crustal wave guide, ultimately forming the Lg wave at long distances (e.g., Frankel, 1989; Xie and Lay, 1994; Vogfjrd, 1997; Xie, Ge, and Lay, 2005). Figure 3a is for a flat free surface and a relatively deep explosion source. As expected, this configuration generates neither noticeable trapped Table 2

Two-Layer Velocity Model

Infinity

5.6 6.5 8.0

3.2 3.8 4.5

2.7 2.9 3.3

45 Infinity

6.5 8.0

3.6 4.5

2.9 3.3

Figure 3. Slowness analysis results in the depth-slowness domain for discrete time intervals as the wave sweeps through the observing array for models with different source depths and free-surface parameters. The frequency band is between 1.5 and 4.5 Hz. All the panels are normalized in the same scale. In each small figure, the horizontal coordinate is horizontal slowness, and the vertical coordinate is depth. The thick vertical line indicates the upper mantle S-wave slowness which separates energy that leaks out of the wave guide (to the left) from energy trapped in the wave guide (to the right) that forms Lg. The PP , PS and PSM are crustal P-wave, S-wave, and upper-mantle S-wave slownesses, with the values of 0.154, 0.28, and 0:22 sec =km, respectively. energy nor clear Rg wave. Figure 3b is for a flat free surface and a shallow explosion source. We now see the Rg wave developed at shallow depth and trapped energy originating from the S wave (Vogfjrd, 1997; Xie, Ge, and Lay, 2005). Figure 3c is for a deeper explosion source and a free surface with 0.15-km rms fluctuation. Compared to 3a, the existence of a random free surface generates significant trapped energy from surface scattering. Although the source is located at a depth of 3.0 km, Rg energy can now be seen at shallow depth. This enhanced Rg wave comes from the freesurface scattering, which can be treated as shallow secondary sources. In Figure 3df, with shallower source or larger rms free-surface fluctuations, a substantial amount of trapped energy can be generated from interactions between the explosion source and the topographic fluctuations. The energy partitioned to specific regions in SPTF domain can be collected and used to estimate the energy input into different regional phases. Using different source/model parameters, the relationship between these parameters and the energypartitioning processes is investigated and the underlying mechanisms revealed.

Contributions of Free-Surface Scattering to Lg-Wave Excitation
We adopt the three-layer crust (Table 1) as our basic velocity model and add random free-surface fluctuations with
Explosion-Source Energy Partitioning and Lg-Wave Excitation Table 3
Source and Model Parameters Used in Numerical Simulations
Group Fixed Parameters Variables
0:5 km Source depth 0:5 km QP infinity QS infinity rms 0:15 km a 0:5 km QP infinity QS infinity rms 0:15 km Source depth 0:5 km QP infinity QS infinity a 0:5 km Source depth 0:5 km QP 100 QS 50
rms free surface fluctuation Source depth
0:0; 0:05; 0:10; 0:15; 0:20; 0:25; 0:30; 0:35; 0:40 km
0:25; 0:50; 0:75; 1:00; 1:50; 2:00; 2:50; 3:00 km 0:4; 0:6; 0:8; 1:0; 2:0; 4:0; 6:0; 8:0; 10:0 km 0:0; 0:05; 0:10; 0:15; 0:20; 0:25; 0:30; 0:35; 0:40 km

Correlation length a

rms free surface fluctuation
different statistical parameters to this basic model. The random topography has an exponential power spectrum. It is located above the source and extends in both directions for 20 km. The model geometry is similar to that shown in Figure 1, with a vertical array composed of receivers located between distances of 30 and 50 km and depths of 0 and 45 km. To isolate the effects of individual factors, we vary individual parameters while keeping other parameters unchanged. The varied parameters are listed in Table 3. Each model is described by a set of parameters including the rms free-surface fluctuation, the correlation length of the random power spectrum, the source depth, and the intrinsic attenuation (quality factor Q). The source
depth is measured locally by taking the vertical distance between the source and the rough surface. This avoids coupled variation of depth and topographic parameters. To characterize the results statistically, we generate 10 realizations for each model. Synthetic seismograms are calculated for each realization and processed separately. We then average the measurements from individual realizations and use their mean value as the final result for a particular case.
Effect of Free-Surface Roughness To investigate the effect of free-surface roughness on the energy partitioning, we use an infinite Q for both P and S

Figure 4. Slowness analysis results for models with different rms topographic fluctuations. The source depth is 0.5 km. In each time frame the horizontal coordinate is the horizontal slowness and vertical coordinate is the depth. The Rg energy is located near the surface, with slowness similar to the S wave.
784 waves, a source depth of 0.5 km, a correlation length of 0.5 km, and a varied rms free-surface fluctuation between 0.0 and 0.4 km (see group 1 in Table 3). The slownessanalysis results are illustrated in Figure 4 with different panels for models with different rms values. A frequencydomain filter between 1.5 and 4.0 Hz was applied to the synthetic seismograms before this slowness analysis. The Rg energy can be seen at depths less than 3 km with a slowness similar to the S wave. In Figure 4a, with rms 0, the Rg wave is directly generated entirely by the explosion source. It arrives at the receiver array at between 10 and 12 sec and is labeled as direct Rg wave. The presence of a rough free surface causes scattering of different waves and redistributes their energy. As shallow secondary sources, the scattering generates scattered Rg waves, which can be observed in all time windows in panels 4b,c. The shallow energy between 8 and 10 sec is from body wave to Rg scattering, as it arrives at the receiver array earlier than the direct Rg. We label this as scattered Rg wave and use it to investigate body- to surface-wave scattering. Because of scattering attenuation, the same surface fluc(a) direct Rg energy
tuations that excite the scattered Rg can also attenuate both direct and scattered Rg, generating scattered body waves. In panels 4d,e, due to strong scattering from a very rugged free surface, both direct- and scattered-Rg waves are very weak. Using panel 4a as a reference, the trapped energy can be quantified as a function of surface roughness. Applying SPTF processing to the slowness analysis, we can separate the energy and estimate the excitation of different phases (Xie, Ge, and Lay, 2005). A series of band-pass filters is used to obtain responses between 1.0 and 4.0 Hz. Figure 5ac illustrates the near-source response functions of direct Rg waves, scattered Rg waves, and Lg waves (summed trapped energy) as functions of frequency and rms surface fluctuations. The Rg-wave energy is obtained by picking energy peaks between depth 0 and 3 km and horizontal slowness 0:180:40 sec =km. The Lg-wave energy is picked from peaks between depths 3 and 45 km, and slowness 0:220:40 sec =km. The vertical coordinate is the square root of normalized energy E=E0 1=2, where E0 is a normalization factor that has the same source-time function and passes through the same frequency filter as that used in

(c) Lg energy

(b) scattered Rg energy

0.05 0.04

0.006 0.005 0.004
0.16 0.14 0.12 0.08 0.06 0.04 0.02 0.00 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.10

(E/E0) 1/2

0.003 0.002 1.0 1.5 0.001 2.0 0.000 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 2.5
0.02 1.0 1.5 0.01 2.0 0.00 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 2.5

1.0 1.5

y H 3.5 z 4.0

nc 3.0 y H 3.5 z 4.0

(d) scattered loss of Rg energy
(e) net scattered Rg energy
(f) scattered gain of Lg energy

-0.05 -0.04 -0.03

0.16 0.14 0.12 0.08 0.06 0.04 0.02 0.00 3.0 3.5 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.10
-0.01 1.0 1.5 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.00 2.0 0.01 2.5 0.02 3.0 3.5
0.002 1.0 1.5 0.001 2.0 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 2.5 -0.001 3.0 3.5
Responses as functions of frequency and rms free-surface fluctuations. The top row shows near-source responses of direct Rg, scattered Rg, and the Lg waves, with (a) RRgdirect , (b) RRgscatt , and (c) RLg. The bottom row shows the contributions of surface scattering to these responses, with (d) RRgdirect RRgdirect , (e) RRgscatt RRgscatt , and (f) RLg RLg. Note that a negative vertical coordinate is used in (d), where the F F F prisms with solid black tops are below the zero plane while prisms with open tops are above the zero plane.

Figure 5.

Explosion-Source Energy Partitioning and Lg-Wave Excitation the simulation. The source spectrum has been removed from the normalized energy, and the results are equivalent to the response functions in equations (1) and (2). Because E f2 A2 , where A is the displacement amplitude, the square root of energy can be compared to the wave amplitude after being scaled with f 1. For the response function RRgdirect , the energy is mainly located at low frequencies and drops quickly with increasing rms values. Note that the histogram labeled with rms 0 indicates the response for the flat-earth model. The net contribution of the surface scattering to direct Rg waves can be obtained by subtracting the flat-earth response from the total response (Fig. 5d). The generally negative values indicate the energy loss from direct Rg due to the scattering attenuation, which is proportional to the rms values. Note that a negative vertical scale is used in Figure 5d, and the prisms with solid black tops are plotted below the plane zero, while prisms with open tops are plotted above the plane. At very low rms values, the Rg wave gains some energy, which appears to be caused by body-tosurface-wave scattering providing slightly more Rg energy than the energy loss from Rg due to scattering. The time windows of direct and scattered Rg waves partially overlap. The application of a narrow-band filter causes broadening of the Rg impulse and results in some leakage of direct Rg-wave energy into the scattered Rg time window (see first row in Fig. 5b). By subtracting this energy, we ob-

Explosion-Source Energy Partitioning and Lg-Wave Excitation listed in group 3 of Table 3. The slowness analysis results are shown in Figure 8, with each panel for a model with different correlation length. The response functions of Rg, scattered Rg, and Lg waves are shown in Figure 9ac. The horizontal coordinate is frequency, and the different rows are for different correlation lengths. From these response functions, we see that the last row (with a correlation length of 10 km) is almost the same as that for a flat-earth model. This indicates that a very smooth, long-wavelength, free-surface fluctuation has almost no effect on these waves. Figure 9de shows net energy loss or gain from scattering. For models with correlation length shorter than 4 km, the surface scattering apparently contributes to the generation of the trapped energy in the wave guide (Fig. 9f). For very long correlation lengths, the random free surface behaves more like a flat free surface. In Figure 9f, the net scattered contribution to the Lg wave, the response function falls with increase of the correlation length at all frequencies. For frequencies used in the simulation (14 Hz) and S-wave velocity in the top layer (3:2 km=sec), the wavelengths are between 0.8 and 3.2 km. We calculate the response as a function of normalized scale length, ka, where k 2=, is the wave length, and a is the correlation length, and we present the behavior in Figure 10. The maximum scattering happens around ka 1 and decreases for larger ka. Extension of the calculation
787 to smaller ka is limited by the grid size used in the boundaryelement calculation and the dimension of the receiver array. The Effect of Intrinsic Attenuation Scattering from topographic fluctuations occurs in the uppermost crust, which is usually a low Q layer. In addition, the scattering increases the propagation distances, especially for high-frequency waves. Attenuation will thus strongly affect the scattering and the overall energy distribution of an explosion source. To assess the effects of shallow attenuation, we use a set of model parameters similar to that used for testing the effect of rms fluctuations, except we replace the infinite Q in the top 10 km with QP 100 and QS 50. These parameters are listed in group 4 of Table 3, and the results are shown in Figure 11. Comparing Figure 11 to Figure 5, two prominent features can be identified. First, compared to the purely elastic case, there is significant energy loss in the model with intrinsic attenuation. For example, the maximum amplitude (square-root energy) drops approximately 35% for direct Rg waves, 35% for scattered Rg waves, and 40% for Lg waves. Second, the short-period waves experience even larger attenuation than long-period waves. This is especially true for the scattered Rg waves and the Lg waves. By using low Q values in the calculation, our results should give fairly extreme characterization of the effect of attenuation.

Figure 8.

Similar to Figure 4 except each panel is for a model with different correlation distance.

0.06 0.05 0.04

0.08 0.07 0.06 0.04 0.03 0.02 0.01 0.00 0.4 0.6 0.8 1.0 2.0 4.0 6.0 8.0 10. 0.05
0.02 1.0 1.5 0.4 0.6 0.8 1.0 2.0 4.0 6.0 8.0 10. 0.01 2.0 0.00 2.5
0.002 1.0 1.5 0.001 2.0 0.00 2.5 0.4 0.6 0.8 1.0 2.0 4.0 6.0 8.0 10.

nc 3.03.5 yH 4.0 z

Fr 2.0 eq 2.5 ue 3.0 nc y H 3.5 4.0 z
-0.025 -0.02 -0.015 -0.01 -0.005 0.00 0.005 0.01 0.015 0.02 0.4 0.6 0.8 1.0 2.0 4.0 6.0 8.0 10.
0.006 0.005 0.004 0.003 0.002 0.001 0.00 0.4 0.6 0.8 1.0 2.0 4.0 6.0 8.0 10.
0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 -0.01 0.4 0.6 0.8 1.0 2.0 4.0 6.0 8.0 10.

-0.001

Figure 9.
Similar to Figure 5, except responses are functions of frequencies and correlation lengths.
Variance of Statistical Results To obtain statistical relationships between the source/ model parameters and the energy-partitioning process, we average the measurements from individual realizations. To
investigate the variance from a group of realizations with the same statistical parameters, we compare results for models with an rms value of 0.15 km, a correlation length of 0.5 km, a source depth of 0.5 km, and an infinite Q. The response functions for different phases are presented in Figure 12ac, where the horizontal coordinates are frequency and different rows are for different realizations. The response curves for different realizations, their mean values, and the standard deviations are also shown in Figure 12df. The primary Rg wave is relatively stable. The scattered Rg wave has large variations at frequencies below about 2 Hz, while the Lg wave shows large variations at frequencies above 2 Hz. Although based on the same statistical parameters, response functions from different realizations show variations in amplitudes and local minima, suggesting that the partitioning is partially affected by deterministic features very close to the source.

(E/E0)1/2

Discussions and Conclusion

scale factor ka.

Figure 10.
Net scattered Lg energy as a function of normalized
We have used the 2D P-SV boundary-element simulation (Ge et al., 2005) and the embedded-array local-slowness analysis method (Xie, Ge, and Lay, 2005) to investigate the effect of topographic scattering on the near-source energy
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
0.02 1.0 1.5 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.01 2.0 2.5 0.00
m sk rm (d) scattered loss of Rg energy
sk rm (e) net scattered Rg energy
m sk rm (f) scattered gain of Lg energy
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00

2.0 2.5

-0.01 1.0 1.5 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.00 0.01 2.0 2.5 0.02 3.0 3.5
0.003 0.002 0.001 2.0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 2.5 0.00 3.0 3.5

Figure 11.

Similar to that shown in Figure 4, except a low Q top layer is used in the simulation. To facilitate comparison, the same vertical scale as in Figure 4 is used here.
partitioning for an explosion source. Random topographic models with different statistical properties, variable source depth, and different Q models are investigated using numerical simulations. The responses of different phases as functions of frequency and source/model parameters are calculated and their energy budget evaluated. The source spectrum has been corrected from these response functions (following equations 1 and 2). To compare the result with the expected phase-amplitude spectrum, we can use the square root of the energy and scale the result with a factor f 1. The results reveal that free-surface scattering has strong effect on near-source energy partitioning. The scattering process can excite the Rg wave for a moderately rugged topography, but dramatically attenuates short-period Rg waves when the surface becomes too rugged. For models with a high velocity shallow crust, the free-surface scattering provides an important mechanism that transfers energy for an explosion source into the Lg wave in the near-source region. At lower frequencies and for a moderately rugged free surface, the Rg-to-Lg transfer is relatively efficient. At higher frequencies and for a very rugged free surface, the bodyto-Lg transfer may dominate the process. The correlation length of the random free-surface fluctuation provides spe-
cific frequency dependence to the transfer function, with maximum coupling near ka 1. Intrinsic attenuation within the uppermost crust has a strong effect on the energy transfer through surface scattering, with high-frequency content losing energy faster than the lower-frequency waves. Recent observations demonstrate that explosiongenerated shear waves have a corner frequency that scales approximately as V S =RE (Xie, 2002; Fisk, 2006), where RE is the elastic radius and V S is the source shear velocity. These results indicate the importance of very near-source processes for the P- and S-wave energy partitioning. The results here show that several mechanisms affect the frequency dependence of the energy-response function. The Rg excitation function, source depth, correlation length, and intrinsic attenuation all impart frequency dependence to the energy transfer functions, either directly or indirectly. Surface scattering (especially Rg-to-Lg) as a mechanism for generating regional-phase Lg signals for explosions in high velocity crust has been proposed and supported by many studies (Gupta et al., 1992; Patton and Taylor, 1995; Gupta et al.,1997; Bonner et al., 2003; Gupta et al., 2005; Myers et al., 2005; Stevens et al., 2005). However, there is some controversy as to whether the Rg-to-Lg scatter-

model the Rg-to-Lg scattering. They suggested that the observed Rg decay rate can be used to constrain the calculation and the result can be compared with regional Lg observations. As a comparison, the method developed here puts the rough topography explicitly into the model and directly calculates the conversion process. Patton and Taylor (1995) and Gupta et al. (1997) investigated the effects of spall on the amplitude spectra of regional seismic signals. They suggested that Lg is generated by near-source scattering of the Rg wave into trapped body waves, and a best-fit spectrum can be explained by an explosion, together with a spall source. Stevens et al. (2005) investigated a group of former Soviet explosions and performed nonlinear source calculations to interpret the observations. They found that the nonlinear calculations can be matched fairly well using linear calculations for a point explosion plus a compensated linear vector dipole (CLVD) source with half the explosion moment, with the CLVD component decreasing slowly with the increasing depth. A CLVD source will efficiently generate S waves, which can be trapped in the crust, and it will also increase Rg wave excitation, which is then available to be scattered into Lg energy through surface scattering. In this article, we focus on the role of the rough free surface in scattering of the surface and body waves into Lg energy, and we have not yet incorporated a composite source model in our simulation. However, any source that intrinsically generates more S waves and Rg waves, should increase the Lg energy.
Explosion-Source Energy Partitioning and Lg-Wave Excitation The numerical simulation and slowness-analysis methods used in this study are for 2D models. The 2D calculation affects the amplitude decay due to geometric spreading. More importantly, scattering in 2D geometry is not the same as in the 3D case. Generally speaking, near-source scattering will probably be more important for 3D cases, and partitioning of energy into the SH component can only be modeled in 3D. Considering these issues, the results presented in this article should be viewed as qualitative, documenting the basic nature of the surface scattering effects, but not the full energy partitioning. Expanding the current analysis to fully 3D models is currently ongoing and will be presented in a future publication.
Frankel, A. (1989). A review of numerical experiments on seismic wave scattering, Pure Appl. Geophys. 131, 639685. Ge, Z., L. Y. Fu, and R. S. Wu (2005). P-SV wavefield connection technique for regional wave propagation simulation, Bull. Seismol. Soc. Am. 95, 13751386. Gupta, I. N., W. W. Chan, and R. A. Wagner (1992). A comparison of regional phases from underground nuclear explosions at East Kazakh and Nevada test sites, Bull. Seismol. Soc. Am. 82, 352382. Gupta, I. N., W. W. Chan, and R. A. Wagner (2005). Regional source discrimination of small events based on the use of Lg wavetrain, Bull. Seismol. Soc. Am. 95, 341346. Gupta, I. N., T. Zhang, and R. A. Wagner (1997). Low-frequency Lg from NTS and Kazakh nuclear explosions: observations and interpretations, Bull. Seismol. Soc. Am. 87, 11151125. Hartse, H. E., S. R. Taylor, W. S. Phillips, and G. E. Randall (1997). A preliminary study of regional seismic discrimination in central Asia with emphasis on western China, Bull. Seismol. Soc. Am. 87, 551568. Hong, T. K., and J. Xie (2005). Phase composition of regional seismic waves from underground nuclear explosions, J. Geophys. Res. 110, B1, 2303, doi 10.1029/2005JB003753. Kim, W. Y., V. Aharonian, A. L. Lerner-Lam, and P. G. Richards (1997). Discrimination of earthquakes and explosions in southern Russia using regional high-frequency three-component data from IRIS/JSP Caucasus Network, Bull. Seismol. Soc. Am. 87, 569588. Kim, W. Y., D. W. Simpson, and P. G. Richards (1993). Discrimination of earthquakes and explosions in the eastern United States using regional high-frequency data, Geophys. Res. Lett. 20, 15071510. McLaughlin, K. L., and R. S. Jih (1988). Scattering from near-source topography: teleseismic observations and numerical simulations, Bull. Seismol. Soc. Am. 78, 13991414. Myers, S. C., J. Wagoner, L. Preston, K. Smith, and S. Larsen (2005). The effect of realistic geologic heterogeneity on local and regional P=S amplitude ratios based on numerical simulations, in Proc. of the 27th Seismic Res. Rev.: Ground-Based Nuclear Explosion Monitoring Technologies, 123132. Myers, S. C., W. R. Walter, K. Mayeda, and L. Glenn (1999). Observations in support of Rg scattering as a source for explosion S waves: regional and local recordings of the 1997 Kazakhstan depth of burial experiment, Bull. Seismol. Soc. Am. 89, 544549. Nuttli, O. W. (1986). Yield estimates of Nevada test site explosions obtained from Lg waves, J. Geophys. Res. 91, 21372151. Patton, H. J. (2001). Regional magnitude scaling, transportability, and Ms: mb discrimination at small magnitudes, Pure Appl. Geophys. 158, 19512015. Patton, H. J., and S. R. Taylor (1995). Analysis of Lg spectral ratios from NTS explosions: implications for the source mechanisms of spall and the generation of Lg waves, Bull. Seismol. Soc. Am. 85, 220236. Stevens, J. L., G. E. Baker, H. Xu, T. J. Bennett, N. Rimer, and S. M. Day (2003). The physical basis of Lg generation by explosion sources, in Proc. of the 25th Seismic Res. Rev., Nuclear Explosion Monitoring: Building the Knowledge Base, 456465. Stevens, J. L., G. E. Baker, H. Xu, and T. J. Bennett (2005). The physical basis of the explosion source and generation of regional seismic phases, in Proc. of the 27th Seismic Res. Rev.: Ground-Based Nuclear Explosion Monitoring Technologies, 663672. Taylor, S. R. (1996). Analysis of high frequency Pn=Lg ratios from NTS explosions and western U.S. earthquakes, Bull. Seismol. Soc. Am. 86, 10421053. Taylor, S. R., and H. E. Hartse (1997). An evaluation of generalized likelihood ratio outlier detection to identification of seismic events in western China, Bull. Seismol. Soc. Am. 87, 824831. Taylor, S. R., M. D. Denny, E. S. Vergino, and R. E. Glaser (1989). Regional discrimination between NTS explosions and western U.S. earthquakes, Bull. Seismol. Soc. Am. 79, 11421176. Vogfjrd, K. S. (1997). Effects of explosion depth and earth structure on excitation of Lg waves: S revistited, Bull. Seismol. Soc. Am. 87, 11001114.

JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 115, B08307, doi:10.1029/2009JB007157, 2010
Seismic Lgwave Q tomography in and around Northeast China
LianFeng Zhao,1 XiaoBi Xie,2 WeiMin Wang,3 JinHai Zhang,1 and ZhenXing Yao1
Received 20 November 2009; revised 8 February 2010; accepted 3 March 2010; published 14 August 2010.
[1] We investigate regional variations in the Lgwave quality factor (Q) in Northeast China and its vicinity with a tomographic method. Digital seismic data recorded at 20 broadband stations from 125 regional events are used to extract Lgwave spectra. Tomographic inversions are independently conducted at 58 discrete frequencies distributed log evenly between 0.05 and 10.0 Hz. We simultaneously invert for the Lgwave Q distribution and source spectra at individual frequencies without using any a priori assumption about the frequency dependence of the Q model and source function. The best spatial resolution is approximately in wellcovered areas for frequencies between 0.4 and 2.0 Hz. The Lg Q shows significant regional variations and an apparent relationship with regional geology. We use a statistical method to investigate the regional variations of Lg Q and their frequency dependence. The average Q0 (1 Hz Lg Q) in the entire investigated region is 414. Sedimentary basins are usually characterized by lower average Q0 values (from 155 to 391), while volcanic mountain areas have relatively high average Q0 values (from 630 to 675). Lg Q generally increases with increasing frequency. However, the frequency dependence has complex nonlinear features on a doublelogarithmic scale, indicating that the commonly adopted powerlaw relationship may be oversimplified in a broad frequency band. The frequency dependence varies in different geological areas, with larger variations seen at lower frequencies.
Citation: Zhao, L.F., X.B. Xie, W.M. Wang, J.H. Zhang, and Z.X. Yao (2010), Seismic Lgwave Q tomography in and around Northeast China, J. Geophys. Res., 115, B08307, doi:10.1029/2009JB007157.

1. Introduction

[2] The Lg wave is one of the most prominent seismic phases in highfrequency seismograms observed over continental paths at regional to teleseismic distances. It is usually understood to be a sum of higher mode surface waves propagating in the crustal waveguide or multiply reflected S waves supercritically bouncing between the free surface and the Moho discontinuity [Knopoff et al., 1973; Bouchon, 1982; Kennett, 1984; Xie and Lay, 1994]. Thus, the Lg wave samples the crust waveguide relatively evenly. It is also sensitive to the characteristics of the free surface and the Moho discontinuity and to the crustal thickness. Given these attributes, much attention has been paid to Lgwave data for investigating the properties of the crust. [3] The Lgwave quality factor Q (QLg) describes the attenuation of Lg signals and is one of the basic parameters useful for characterizing the Earths crust. Scattering losses and intrinsic attenuation are both responsible for the ampli1 Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, China. 2 Institute of Geophysics and Planetary Physics, University of California, Santa Cruz, Santa Cruz, California, USA. 3 Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing, China.
Copyright 2010 by the American Geophysical Union. 01480227/10/2009JB007157
tude decay of Lg waves. Lgwave scattering is closely related to the heterogeneities of all scales in the crustal waveguide. Thus, Lg waves and Lg coda have been used successfully for the past three decades to infer and measure the attenuation and scattering structure of the crust [Nuttli, 1973; KadinskyCade et al., 1981; Xie and Nuttli, 1988; McNamara et al., 1996; Mitchell et al., 1997; Rodgers et al., 1997; Cong and Mitchell, 1998; Mellors et al., 1999; Philips et al., 2000; Wu et al., 2000; Sandvol et al., 2001; Fan and Lay, 2002, 2003a, b; Ottemller, 2002; Ottemller et al., 2002; Xie et al., 2004, 2006; Phillips et al., 2005; Zor et al., 2007; Mitchell et al., 2008; Phillips and Stead, 2008]. QLg depends on the types of material, thermal status, and degree of heterogeneities in the crust. In general, high QLg values correlate well with stable ancient crust, while relatively low QLg values correlate well with recently deformed crust and active tectonic environments [e.g., Fan and Lay, 2002, 2003a, b; Xie et al., 2004, 2006; Phillips et al., 2005]. Several studies have also suggested that crustal thickness or undulation of the Moho discontinuity can result in significant Lgwave attenuation and that both the oceanic crust and the sedimentary basins can strongly attenuate or even block Lg waves by disrupting their underlying mode structures [e.g., Kennett, 1986; Campillo, 1987; Campillo et al., 1993; Zhang and Lay, 1995; Shapiro et al., 1996; Wu et al., 2000]. [4] Northeast China is a complex convergence zone between several largescale geological units including the West Pacific tectonic belt, the Siberian plate, the North China

a CNDSN, China National Digital Seismic Network; GSN, Global Seismographic Network.
craton, and the Yangzi plate. During the Cenozoic, tectonic events frequently occurred in this region [Chi, 1988; Cai et al., 2002; Ma et al., 2003; Ge and Ma, 2007]. Figure 1a shows a geological map of Northeast China including the major faults, folding belts, and basins. Figure 1b illustrates the investigated area (highlighted by colored topography), with the two major mountain areas, the Changbaishan Mountains (CMs) and Daxinganling Mountains (DMs), labeled. [5] Investigating the QLg distribution in Northeast China and its vicinity can provide new insight into the regional geology. Earlier works addressing Lgwave attenuation in this region were conducted by Jin and Aki [1988], who obtained a 1 Hz coda Q map in China. By using analog records Ge et al. [1989] and Huang et al. [1990] obtained an average QLg of 500 in Northeast China. Mitchell et al. [1997, 2008] investigated the Lg coda Q variations across Eurasia and obtained a largescale 1 Hz coda Q image. Using the powerlaw attenuation model Q( f ) = Q0f h, where Q0 is the 1 Hz Q and h is an index, Xie et al. [2006] measured the QLg in eastern Eurasia and derived a tomographic model for Q0 with a resolution of between 4 and 10. They revealed the relation between the QLg distribution and the deformed regions of the regional fault systems. Phillips et al. [2005] used the amplitude ratio technique to image the QLg in central and eastern Asia. On average, their image is resolved to 2.5, with the resolution peaking at 1.5 in the best covered areas. Using the ML amplitude catalog, Pei et al. [2006] obtained a laterally varying 1 Hz Q model within the crust of northern China. Chung et al. [2007] studied the 1 Hz QLg around the Korean Peninsula including part of Northeast China and the Sea of Japan. Ford et al. [2009] compared the QLg images in the Yellow Sea and North Korea region obtained using different methods. In this study, we use a larger regional data set to constrain the QLg and focus our attention on Northeast China and its vicinity. Without applying any a priori assumptions on the frequency dependence of the attenuation model or source spectra, we invert

3. Methodology

[9] Commonly used methods for Lgwave Q inversion include simultaneous inversion of source and attenuation, the twostation or reversed twostation method, and the source pair/receiver pair method. For Lg coda waves the coda normalization and codasource normalization [e.g., Ford et al.,

3 of 22

Table 2. Parameters of the Events Used in This Study
Epicentral parameter Contributor ISC ISC ISC ISC ISC ISC ISC ISC ISC ISC ISC ISC BJIa ISC ISC ISC ISC ISC BJI BJI BJI NEIC ISC NEIC NEIC NEIC NEIC SKHLb NEIC ISC ISC ISC ISC ISC BJI BJI ISC ISC ISC ISC ISC BJI ISC NEIC NEIC BJI BJI ISC ISC BJI ISC ISC ISC BJI SKHL ISC ISC ISC ISC BJI ISC BJI BJI BJI ISC BJI ISC ISC ISC Date (yyyy/mm/dd) 1995/10/05 1998/04/14 1998/07/24 1998/07/24 1998/09/26 1998/11/11 1999/01/20 1999/01/20 1999/01/20 1999/04/01 1999/08/13 1999/08/14 1999/09/03 1999/12/27 2000/01/11 2000/01/12 2000/05/14 2000/05/23 2000/11/08 2000/11/14 2000/12/13 2001/03/02 2001/04/20 2001/05/21 2001/06/05 2001/08/25 2001/09/01 2001/09/03 2001/09/19 2001/10/05 2002/01/27 2002/04/16 2002/04/21 2002/04/30 2002/05/05 2002/06/13 2002/06/16 2002/08/29 2002/10/20 2002/10/20 2002/10/21 2003/03/30 2003/03/30 2003/04/23 2003/04/23 2003/05/03 2003/05/07 2003/05/22 2003/06/01 2003/06/05 2003/06/05 2003/06/08 2003/06/10 2003/06/12 2003/06/14 2003/06/18 2003/06/21 2003/07/02 2003/07/18 2003/08/03 2003/08/16 2003/08/16 2003/08/17 2003/08/17 2003/08/25 2003/09/09 2003/09/16 2003/10/07 2003/10/09 Time (UTC) 22:26:55.39 02:47:50.01 23:12:16.10 23:19:39.92 14:16:43.93 17:27:55.36 13:39:52.10 13:27:58.02 12:36:47.21 01:59:47.39 18:36:22.38 00:04:37.69 11:23:32.70 11:27:19.01 23:43:56.00 05:00:37.10 15:48:50.23 23:44:36.59 13:06:37.10 08:41:52.50 09:59:31.20 04:29:06.26 04:35:37.75 15:35:12.75 14:56:45.61 13:21:26.11 13:08:11.90 03:09:34.00 08:07:26.07 09:20:57.05 18:20:31.76 22:52:38.63 19:34:38.61 03:22:43.56 10:02:24.90 07:36:03.70 21:58:40.16 18:32:16.26 15:46:19.95 15:52:10.54 00:10:21.67 11:00:44.40 11:10:55.46 13:46:08.32 18:39:19.17 12:35:24.70 11:34:16.50 08:48:46.88 02:49:17.28 10:35:55.40 23:18:42.55 01:56:48.46 03:23:20.18 09:33:59.40 14:10:09.70 14:24:31.50 12:13:09.36 07:44:08.09 14:03:06.81 23:49:10.00 10:58:40.76 11:59:51.30 05:24:23.70 16:33:22.50 18:48:32.14 05:09:29.60 11:24:53.67 15:27:27.80 15:53:29.88 Latitude (N) 39.727 39.595 48.889 48.903 42.296 48.189 57.517 57.569 57.513 39.650 48.468 48.114 48.960 40.541 40.546 40.676 48.982 40.645 45.590 39.150 40.850 40.496 48.289 36.817 40.665 35.047 47.360 48.260 38.043 45.095 53.681 40.658 37.374 40.692 40.060 49.260 40.601 49.451 44.600 44.914 44.605 42.050 37.625 39.894 39.532 42.100 37.870 39.446 49.800 36.300 36.477 35.409 40.732 36.360 49.110 47.533 39.691 36.995 53.953 41.350 43.813 43.770 44.040 37.930 38.642 44.060 56.075 45.240 41.435 Longitude (E) 118.531 118.577 131.484 131.241 123.668 133.145 120.479 120.476 120.601 125.161 128.538 128.236 130.410 123.070 123.095 122.850 129.924 122.855 118.140 125.370 125.530 115.033 117.117 106.507 108.322 135.654 142.523 133.150 119.531 105.585 125.376 128.652 114.697 122.964 127.240 122.940 123.007 123.067 117.470 117.142 117.801 123.630 123.856 117.328 117.705 124.040 121.450 118.060 130.794 120.580 119.940 111.605 111.340 120.370 131.830 116.985 118.368 103.944 134.323 130.580 119.658 119.670 119.800 120.610 112.528 119.350 111.296 133.597 125.982 Depth (km) 8.6 23.0 14.3 11.8 10.0 9.0 28.6 33.9 30.8 10.4 8.0 44.5 29.0 21.0 5.0 36.3 5.5 30.6 31.0 15.0 27.0 33.0 10.0 33.0 33.0 23.2 10.0 13.0 33.0 32.0 16.0 10.0 24.0 15.0 15.0 20.0 38.7 16.6 33.0 10.0 10.0 28.0 10.0 33.0 47.5 21.0 10.0 10.0 7.0 28.0 33.0 14.0 22.8 15.0 7.0 16.2 16.0 34.4 10.0 15.0 8.8 16.0 20.0 17.0 10.0 19.0 18.3 49.4 8.0 Magnitude (mb) 4.0 4.5 4.6 4.9 5.1 4.2 4.6 4.3 5.1 3.9 4.1 4.2 4.3 3.7 4.1 4.4 4.5 3.8 4.2 3.5 4.0 4.1 4.1 4.4 4.6 4.7 4.8 4.2 4.6 4.8 4.9 5.0 3.8 4.1 4.5 4.7 3.4 3.5 4.1 4.4 4.5 4.6 3.7 3.6 4.5 4.7 4.8 4.9 5.1 4.5 3.7 4.1 4.2 4.4 3.9 4.1 4.3 4.1 3.5 3.6 4.3 4.4 5.1 5.7 4.2 4.4 4.1 4.3 4.5 4.2 4.1 4.8 5.0 3.7 3.9 4.0 3.5 3.6 4.6 4.8 3.8 4.2 3.4 4.3 4.7 4.1 4.3 3.6 3.7 3.8 4.0 3.7 3.8 3.7 4.3 4.0 3.9 4.4 3.7 4.3 4.5 4.6 3.9 4.0 5.4 5.5 5.8 3.8 4.1 5.2 5.3 4.4 4.5 Inverted source parameter Seismic moment M0(N m) 7.36 .37E+15 4.67 .24E+15 1.62 .23E+15 2.78 .25E+15 3.57 .41E+14 9.84 1.27E+14 2.09 .26E+15 5.50 .31E+15 2.30 .14E+15 1.33 .19E+14 3.05 .20E+15 1.64 .23E+14 5.50 .67E+14 5.17 .45E+14 1.01 .12E+16 9.84 .89E+14 1.10 .12E+15 6.74 .98E+14 2.31 .42E+14 1.13 .15E+14 3.71 .41E+13 1.06 .05E+15 1.71 .29E+14 1.86 .17E+15 3.10 .18E+15 4.00 1.72E+16 6.96 1.37E+16 2.57 .40E+14 1.66 .15E+15 1.97 .26E+15 2.25 .49E+14 1.17 .09E+15 5.23 .59E+15 4.27 .26E+14 4.00 .52E+13 9.83 1.02E+13 3.50 .26E+14 6.21 .60E+14 4.08 .41E+15 1.51 .15E+15 2.61 .35E+14 2.24 .23E+14 4.70 .39E+15 3.27 .21E+14 7.44 .58E+14 1.79 .22E+14 2.43 1.11E+14 2.79 .22E+14 5.10 .51E+14 6.75 .43E+14 1.16 .33E+14 1.31 .22E+14 7.67 .70E+14 5.37 1.38E+14 2.13 .16E+14 1.26 .07E+15 3.14 .12E+14 5.24 .66E+15 2.12 .14E+15 1.54 .16E+14 5.23 .49E+16 1.50 1.01E+14 1.26 .08E+14 2.46 .25E+14 5.56 .60E+14 1.67 .17E+14 1.23 .23E+17 3.20 .22E+15 1.65 .23E+14 Corner frequency fc (Hz) 1.08 1.01 1.17 1.01.88.92.64.50.55 1.92.96 2.22 1.13.91.81.93 1.42 1.43 1.93 1.55 2.71 1.86 1.88 1.17 1.39.18.39.92 1.38 1.50.78.73 1.20 1.06 1.74 1.69 1.40.97 1.08 1.28 1.76 1.69.79 2.61 1.46.78.73 2.39 1.27.94 1.39 1.55 2.14.75.86 1.52 2.12.76.80 1.38.68 1.23 2.43.89 1.82 1.73.32.54 1.60 .07.07.12.10.11.08.08.03.04.31.09.36.11.09.09.08.17.19.23.22.29.09.23.11.11.03.06.12.14.20.11.05.11.07.21.15.12.07.09.12.23.19.07.24.09.06.21.18.11.06.28.23.29.11.05.13.14.06.06.11.03.56.18.07.20.22.04.04.18

4 of 22

B08307 Table 2. (continued)
Epicentral parameter Contributor BJI ISC ISC ISC ISC CENCc ISC ISC ISC ISC ISC ISC BYKLd ISC ISC NEIC NEIC NEIC NEIC ISC NEIC BYKL NEIC NEIC NEIC NEIC ISC CENC NEIC NEIC NEIC NEIC NEIC CENC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC NEIC
Inverted source parameter Longitude (E) 125.650 134.389 119.820 100.956 101.058 103.930 118.704 111.662 129.032 126.702 118.982 128.883 118.640 118.209 125.723 111.714 111.643 129.928 142.120 117.944 131.727 117.930 127.979 117.753 122.584 107.909 141.879 104.450 121.966 131.473 109.017 109.260 125.002 125.030 104.144 121.003 120.537 120.751 116.413 124.122 115.591 121.049 118.153 118.469 116.153 102.352 141.908 129.094 119.558 120.726 110.075 110.300 141.733 141.722 141.798 113.381 Depth (km) 15.0 10.0 10.0 9.0 0.8 10.0 10.0 10.0 10.0 9.0 23.4 9.0 16.0 18.0 10.0 10.0 10.0 29.2 13.7 4.6 10.0 12.0 10.0 13.4 10.0 10.0 10.0 10.0 10.0 10.0 6.0 17.8 10.0 25.0 16.0 10.0 6.0 20.0 10.0 10.0 10.0 35.4 10.0 10.0 10.0 10.0 14.0 0.0 10.0 10.0 10.0 10.0 10.0 10.0 5.0 19.0 Magnitude (mb) 4.3 4.4 4.5 3.8 4.3 5.0 5.1 5.2 4.0 4.6 4.8 4.9 4.0 4.2 4.3 4.4 5.2 3.8 4.6 4.5 4.0 4.2 4.8 4.4 5.5 5.6 5.7 5.8 4.0 4.1 4.3 4.9 4.6 4.1 5.3 4.9 4.6 4.7 4.8 3.9 4.0 4.1 4.4 4.0 5.2 4.5 4.0 4.9 5.2 5.4 4.9 4.0 4.3 5.5 5.1 5.0 5.2 4.9 5.8 4.2 5.6 5.7 4.6 4.9 4.4 3.7 3.9 4.1 5.0 5.1 6.0 3.6 4.2 4.3 4.7 4.8 5.1 5.2 5.3 5.3 5.6 5.3 5.0 Seismic moment M0(N m) 2.82 .65E+14 1.72 .15E+15 6.87 .46E+14 4.79 .19E+17 7.51 2.90E+16 8.30 1.95E+16 4.21 .21E+14 5.44 .57E+15 1.21 .12E+16 3.01 .42E+14 2.04 .14E+15 2.84 .25E+15 2.05 .36E+14 5.17 .34E+16 9.05 .54E+14 4.30 .46E+15 6.20 .97E+14 5.30 .75E+16 7.45 .66E+15 8.66 1.28E+15 3.64 .31E+14 3.89 .28E+14 6.57 .56E+14 3.39 .54E+16 1.12 .08E+15 1.13 .11E+15 2.00 .83E+17 1.80 .34E+16 5.29 .33E+14 8.99 .66E+14 4.12 .43E+16 1.89 .16E+16 1.07 .13E+16 3.69 5.73E+14 4.89 1.24E+16 3.12 .20E+15 4.01 .54E+17 1.78 .41E+17 2.70 .22E+15 5.62 .60E+15 3.12 .31E+15 7.87 .44E+14 5.87 .29E+14 9.03 .44E+14 6.29 .40E+15 3.03 .40E+16 1.55 .14E+17 6.71 .81E+14 1.41 .10E+15 2.93 .35E+15 1.02 .25E+17 1.09 .17E+17 3.40 .25E+16 2.80 .46E+17 4.41 .45E+17 1.28 .21E+16 Corner frequency fc (Hz).82.95 1.67.29.24.23 2.16 1.00.62 1.09 1.25.82.91.69 1.09.62 1.08.29.73.63.90.80 1.15.48 1.55.55.14.40 1.42 1.06.67.72.91 1.02.40 1.27.23.29 1.27 1.00 1.23 1.72 1.45 1.54 1.04.50.27.60 1.27.48.25.27.44.15.13.66 .13.11.11.02.05.03.14.09.08.12.10.08.12.03.09.06.17.03.07.08.06.05.10.07.13.04.03.05.08.11.05.05.09.47.07.10.02.04.12.07.14.17.07.08.06.05.03.05.10.04.04.04.02.02.01.08
Date (yyyy/mm/dd) 2003/10/10 2003/10/16 2003/10/17 2003/10/25 2003/10/25 2003/11/13 2003/11/14 2003/11/25 2004/01/16 2004/01/17 2004/01/20 2004/01/25 2004/03/22 2004/03/24 2004/03/24 2004/05/26 2004/05/27 2004/05/29 2004/05/30 2004/06/28 2004/09/16 2004/10/24 2004/12/16 2005/01/02 2005/01/24 2005/02/27 2005/03/08 2005/04/09 2005/05/09 2005/07/06 2005/07/20 2005/07/20 2005/07/25 2005/07/25 2005/07/26 2005/09/19 2005/11/10 2005/12/11 2006/01/06 2006/03/31 2006/04/09 2006/05/03 2006/05/03 2006/05/03 2006/07/04 2006/07/07 2006/08/17 2006/10/09 2006/11/03 2006/11/20 2006/12/04 2007/07/04 2007/08/02 2007/08/02 2007/08/02 2007/08/23

Time (UTC) 13:34:27.80 21:19:20.01 01:38:30.75 12:41:35.20 13:25:21.87 02:35:10.00 21:43:18.57 05:40:30.54 19:08:31.92 05:09:30.54 08:34:12.05 19:46:27.44 03:49:59.60 01:53:47.50 19:55:49.18 23:56:52.11 00:36:42.74 10:14:28.44 02:52:12.28 14:22:44.47 17:14:37.47 04:31:38.60 18:59:14.60 00:24:37.88 12:22:46.53 13:27:36.57 23:58:36.94 02:44:18.90 11:02:24.21 23:10:16.68 21:54:05.72 18:06:57.44 15:43:36.30 15:57:14.80 12:16:08.93 03:27:53.20 19:29:54.00 15:54:15.90 01:56:38.94 12:23:17.86 09:23:59.66 00:26:37.50 13:53:42.61 14:02:25.93 03:56:26.94 14:12:07.14 15:20:35.02 01:35:28.00 06:21:39.28 00:07:28.37 09:14:04.60 01:23:24.44 08:06:28.81 05:22:17.69 02:37:42.38 04:49:19.80
Latitude (N) 41.430 53.927 43.534 38.396 38.340 34.780 39.887 36.164 53.184 54.286 39.858 53.190 56.680 45.349 54.316 54.124 54.106 36.634 47.311 56.592 45.137 56.630 41.804 56.723 51.755 40.738 52.163 34.040 37.686 48.295 43.066 43.045 46.892 47.140 42.541 49.878 57.445 57.431 51.716 44.624 35.752 48.783 39.991 39.690 39.071 44.551 46.542 41.294 43.469 57.287 55.769 55.474 46.743 46.714 47.116 55.950
Beijing regional network. Sakhalin regional network. China earthquake network center. d Baykal regional network.
2008] methods are often used. Limited by the available data, tradeoffs are often adopted to simplify the QLg inversion. For example, Xie [1993] addressed the tradeoffs among QLg, the source radiation pattern, and the site response. Ottemller [2002] and Ottemller et al. [2002] imaged the
apparent crustal QLg in Central America based on similar considerations. However, many researchers [e.g., Phillips et al., 2005] have emphasized the importance of source coupling and site response on attenuation measurement. Here we apply the tomographic method to investigate QLg

5 of 22

Figure 2. Statistics of the regional data set, with (a) number of rays versus distances and (b) number of records per event. s M P mean square method, Af 1=M 1 Afi 2 , where A( fi) is the spectrum directly from the Fourier transform, fi 2[ f1, f2], f1, and f2 can be determined by log10( f ) log10( f1) = log10( f2) log10( f ) = 0.02, and M is the number of data points between f1 and f2. Figure 5e shows 58 discrete amplitude values calculated between 0.05 and 10.0 Hz for both Lg waves and the noise series. [14] To conduct the noise correction we assume that the recorded seismogram in the Lg window is a superposition of the Lg wave and noise. By further assuming that the two are uncorrelated, we can obtain the Lgwave spectral amplitude from the raw data and the noise amplitude using A2( f ) = S A2 ( f ) A2 ( f ), where A is the amplitude and the subscripts O N S, O and N denote the true signal, the observed raw data, and the noise, respectively. From spectral amplitudes of the signal and the noise we can calculate signaltonoise ratios at individual frequencies, which are shown in Figure 5f. By introducing a threshold of 2.0 for the signaltonoise ratio (dashed line in Figure 5f), we can judge the quality of the Lgwave spectrum. Illustrated in Figure 5g is the noise corrected Lgwave spectrum where points below the threshold have been dropped. [15] Owing to the strong attenuation of highfrequency signals in the crustal waveguide, the amplitude of the high frequency Lg decays rapidly and tends to be affected by noise at large distances. To examine the quality of these data, amplitudedistance curves for individual frequencies are calculated. From these curves we find that spectral amplitudes beyond certain epicentral distances are sometimes suspicious. To eliminate their effect we set up a truncation distance depending on the frequency (Table 3). For a given frequency the data collected from stations beyond the truncation distance are excluded. 3.2. Lg Q Tomographic Scheme [16] We design a tomographic method to simultaneously invert the QLg and the source function frequency by frequency. This method is similar to those used in the seismic velocity tomography and is reported in the Appendix. To construct the inversion system we assume that the Lg wave propagates in the crustal waveguide along the great circle path between the source and the station. Additionally, we assume that the QLg can be expressed as Q = Q(x,y, f ), where (x,y) is the surface location. The inversions are independently conducted for individual frequencies. We do

in Northeast China and the surrounding area. Limited by the available data, the frequencydependent site response and source coupling are not considered in our inversion. For the same reason, we neglect the effect of radiation pattern and assume that the source is isotropic. 3.1. Data Preprocessing [10] Following Xie and Mitchell [1990], we conduct the following data preprocessing for all available Lg waveforms: (1) setting windows for the Lg waves and the noise, (2) calculating Fourier spectra for the Lg waves and the noise, (3) sampling the spectral amplitudes, and (4) correcting the effect of noise on the signal amplitudes. [11] Xie and Mitchell [1990] suggested a recursive zero phase, fourthorder Butterworth filter with a 1 Hz corner frequency to isolate the corresponding Lgwave train and used it to determine the Lgwave groupvelocity window. Nuttli [1986] determined the Lg window directly from the World Wide Standard Seismogram Network (WWSSN) shortperiod (SP) seismograms. Analogously, Patton and Schlittenhardt [2005] convolved the broadband waveforms with a WWSSN response to simulate shortperiod seismograms from which they visually selected the Lg groupvelocity window. Both methods emphasize the highfrequency content in the Lg wave train. Here we follow Patton and Schlittenhardt [2005]. Figure 4 illustrates simulated WWSSN SP seismograms for the same event as shown in Figure 3. After visually inspecting all of the seismograms, we set the Lg groupvelocity window to be either 3.63.0 or 3.52.9 km/s. The noise series are picked from a time window that has the same length as the Lg window and is located immediately before the first arriving P wave [Zhao et al., 2008]. [12] To calculate Fourier spectra for both the Lg wave and preP noise, we add 10% duration intervals both before and after the time windows and apply cosine tapers on the extended portions. A fast Fourier transform (FFT) is then performed on the windowed signals to calculate the spectra. Figure 5 illustrates this process, where Figures 5a and 5b are broadband seismograms with time windows for the Lg wave and preP noise, Figures 5c and 5d are windowed wave trains for the Lg wave and the noise series, and Figure 5e shows the amplitude spectra for the Lg wave and ambient noise. [13] To obtain a broadband QLg model, we measure the Lg spectral amplitudes at 58 frequencies distributed log evenly between 0.05 and 10.0 Hz. At each frequency f the corresponding spectral amplitude is calculated using a root

6 of 22

Figure 3. Sample records from earthquake 20040324. Shown here are normalized vertical ground velocities ordered according to their epicentral distances. Station names and maximum amplitudes as micrometers per second are listed at the left. Numbers on the waveforms indicate apparent group velocities. Inset: map showing the great circle paths from the epicenter to stations.

7 of 22

Figure 4. Simulated WorldWide Standardized Seismograph Network shortperiod vertical displacement seismograms for the same event shown in Figure 3. Maximum amplitudes are in micrometers.

8 of 22

not assume any a priori frequency dependence for QLg or the source function. The relationship between the Q value and the observed spectral amplitude is nonlinear. We use perturbation theory to linearize the relationship among the Q model, the source function, and the observed Lgwave spectrum. Finally, an iterative method is used to solve the problem. [17] In the Appendix, we obtain a linear system,

H A Q E U; 1

where H is a vector composed of residuals between the observed and the synthesized Lg spectra, dQ is a vector composed of the perturbations of the Q model, matrix A sets up the relationship between Q perturbations and the observed Lgwave spectra, dU is a vector composed of the perturbations of source terms, and matrix E sets up the relationship between the source perturbation and the observed Lgwave spectra. Detailed expression of equation (1) is provided in the Appendix. [18] To start the calculation we use a unit source function and a constant initial QLg model. On the basis of previous work in this region, for example, Xie et al. [2006] and Zhao et al. [2008], the initial Q model is prescribed using Q( f ) = Q0 f h, with Q0 = 420 and h = 0.15. At each step we solve for perturbations dQ and dU by minimizing the vector H. These perturbations are used to update the QLg model and source function. We iterate this process until satisfactory convergence is obtained. By inverting independently for all frequencies we obtain the frequencydependent QLg model and the source functions. 3.3. Resolution Test [19] The resolution and covariance matrixes are often used for assessing the resolving power of an inversion system [e.g., Crosson, 1976; Phillips et al., 2005; Phillips and Stead, 2008]. However, for a largescale inversion, the calculation is often timeconsuming. Alternatively, the checkerboard test method has been used in traveltime tomography [e.g., Zelt, 1998; Morgan et al., 2002; Pei et al., 2007], waveform tomography [Rao et al., 2006], and imaging of regional seismic propagation efficiency [e.g., Calvert et al., 2000; AlDamegh et al., 2004; Pei et al., 2006]. Xie [2006] discussed the resolving power of Q tomography and argued that the checkerboard resolution test used in velocity tomography yields a less satisfactory result in Q tomography. We adopt the checkerboard test method in our QLg tomography for resolution analysis. Bearing in mind that the signaltonoise ratios vary between frequencies, the available data points are different, resulting in different resolutions at different frequencies. Therefore, the resolution analyses have to be conducted for individual frequencies. [20] We first create a QLg model by superposing checkerboardshaped positive and negative perturbations on a Figure 5. Data processing procedure for Lg waves. (a) Original seismogram, (b) velocity record after deconvolving with the instrument response, (c, d) windowed Lg phase and preP noise, (e) Lgwave and noise spectra, (f) signalto noise ratio, and (g) Lgwave spectrum. Note that the data points have been dropped where the signaltonoise ratio is below the threshold of 2.0.

Table 3. Truncation Distances
Frequency (Hz) <2.0 2.02.5 2.55.0 56.5 >6.5 Distance (km) 1000

9 of 22

Figure 6. Resolution analysis at 1 Hz frequency. (a) Event and station geometry and ray paths for 1 Hz Lg data, where the crosses and triangles are epicenters and stations, respectively. (b) Initial checkerboard Q model with the alternating positive and negative perturbations. (c, d) Input and retrieved source functions. (e) Retrieved checkerboard Q perturbation. background Q. The Q perturbation has a constant percentage relative to the background Q. Next, we use equation (A1) to generate a synthetic Lg spectrum data set. The epicenters and station locations are adopted from the actual observation geometry. The source functions are calculated using equation (A2), with the seismic moment M0 and the corner frequency fc estimated from the observed magnitude using empirical relations. At each frequency, only events and stations that actually provide data above the signaltonoise ratios are used to generate the checkerboard test data. To simulate the noise in the real data a 5% root mean square fluctuation is added to the test data. This synthetic data set is input to the inversion system and the result is compared to the checkerboard model. [21] Figure 6 demonstrates the analysis for 1 Hz data. Shown in Figure 6a are 1340 ray paths from all available regional events (crosses) and stations (open triangles). There is denser ray coverage in Northeast China than in surrounding areas. Illustrated in Figure 6b is the checkerboard test model, which is formed by superimposing 7% Q perturbations on a constant background model of Q0 = 420. Figure 6c gives the 125 preassigned source functions at 1 Hz and Figure 6d shows the retrieved source functions after inversion. The two groups of numbers are highly correlated, indicating that the source functions are properly retrieved.

10 of 22

grid sizes vary from 0.5 0.5 to 2 2, with an increment of 0.1. Generally, a fine grid size provides a higher resolution, which is preferred for revealing the correlation between structures and Lgwave attenuation. In contrast, a coarse grid leads to a reliable result but often lacks the details required for characterizing regional geology. Thus, a trade off must be accepted to balance inversion resolution with reliability. As an example, Figure 7 shows the retrieved 1 Hz checkerboard models with grid sizes of 0.7 0.7, 1.5 1.5, and for comparison. Clearly, the 0.7 0.7 grid causes many regions, even in the central area, to lack resolution. The 1.5 1.5 and grids provide excellent resolving power for almost the entire region except some border areas. However, these coarse grids do not provide enough resolution to reveal attenuation distributions in some basinmountain transition areas. The grid, shown in Figure 6e, appears to have the optimal tradeoff between the resolution and the reliability for the 1 Hz Lg spectrum data. We perform this resolution analyses for all 58 individual frequencies. Figure 8 summarizes the numbers of available rays versus frequency, with the shading illustrating the estimated resolution for particular frequencies obtained by the resolution analyses.

4. Results

[23] Our inversion estimates the QLg distributions in Northeast China and its vicinity at 58 discrete frequencies, along with the source spectra for 125 selected events. 4.1. QLg Distributions [24] Figures 9a, 9c, and 9e illustrate examples of QLg distributions at 0.5, 1.0, and 2.0 Hz, along with regional geological structures including major fault systems, sutures, and basins. Shown in Figures 9b, 9d, and 9f are ray coverage and resolution tests at these frequencies. The most prominent feature in the QLg models is that the highfrequency QLg is generally higher than the lower frequency values. Within the investigated region the QLg shows a general trend to increase from south to north and the basins are often characterized by low Q values compared to the relatively high Q values in volcanic mountain areas. The active faults delimiting the basin and mountain areas are often related to strong gradients in Q. For example, the TanLu Fault is located on a belt of
Figure 7. Comparisons of 1 Hz spatial resolutions for different grid sizes. Grid sizes of (a) 0.7 0.7, (b) 1.5 1.5, and (c) 2 2. Figure 6e illustrates the reconstructed Q model. Comparing Figure 6e with 6b we see that the original checkerboard pattern is mostly resolved in Northeast China compared to surrounding areas and it is consistent with the ray coverage. [22] For an individual frequency we conduct resolution tests by using a series of checkerboard models. The model
Figure 8. Numbers of rays for frequencies between 0.05 and 10.0 Hz, along with results from resolution analyses.

11 of 22

Figure 9. Lg Q distributions at frequencies of 0.5, 1.0, and 2.0 Hz, along with their ray coverage and resolution analyses.

12 of 22

Figure 10. Frequencydependent Lg Q for selected subregions. Gray crosses are directly measured QLg. Filled circles with error bars are average values. Q0 values are labeled in each plot. Standard deviations on a logarithmic scale are converted to Q values and listed in parentheses. strong Q variation. To the east of this belt the CMs is a volcanic area that has relatively high Q values. To the west of the fault is the Songliao Basin (B4), which bears the largest hydrocarbon deposit in China. During development the basin has experienced a complex process including asthenosphere upwelling, rifting, postrift thermal subsidence, and structural inversion [Feng et al., 2010]. Relatively low Q values are observed in Songliao Basin (also refer to Figure 1). [25] To characterize the Lgwave attenuation for different geological formations we investigate regional variations and frequency dependence of QLg in different geological units. Figure 10 shows inverted broadband QLg versus frequency

13 of 22

Figure 11. Lg Q versus frequency for different geological units. Parameters for powerlaw Q models obtained by fitting the Lg Q between 0.5 and 1.5 Hz are listed at the right. in selected subregions. The light gray crosses are directly read from inverted Q(x,y, f ) maps (e.g., Figures 9a, 9c, and 9e) within the interested geological formations. We calculate the mean QLg values as functions of frequency in different regions. At frequency f the mean value is obtained by averaging all observations within the specific geological unit and between p p frequencies f/ 2 and f 2. The mean QLg values are shown in each plot by filled circles, along with their standard deviations. Note that the 1 Hz QLg in Figure 9c is directly inverted from the 1 Hz Lg data, while here the Q0 is an average within a frequency band and within a chosen geology unit. Thus the latter is more robust and representative. [26] Figure 10a shows the observed QLg for the entire investigated region, where the average Q0 is 414, but with a large scatter. The standard deviation is equivalent to Q0 values of from 232 to 739; this is listed in parentheses. Surrounded by largescale plates the investigated region is composed of three major units (refer to the caption to Figure 1): Siberian Craton (I), Northeast China Collage plate (II), and North China Craton (III). Figures 10b and 10c show that the average Q0 for the Northeast China Collage plate is 469 (306717), and that for the North China Craton is 338 (170670). The current investigation covers only part of the Siberian Craton, within which the average Q0 is 540 (321911) (not shown here). The average Q0 values reveal a significant difference in crustal attenuation between these firstorder geological units, while their relatively large standard deviations indicate strong variations within these units. [27] We focus our attention on the Northeast China Collage plate, where smaller scale structures include several main sedimentary basins (B1B6) and mountain areas (DMs and CMs). Figures 10d to 10g show some examples from basin and mountain areas. Figure 9 reveals that basins are generally characterized by low Q values, which may result from the strong attenuation in sedimentary stratification or because of the blockage of Lg waves at the edge of basins. The lowest mean Q0, 155 (72336), appears in Bohai Bay basin (B5). On the contrary, high Q0 values are found in mountain areas, with 675 (559814) for DMs and 630 (459 to 864) for CMs. The Northeast China Collage plate is composed of approximately a dozen massifs including folding belts and microplates. Shown in Figures 10h and 10i are two examples of the average QLg in these massifs. The standard deviations at this level of geology units are usually smaller than those for larger units but their frequency dependence appears to be more complex. 4.2. The Frequency Dependence of QLg [28] Figure 11 summarizes the average QLg versus frequency for different geological units in Northeast China and the Korean Peninsula (refer to Figure 1). In general, the mean QLg values increase gradually with increasing frequencies of between 0.05 and 1 Hz but rise more steeply above 1 Hz. The lowfrequency QLg appears to have larger regional variations and a more complex frequency dependence than at high frequencies. By comparing the Q distributions in Figure 9 with the crustal model (Laske, G., G. Masters, and C. Ref (2001), CRUST 2.0: A new global crustal model at 2 2; available at http://igppweb.ucsd.edu/gabi/crust2.html) in this region, there appears to be some correlation between QLg and the crustal thickness. Zhang and Lay [1995] investigated the effect of crust thickness on Lgwave propagation efficiency. They found the crust thickness to be closely related to the maximum number of overtone modes a waveguide can carry, and this number is a dominant factor controlling the Lgwave propagation efficiency. The number of overtones is proportional to fH [Zhang and Lay, 1995], where H is the crust thickness. Although those authors focused on highfrequency Lg waves propagating in continental and oceanic crusts, based on their theory the lowfrequency Lg waves should be more sensitive to the crustal thickness. The actual Lgwave attenuation should be dependent on both the geometrical parameters of the crust waveguide and its physical properties.

15 of 22

region. Compared to the lowfrequency results the h values in basins are less scattered but still distributed between 0.35 and 1.1. The values for massifs of the Northeast China Collage plate and mountain area are distributed from 0.35 to 0.7 and from 0.2 to 0.45, respectively. For the high frequencies (1.04.0 Hz) in Figure 13c the distribution of h appears to be more concentrated and tends to decrease with increasing QLg. In general, with increasing frequency the h values are shifted upward, a phenomenon shown in Figure 11, where the QLgfrequency curves increase in steepness at high frequencies. Our results demonstrate that for a large region, for example, the entire investigated region or the three subregions, the QLg h values at 1 Hz are between 0.5 and 0.8 and are relatively stable. However, for a small area and within a limited frequency band, the h values may vary to reflect the complex frequency dependence. The aforementioned behavior of h and QLg values and regional geology may indicate that the seemingly unstable observations of h values by previous authors [e.g., Xie et al., 2006; Chung et al., 2007] were not simply caused by insufficient data but actually revealed the complex frequency dependence of QLg. Wu and Aki [1985] reported observations of h and investigated their relation to Q0 based on the fractal nature of heterogeneities in random models. Campillo [1990] summarized several highfrequency QLg observations and found that the tectonic regions are often associated with low Q values and strong frequency dependence, while stable areas have high Q values and weak frequency dependence. These results are consistent with ours, particularly at high frequencies. 4.3. LgWave Source Spectral Functions [32] The simultaneous inversion also provides us with the Lgwave excitation functions at discrete frequencies. Figure 14 illustrates source functions for selected events. The crosses represent results for individual frequencies. To obtain the scalar seismic moment M0 and the corner frequency fc, we fit the Lgwave excitation spectrum with the w2 source model [Brune, 1970; Street et al., 1975; Sereno et al., 1988]. With this model, equation (A2) in the Appendix becomes
S f M 0 ; 2 4vf 2 =fC S 2
Figure 13. Comparison of distributions of h values (local slopes on QLgfrequency curves) for different geological units and for different frequencies: (a) 0.21.0 Hz, (b) 0.5 1.5 Hz, and (c) 1.04.0 Hz. behavior shown in Figure 11, where the QLg curves from basins often have turning points or even undulations. In contrast, for massifs of the Northeast China Collage plate (excluding those overlapping with basins), h values are mostly distributed between 0.4 and 0.6. For mountain areas h values are between 0.4 and 0.75, although there are only two data points. Figure 13b shows the results for middle frequencies, between 0.5 and 1.5 Hz, which are centered at 1 Hz, and the statistics should be comparable to the conventional narrow band powerlaw Q model using a 1 Hz reference frequency. The nominal h is close to 0.55 for the entire investigated

References

AlDamegh, K., E. Sandvol, A. AlLazki, and M. Barazangi (2004), Regional seismic wave propagation (Lg and Sn) and Pn attenuation in the Arabian Plate and surrounding regions, Geophys. J. Int., 157, 775795. Bouchon, M. (1982), The complete synthesis of seismic crustal phases at regional distances, J. Geophys. Res., 87, 1,7351,741. Brune, J. N. (1970), Tectonic stress and the spectra of seismic shear waves from earthquakes, J. Geophys. Res., 75, 49975009. Cai, X. L., J. S. Zhu, J. M. Cao, Z. Q. Yan, Z. X. Yang, and X. H. Hong (2002), Structure and dynamics of lithosphere and asthenosphere in the gigantic East AsianWest Pacific rift system (in Chinese with English abstract), Geol. China, 29, 234245. Calvert, A., E. Sandvol, D. Seber, M. Barazangi, F. Vidal, G. Alguacil, and N. Jabour (2000), Propagation of regional seismic phases (Lg and Sn) and Pn velocity structure along the AfricaIberia plate boundary zone: tectonic implications, Geophys. J. Int., 142, 384408. Campillo, M. (1987), Lg wave propagation in a laterally varying crust and the distribution of the apparent quality factor in central France, J. Geophys. Res., 92, 1260412614. Campillo, M. (1990), Propagation and attenuation characteristics of the crustal phase Lg, Pure Appl. Geophys., 132, 120. Campillo, M., B. M. Feignier, M. Bouchon, and N. Bthoux (1993), Attenuation of crustal waves across the Alpine range, J. Geophys. Res., 98, 19871996. Chi, G. S. (1988), The study of Cenozoic basalts and upper mantle beneath Eastern China (in Chinese with English abstract), pp. 1216, China University of Geosciences Press, Wuhan. Chung, T. W., M.H. Noh, J.K. Kim, Y.K. Park, H.J. Yoo, and J. M. Lees (2007), A study of the regional variation of lowfrequency Q1 around the Lg Korean Peninsula, Bull. Seismol. Soc. Am., 97, 21902197, doi:10.1785/ 0120070066.
Cong, L., and B. Mitchell (1998), Lg coda Q and its relation to the geology and tectonics of the Middle East, Pure Appl. Geophys., 153, 563585. Cong, L., J. Xie, and B. J. Mitchell (1996), Excitation and propagation of Lg from earthquakes in central Asia with implications for explosion/ earthquake discrimination, J. Geophys. Res., 101, 2777927789. Crosson, R. S. (1976), Crustal modelling of earthquake data, 1. Simultaneous least squares estimation of hypocenter and velocity parameters, J. Geophys. Res., 81, 30363046. Drouet, S., A. Souriau, and F. Cotton (2005), Attenuation, seismic moment, and site effects for weakmotion events: Application to the Pyrenees, Bull. Seismol. Soc. Am., 95, 17311748, doi:10.1785/0120040105. Efron, B. (1983), Estimating the error rate of a prediction rule: Improvement on crossvalidation, J. Am. Stat. Assoc., 78, 316331. Fan, G. W., and T. Lay (2002), Characteristics of Lg attenuation in the Tibetan Plateau, J. Geophys. Res., 107(B10), 2256, doi:10.1029/ 2001JB000804. Fan, G. W., and T. Lay (2003a), Strong Lg wave attenuation in the northern and eastern Tibetan Plateau measured by a twostation/twoevent stacking method, Geophys. Res. Lett., 30(10), 1530, doi:10.1029/ 2002GL016211. Fan, G. W., and T. Lay (2003b), Strong Lg attenuation in Tibetan Plateau, Bull. Seismol. Soc. Am., 93, 22642272. Feng, Z., C. Jia, X. Xie, S. Zhang, Z. Feng and T. A. Cross (2010), Tectonostratigraphic units and stratigraphic sequences of the nonmarine Songliao basin, northeast China, Basin Res., 22, 7995, doi:10.1111/ j.1365-2117.2009.00445.x. Fisk, M. D. (2006), Source spectral modeling of regional P/S discriminants at nuclear test sites in China and the former Soviet Union, Bull. Seismol. Soc. Am., 96, 23482367, doi:10.1785/0120060023. Ford, S. R., D. S. Dreger, K. Mayeda, W. R. Walter, L. Malagnini, and W. S. Phillips (2008), Regional Aattenuation in Northern California: A comparison of five 1D Q methods, Bull. Seismol. Soc. Am., 98, 20332046, doi:10.1785/0120070218. Ford, S. R., W. S. Phillips, W. R. Walter, M. E. Pasyanos, K. Mayeda, and D. S. Dreger (2009), Attenuation tomography of the Yellow Sea/Korean Peninsula from codasource normalized and direct Lg amplitudes, Pure Appl. Geophys., doi:10.1007/s00024-009-0023-2. Gao, X., W. M. Wang, and Z. X. Yao (2005), Crustal structure of China mainland and its adjacent regions (in Chinese), Chin. J. Geophys., 48, 591601. Ge, H., Y. Ni, C. Huang, J. Li, Y. Li, Y. Jiang, and F. Lu (1989), Lg velocity, attenuation and magnitude of the continent area in China (in Chinese), Chin. Sci. Bull., 24, 18891892. Ge, X. H., and W. P. Ma (2007), MesozoicCenozoic tectonic framework of southern Northeast Asia (in Chinese with English abstract), Geol. China, 34, 212228. Herrmann, R. B., and A. Kijko (1983), Modeling some empirical Lg relations, Bull. Seismol. Soc. Am., 73, 157172. Huang, C., Z. Gao, L. Shi, and Y. Ni (1990), Lg wave attenuation and magnitude in Northeast China region (in Chinese), Northeast. Seismol. Res., 6, 4149. Jih, R. S. (1998), Location calibration effects in China, in Proceedings, 20th Annual Seismic Research Symposium on Monitoring a Comprehensive Test Ban Treaty, 2123 September 1998, pp. 4445. Jin, A., and K. Aki (1988), Spatial and temporal correlation between coda Q and seismicity in China, Bull. Seismol. Soc. Am., 78, 741769. KadinskyCade, K., M. Barazangi, J. Oliver, and B. Isacks (1981), Lateral variation in high frequency seismic wave propagation at regional distances across the TurkishIranian plateaus, J Geophys. Res, 86, 93779396. Kennett, B. L. N. (1984), Guided wave propagation in laterally varying media; I. Theoretical development, Geophys. J., 79, 235255. Kennett, B. L. N. (1986), Lg waves and structural boundaries, Bull. Seismol. Soc. Am., 76, 11331141. Knopoff, L., F. Schwab, and E. Kausel (1973), Interpretation of Lg, Geophys. J. R. Astron. Soc., 33, 389404. Li, X., and X. Yuan (2003), Receiver functions in northeast China Implications for slab penetration into the lower mantle in Northwest Pacific subduction zone, Earth Planet. Sci. Lett., 216, 679 691. Ma, W. P., S. N. Lu, and G. S. Wang (2003), Tectonic positioning of the islands of Japan in the Mesozoic Asian frontier and its relation to the regional geology of eastern China (in Chinese with English abstract), Geol. Bull. China, 22, 192199. Mayeda, K., and W. R. Walter (1996), Moment, energy, stress drop, and source spectra of western United States earthquakes from regional coda envelopes, J Geophys. Res, 101, 11,19511,208. McNamara, D. E., T. J. Owens, and W. R. Walter (1996), Propagation characteristics of Lg across the Tibetan Plateau, Bull. Seismol. Soc. Am., 86, 457469.

21 of 22

Mellors, R., F. Vernon, V. Camp, A. AlAmri, A. Ghalib, and M. AlDail (1999), Regional wave propagation in the Saudi Arabian peninsula, Geophys. Res. Lett., 104, 2022120235. Menke, W., R. C. Holmes, and J. Xie (2006), On the nonuniqueness of the coupled origin timevelocity tomography problem, Bull. Seismol. Soc. Am., 96, 11311139. Mitchell, B. J., Y. Pan, J. Xie, and L. Cong (1997), Lg coda Q variation across Eurasia and its relation to crustal evolution, J. Geophys. Res., 102, 22276722780. Mitchell, B. J., L. Cong, and G. Ekstrm (2008), A continentwide map of 1Hz Lg coda Q variation across Eurasia and its relation to lithospheric evolution, J. Geophys. Res., 113, B04303, doi:10.1029/2007JB005065. Morgan, J. S., G. L. Christeson, and C. A. Zelt (2002), Testing the resolution of a 3D velocity tomogram across the Chicxulub crater, Techonophysics, 355, 215226. Mueller, R. A., and J. R. Murphy (1971), Seismic characteristics of underground nuclear explosions, Part I, Seismic spectrum scaling, Bull. Seismol. Soc. Am., 61, 16751692. Nuttli, O. W. (1973), Seismic wave attenuation and magnitude relations for eastern North America, J. Geophys. Res., 78, 876885. Nuttli, O. W. (1983), Average seismic sourceparameter relations for midplate earthquakes, Bull. Seismol. Soc. Am., 73, 519535. Nuttli, O. W. (1986), Yield estimates of Nevada Test Site explosions obtained from seismic Lg waves, J. Geophys. Res., 91, 21372151. Ottemller, L. (2002), Lg wave Q tomography in Central America, Geophys. J. Int., 150, 295302. Ottemller, L., N. M. Shapiro, S. K. Singh, and J. F. Pacheco (2002), Lateral variation of Lg wave propagation in southern Mexico, J. Geophys. Res., 107(B1), 2008, doi:10.1029/2001JB000206. Paige, C. C., and M. A. Saunders (1982), LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Softw., 8, 4371. Patton, H. J., and J. Schlittenhardt (2005), A transportable mb (Lg) scale for central Europe and implications for lowmagnitude Msmb discrimination, Geophys. J. Int., 163, 126140. Pei, S., J. Zhao, C. A. Rowe, S. Wang, T. M. Hearn, Z. Xu, H. Liu, and Y. Sun (2006), ML amplitude tomography in North China, Bull. Seismol. Soc. Am., 96, 15601566, doi:10.1785/0120060021. Pei, S., J. Zhao, Y. Sun, Z. Xu, S. Wang, H. Liu, C. A. Rowe, M. N. Toksz, and X. Gao (2007), Upper mantle seismic velocities and anisotropy in China determined through Pn and Sn tomography, J. Geophys. Res., 112, B05312, doi:10.1029/2006JB004409. Phillips, W. S., and R. J. Stead (2008), Attenuation of Lg in the western US using the USArray, Geophys. Res. Lett., 35, L07307, doi:10.1029/ 2007GL032926. Phillips, W. S., H. E. Hartse, and S. R. Taylor (2000), 1 Hz Lg Q tomography in central Asia, Geophys. Res. Lett., 20, 34253428. Phillips, W. S., H. E. Hartse, and J. T. Rutledge (2005), Amplitude ratio tomography for regional phases Q, Geophys. Res. Lett., 32, L21301, doi:10.1029/2005GL023870. Rao, Y., Y. Wang, and J. V. Morgan (2006), Crosshole seismic waveform tomography: II. Resolution analysis, Geophys. J. Int., 166, 12371248. Rodgers, A., J. Ni, and T. Hearn (1997), Propagation characteristics of short period Sn and lg in the Middle East, Bull. Seismol. Soc. Am., 87, 396413. Sandvol, E., K. AlDamegh, A. Calvert, D. Seber, M. Barazangi, R. Mohamad, R. Gok, N. Turkelli, and C. Gurbuz (2001), Tomographic imaging of observed regional wave propagation in the Middle East, Pure Appl. Geophys., 158, 11211163. Sereno, T. J., Jr., S. R. Bratt, and T. C. Bache (1988), Simultaneous inversion of regional wave spectra for attenuation and seismic moment in Scandinavia, J. Geophys. Res., 93(B3), 20192035. Shapiro, N., N. Bthous, M. Campillo, and A. Paul (1996), Regional seismic phases across the Ligurian sea: Lg blockage and oceanic propagation, Phys. Earth Planet Inter., 93, 257268.