# Seismological Research Letters

- © 2014 by the Seismological Society of America

*Online Material: *Code for computing *V*_{S} by the r/z ratio method; seismogram plots.

## INTRODUCTION

In contrast to the western United States where seismicity is high and seismic networks are dense, the central and eastern United States (CEUS) is covered with sparse seismic networks. The relatively low seismicity in the CEUS has also led to a limited dataset of waveform data from local earthquakes, hindering the development of reliable ground‐motion prediction equations in the CEUS. This situation has been changing with the Earthscope Transportable Array (USArray) program that started installing about 1000 broadband and short‐period stations in the CEUS starting in 2009, and will fully cover the CEUS soon. Such large numbers of stations provide a rich dataset for ground‐motion studies in the United States. Also, the increasing number of Advanced National Seismic System (ANSS) broadband stations with high‐sampling rate substantially augments the waveform data.

Site responses at each seismic station need to be obtained before the waveform data can be used for ground‐motion studies. For example, Boore (2003) and Atkinson and Boore (2006) find that site amplification needs to be taken into account for ground‐motion modeling in eastern North America, especially for high frequencies. Site response is controlled by the subsurface shear‐velocity profile, and the velocity structure of the top ten to hundreds of meters is particularly important (Wald and Mori, 2000; Boore, 2006). Therefore an essential part of site‐response modeling is determining subsurface velocity structure.

A detailed review of methods for determining subsurface velocity structure was made by Boore (2006). Theoretically, invasive methods such as borehole logging provide the most accurate measurement of velocity structure, but their cost prohibits extensive application in site‐response studies. Noninvasive geophysical exploration methods are more often used in determining shallow velocity structure, and these include spectral analysis of surface waves (SASW), multichannel analysis of surface waves (MASW), and refraction microtremor (ReMi) (Louie 2001; Stephenson *et al.*, 2005; Odum *et al.*, 2010). These noninvasive methods typically employ arrays of seismic receivers in deriving either travel‐time curves of body waves or dispersion curves of surface waves excited by a man‐made source or by natural ambient noise. Because of the close spacing between seismic receivers, array‐based methods can provide high‐resolution velocity profiles at relatively low cost.

However, array‐based geophysical methods can only be conveniently applied to a relatively small number of sites because they are costly and time consuming. It would be very challenging to perform subsurface structure studies using array‐based geophysical exploration methods for all the 400+ stations in the USArray and ANSS. Instead, methods such as Horizontal‐to‐Vertical (H/V) ratio of ambient noise and Horizontal‐to‐Vertical ratio of *S*‐wave spectral Ratio (HVSR) only require waveforms at single stations from recordings of ambient noise or local *S* waves (Nakamura, 1989; Lermo and Chávez‐García, 1993). Because all the USArray and ANSS stations provide continuous waveforms, these single station‐based methods can be readily applied to infer site responses.

Both the H/V and HVSR methods use the spectral ratio between horizontal and vertical components as a proxy of site response. The difference between H/V and HVSR is that the H/V method is applied to the waveform of ambient noise whereas HVSR is applied to *S* waves from local earthquakes. Though it is associated with some controversy (Lachet and Bard, 1994; Al Yuncha and Luzon, 2000; Langston *et al.*, 2009), the H/V method is widely applied because seismic stations record ambient noise most of the time (Lozano *et al.*, 2008). The H/V method assumes that the ambient noise consists mostly of Rayleigh waves, and therefore the ellipticity of particle motion can be used to infer subsurface response (Nakamura, 1989, 2000). Some studies directly use the spectral ratio of H/V as a proxy for site amplification (Read *et al.*, 2008). There are also efforts to invert the subsurface velocity profile from the H/V ratio (Arai and Tokimatsu, 2004). The HVSR method is also widely adopted in modeling site response, but requires *S* waveforms from local earthquakes (Lermo and Chávez‐García, 1993). For example, earthquakes in the magnitude range of 2.0–3.8 have been used in modeling site responses for POLARIS stations (Read *et al.*, 2008).

In this paper, we introduce a new method of characterizing site effects by modeling the H/V ratio of initial portion of local *P* waves. We first present the theoretical basis for the method, then demonstrate its feasibility with both synthetic and observed seismograms, and last, validate this method by comparing subsurface shear‐velocity models inverted from initial *P* waves with those measured from array methods. We also provide a package of computer scripts and codes for measuring subsurface shear velocity as supplementary material. Detailed validation of the method with other approaches and its application to the estimation of site characterization for ANSS seismic station CBN, recording the 2011 Virginia earthquake, are presented in Kim *et al.* (2013) and Li *et al.* (2013).

## THE THEORY OF ESTIMATING SHALLOW SHEAR VELOCITY FROM LOCAL *P* WAVEFORMS

As for *S* waves, *P* waves should be able to provide constraints on subsurface structure. For the case of a layer of unconsolidated sediments over bedrock, various seismic phases are observed due to conversion of waves at the interface (*SP*, incident *S* wave converted to *P* wave; or *PS*, incident *P* wave converted to *S* wave), and are used to improve estimates of the velocity structure in the sediments (Chen *et al.*, 1996; Langston, 2003). The same approach can be used to model near‐surface structure (tens to hundreds of meters), when the first tenths of a second of the *P* waveforms are examined. As most earthquakes occur in crystalline basement where seismic velocities are high, and most stations are situated on material of lower velocities, the incidence angles of *P* and *PS* near the seismic station are typically much smaller than the takeoff angle at the earthquake. Thus intuitively, the *P* wave is weaker on the radial component than on the vertical component because of the steeper ray path of the *P* wave, and *PS* is stronger on the radial component because its polarization is perpendicular to the ray path and almost horizontal.

Because of the interaction of incoming *P* or *S* waves with the free surface, the actual particle motion recorded by a seismic station is different from that of free body waves. From Aki and Richards (2002), the particle motions of incoming *P* and *S* waves are described by equations (1) and (2), respectively, in which *U*_{R}, *U*_{T}, *U*_{Z} are radial, tangential, and vertical components of the particle velocity at the free surface; *α* is the subsurface *P* velocity; *β* is the subsurface shear velocity; *p* is the ray parameter (horizontal slowness); *i* is the angle between the *P* rays and the vertical axis; and *j* is the angle between the *S* ray and the vertical axis (Fig. 1a): (1)(2)

Then from equation (1), it is straightforward to obtain the ratio between the radial and vertical components of *P* waves, (3)

And similarly for the *S* wave, but for *U*_{Z} over *U*_{R}, (4)

For the incoming *P* wave, the ratio of radial to vertical component of particle motion (as displayed in equation 3) depends mostly on the shear velocity *β* and ray parameter *p* because cos*j* is close to one due to the small shear velocity near the free surface. From equation (4), it follows that *PS* is much stronger on the radial component.

Because the ratio of the *P* wave is proportional to the subsurface shear velocity, this is then a good way to measure subsurface shear velocity when the ray parameter *p* is known. Fortunately, the crustal structure in the CEUS is relatively well studied (Ou and Herrmann, 1990), and *p* can be determined once the location and depth of an earthquake is determined.

## NUMERICAL SIMULATION OF THREE‐COMPONENT *P* WAVES FOR DIFFERENT SUBSURFACE STRUCTURES

In order to confirm that the ratio of the *P* wave is mostly sensitive to subsurface shear velocity, we computed synthetic seismograms of three‐component *P* waves. For 1D horizontally layered velocity models, two approaches can be taken for calculating synthetic waveforms. One approach is the wave number–frequency integral method that takes into account full waves for a point source with focal mechanism (Zhu and Rivera, 2002). This approach is time consuming when high frequencies (for example, 20 Hz and above) is involved. The other approach assumes that the initial portion of *P* waves is generated by incident plane *P* waves, thus an integral only in the frequency domain is required and computation is much faster (Randall, 1989). As displayed in Figure 2, synthetic seismograms from the two approaches agree well. Hereafter, we adopt the plane‐wave approximation for computing synthetic seismograms.

To demonstrate that the radial and vertical ratio of the *P* wave is sensitive to subsurface shear‐velocity structure, we compute synthetic seismograms for five velocity models, in which the top 30 m shear velocity (*V*_{S30}) corresponds to the five National Earthquake Hazards Reduction Program (NEHRP) site classes (A, B, C, D, E), respectively, (Fig. 3a). In Figure 3b, the radial component for each subsurface model is displayed, and the source duration is chosen to be very short (0.05 s) so that individual seismic phases are isolated from each other. It is obvious that for the radial component of synthetic seismograms, smaller subsurface shear velocity is always associated with smaller *P* amplitude and larger amplitude of secondary arrivals, which could consist of *PS* and its multiples. When the source time function is chosen to be 0.1 s (appropriate for **M** 3 earthquakes), a similar pattern can be observed in Figure 3c. In both Figure 3b and Figure 3d, the NEHRP site class A model shows the strongest *P* amplitude and the smallest secondary arrivals on the radial component. Although the subsurface shear‐velocity models affect the radial component of *P* waves, the vertical component (black traces) is not affected much (Fig. 3d). From this figure, it is clear that the ratio of the *P* wave is indeed very sensitive to subsurface velocity. For NEHRP sites C and E, the stronger secondary arrivals lag the direct *P* wave by 0.15 and 0.3 s, respectively. Such a lag can be readily observed on broadband or short period seismograms with sampling rate of 20 samples per second or higher.

We also test whether the radial component of the *P* wave depends on the subsurface *P*‐wave velocity. In Figure 4a, a 100 m layer with shear velocity of 300 m/s situated on half‐space and two *P*‐velocity models are tested (*V*_{P} of 0.75 and 1.6 km/s, respectively). Synthetic seismograms for the two *P*‐velocity models are almost identical, indicating that subsurface *P* velocity does not affect the radial component of initial *P* waves. Synthetic seismograms for the models with a top layer with linear velocity gradient (Fig. 4b) again confirms that the radial component of the initial *P* wave is not sensitive to subsurface *P* velocity.

As the ratio of radial and vertical initial *P* waves has to be measured over a finite time window, the inferred subsurface shear velocity is not necessarily the shear velocity for depth of zero. Instead, the inferred shear velocity should represent the average from surface to some depth. Usually the resolving capability of seismic‐body waves is limited by wavelength. Therefore, the inferred subsurface shear velocity could be the averaged velocity for the depth range of one wavelength of the *S* wave. We test the depth range that can be sampled with the vertical‐to‐horizontal ratio of *P* waves in Figure 5. The thickness of the top layer varies from 10 to 80 m with increments of 10 m. The earthquake source time function is set to be a 0.1 s wide pulse (typical for **M** 3 earthquakes). From the synthetic seismograms, it is observed that only when the thickness is larger than the wavelength of *S* waves (0.1 s times the shear velocity), does the radial component of initial *P* waves have the same waveform as the vertical component and show the expected amplitude. That is to say, the depth extent resolvable from the ratio of radial to vertical *P*‐wave amplitudes depends on the earthquake source duration.

## OBSERVED THREE‐COMPONENT *P* WAVES FROM LOCAL EARTHQUAKES AND VALIDATION WITH MEASURED SHEAR‐WAVE VELOCITY PROFILES

An observation from the sequence of the 18 April 2008 Mt. Carmel earthquake is displayed in Figure 6 for the ANSS station OLIL. For each aftershock, a clear *P* wave can be observed and the radial component of the *P* wave is much weaker than the vertical component for the first 0.2 s. Moreover, on the radial component there are strong secondary arrivals after the direct *P* wave, with features similar to those in Figure 3 for the case of NEHRP site class C. On the tangential component ( Figure S1 available in the electronic supplement to this article), *P* waves are much weaker, suggesting that 1D modeling is applicable for this site.

Odum *et al.* (2010) performed detailed studies of five ANSS stations in the CEUS, and found that *V*_{S30} is 480 m/s for station OLIL from the ReMi technique. We take the subsurface velocity model from their ReMi result and compute synthetic seismograms of *P* waves. The synthetic seismograms (Fig. 7) explain major features on the observed *P* waveforms as displayed in Figure 6, with smaller *P* amplitude on the radial component, *PS* lagging *P* by 0.1 s, and a fairly strong arrival about 0.5 s after the *P* wave. Hence, the initial few tenths of a second of three‐component *P* waveforms have enough resolving power to estimate subsurface shear‐velocity profiles.

The waveforms of the radial components (red) in Figure 6 show some variability, and this could be related to the different waveforms of the vertical components (green). In a way analogous to teleseismic receiver functions, we can compute a local receiver function by deconvolving the radial component from the vertical component. The teleseismic receiver function is defined as the radial component of *P* waves that is deconvolved from the vertical component either in the frequency or the time domain (Ligorria and Ammon, 1999). Receiver function is mostly sensitive to velocity structures beneath the seismic station, because earthquake source contamination is removed after the deconvolution processes. The deconvolution is performed with an iterative procedure in the time domain (Ligorria and Ammon, 1999), and the results are displayed in Figure 8a. The deconvolved waveforms are very similar for different events, showing that the initial portion of *P* waves is mostly sensitive to structure beneath the receivers, and not very sensitive to the location and depth of the events.

The approach of deconvolving the radial from the vertical component is necessary when waveform data of larger earthquakes are used to study subsurface structure, because the source‐function time could be much longer than the time window (typically less than 1 s) useful for subsurface studies. For example, the local receiver function of the *M*_{w} 5.2 2008 Mt. Carmel mainshock (red trace in Fig. 8b) is very similar to those of the aftershocks (black traces), demonstrating the feasibility of modeling subsurface structure with waveform data from earthquakes up to **M** 5. There are subtle differences between the red and black traces, which may be caused by source complexities or differences in location and depth.

As mentioned above, the vertical *P* waveform is not much modified by subsurface structure. Therefore, resolving subsurface structure from local waveforms is very similar to resolving crustal structure with teleseismic *P*‐receiver functions (Langston, 1989), where it is assumed that vertical *P*‐wave propagation is simple. The major difference between resolving subsurface velocity structure and resolving crustal structure is temporal precision. Teleseismic *P* waves are only capable of providing robust planar *P* waves with frequency below 2 Hz (due to mantle attenuation and scattering of 3D crustal heterogeneity). In contrast, subsurface velocity structure can be resolved at frequencies as high as 10 Hz (for **M** 3 events) or 20 Hz (for **M** 2 events). Also, 3D crustal heterogeneity is less of a problem because only tens of meters to hundreds of meters in the region near the site are involved.

Consequently, algorithms used for resolving subsurface velocity structure with the initial portion of the local *P* wave can be the same as those used in teleseismic *P*‐receiver functions. Among many inversion schemes, global optimization methods such as the Neighborhood Algorithm and the Differential Evolution (DE) method have the advantage of finding the global minimum, and are particularly appropriate for highly nonlinear problems such as *P*‐receiver functions (Sambridge 1999a,b; Li *et al.*, 2010). We adopt the technique used by Li *et al.* (2010), who obtained a high‐quality crustal model at station BJT in China using teleseismic *P*‐receiver functions. To demonstrate that the DE inversion algorithm is effective for inverting shallow velocity structure, we apply the algorithm to local receiver functions recorded at station OLIL. In order to enhance signal‐to‐noise ratio and suppress effects due to 3D lateral variation, we stacked the local receiver function traces and inverted subsurface structure with the stacked trace. In the inversion, the subsurface velocity model is parameterized to have three layers. The thickness and seismic velocities (*V*_{P} and *V*_{S}) of each layer are chosen to be within wide ranges that are set *a priori*, and the optimal values of thickness and values are found through DE inversion algorithm. After 1000 iterations, a minimum in misfit between observed and predicted local receiver functions is achieved (Fig. 9). Overall, the shear‐velocity model for the top 60 m agrees with those measured from ReMi and Refraction/Reflection methods (Fig. 9d).

## SUMMARY

In summary, the initial portion of three component local *P* waves can be used to constrain subsurface shear‐velocity structure in a similar way to teleseismic *P*‐receiver functions, only at much higher temporal resolution and smaller spatial scale. More detailed structure can be achieved via modeling of local *P* and *S* waves involving whole waveform inversion or backprojection methods (Chiu and Langston, 2009, 2011). Instead, the algorithm presented in this study is straightforward. The shear velocity just beneath the free surface can be readily estimated from the ratio of the *P* wave after ray parameter *p* is determined (equation 3). Of course, the depth extent of subsurface shear velocity determined this way depends on the source duration of vertical *P*. For a local earthquake with duration of 0.1 s (**M** 3 earthquake), the shear velocity estimated from of *P* is actually the average velocity *β* from free surface to depth*h*=0.1 s×*β*. For NEHRP site class C that corresponds to a depth extent of 50 m. With additional constraints from noise H/V ratio, the subsurface shear‐velocity structure can be resolved better.

## ACKNOWLEDGMENTS

Seismic waveform data were obtained from Incorporated Research Institutions for Seismology Data Management Center (IRIS‐DMC). Supported Chinese Academy of Sciences (CAS) program KZCX2‐EW‐121, USGS National Earthquake Hazards Reduction Program (NEHRP), and Chinese Academy of Sciences (CAS) 40674027, 41204044.