Estimating near-surface shear-wave velocities in Japan by applying ...
CWP-732
Estimating near-surface shear-wave velocities in Japan by applying seismic interferometry to KiK-net data Nori Nakata & Roel Snieder Center for Wave Phenomena, Geophysics Department, Colorado School of Mines, Golden, Colorado 80401
ABSTRACT
We estimate shear-wave velocities in the shallow subsurface throughout Japan by applying seismic interferometry to the data recorded with KiK-net, a strongmotion network in Japan. Each KiK-net station has two receivers; one receiver on the surface and the other in a borehole. By using seismic interferometry, we extract the shear wave that propagates between these two receivers. Applying this method to earthquake-recorded data at all KiK-net stations from 2000 to 2010 and measuring the arrival time of these shear waves, we analyze monthly and annual averages of the near-surface shear-wave velocity all over Japan. Shear-wave velocities estimated by seismic interferometry agree well with the velocities obtained from logging data. The estimated shear-wave velocities of each year are stable. For the Niigata region, we observe a velocity reduction caused by major earthquakes. For stations on soft rock, the measured shear-wave velocity varies with the seasons, and we show negative correlation between the shear-wave velocities and precipitation. We also analyze shear-wave splitting by rotating the horizontal components of the surface sensors and borehole sensors and measuring the dependence on the shear-wave polarization. This allows us to estimate the polarization with the fast shear-wave velocity throughout Japan. For the data recorded at the stations built on hard-rock sites, the fast shearwave polarization directions correlate with the direction of the plate motion. Key words: seismic interferometry, shear-wave splitting, shear-wave velocity, precipitation, near surface, KiK-net
1
INTRODUCTION
Seismic interferometry is a powerful tool to obtain the Green’s function that describes wave propagation between two receivers (e.g., Aki, 1957; Claerbout, 1968; Lobkis and Weaver, 2001; Roux and Fink, 2003; Schuster et al., 2004; Wapenaar et al., 2004; Bakulin and Calvert, 2006; Snieder et al., 2006; Wapenaar and Fokkema, 2006). Seismic interferometry is applied to ambient noise (e.g., Hohl and Mateeva, 2006; Draganov et al., 2007, 2009; Brenguier et al., 2008; Stehly et al., 2008; Lin et al., 2009), traffic noise (e.g., Nakata et al., 2011), production noise (e.g., Miyazawa et al., 2008; Vasconcelos and Snieder, 2008), earthquake data (e.g., Larose et al., 2006, Sens-Sch¨ onfelder and Wegler, 2006; Snieder and S ¸ afak, 2006; Ma et al., 2008; Ruigrok et al., 2010), and active sources (e.g., Bakulin and Calvert, 2004; Mehta et al., 2008). In Japan, large seismometer networks, such as Hinet, F-net, K-NET, and KiK-net (Okada et al., 2004),
are deployed. By using these networks for seismic interferometry, Tonegawa et al. (2009) extract the deep subsurface structure of the Philippine Sea slab. These data have also been used to observe time-lapse changes in small regions (Wegler and Sens-Sch¨ onfelder, 2007; Sawazaki et al., 2009; Wegler et al., 2009; Yamada et al., 2010). Each KiK-net station has two receivers, one on the ground surface and the other at the bottom of a borehole. One can estimate the body-wave velocity between two receivers by using seismic interferometry (Trampert et al., 1993; Snieder and S ¸ afak, 2006; Mehta et al., 2007b; Miyazawa et al., 2008). By applying seismic interferometry to KiK-net data, we analyze near-surface velocities throughout Japan. Because KiK-net has recorded strong-motion seismograms continuously since the end of 1990s, the data are available for time-lapse measurements. Measuring time-lapse changes of the shallow subsurface is important for civil engineering and for estimating the site response to earthquakes. Previous studies ex-
316
N. Nakata & R. Snieder
a) 48 oN
b)
o 39 N
Jul. 16, 2007 M6.8
o
44 N o 38 N
o
Oct. 23, 2004 M6.8
NIGH13
40 oN o 37 N
36 N
Jan. 4, 2001 M5.3
32 oN
o 36 N o 137 E o o 128oE 132oE 136oE 140oE 144 E 148 E
o 138 E
o 139 E
o 140 E
Figure 1. KiK-net stations (December, 2010). The black dots on the map represent the locations of the stations. The dark gray shows the area analyzed in section 6.1. The light gray illustrates the area where we apply the analysis for seasonal change (section 6.2). The rectangle area in panel (a) is magnified in panel (b). The large black circle indicates station NIGH13, which we use for examples of analysis in Figures 2, 4, 5, 6, 9, and 10. The black crosses depict the epicenters of three significant earthquakes that occurred in the vicinity.
tracted time-lapse changes caused by earthquakes (Li et al., 1998; Vidale and Li, 2003; Schaff and Beroza, 2004; Wegler and Sens-Sch¨ onfelder, 2007; Brenguier et al., 2008). Interferometry applied to a single KiKnet station has also been used to measure time-lapse change due to earthquakes (Sawazaki et al., 2009; Yamada et al., 2010). Interferometric studies have shown changes in the shear-wave velocity caused by precipitation (Sens-Sch¨ onfelder and Wegler, 2006) and have measured shear-wave splitting (Bakulin and Mateeva, 2008; Miyazawa et al., 2008). We study the annual and monthly averages of the shear-wave velocity and the fast shear-wave polarization directions for stations all over Japan, and the temporal change in shear-wave velocity in the Niigata prefecture for three major earthquakes. This paper presents data processing of KiK-net data with seismic interferometry. We first introduce the properties of KiK-net. Next, we show the data analysis method. Then we present near-surface shear-wave velocities in every part of Japan. Finally, we interpret these velocities to study time-lapse changes, which are related to major earthquakes and precipitation, and present measurements of shear-wave splitting.
2
KIK-NET DATA
About 700 KiK-net stations are distributed in Japan (Figure 1). The stations are operated by the National
Research Institute for Earth Science and Disaster Prevention (NIED). Each station has a borehole and two seismographs which record strong motion at the bottom and top of the borehole. Each seismograph has three components: one vertical component and two horizontal components. Although the two horizontal components of the surface seismograph are oriented in the northsouth and east-west directions, respectively, the horizontal components of the borehole seismograph are not always aligned with the north-south and east-west directions because of technical limitations. Therefore, we rotate the directions of the borehole seismograph northsouth and east-west directions before data processing. The depth of about 25% of the boreholes is 100 m, and the other boreholes are at greater depth. Since our target is the near surface, we use the stations with a depth less than 525 m, which accounts for 94% of the stations. The sampling interval is either 0.005 or 0.01 s, depending on the station and the recording date. We show example records of an earthquake in Figure 2. Figure 2a illustrates bandpass-filtered time series, and Figure 2b the power spectra of the unfiltered data. As shown in Figure 2b, most energy is confined to 1-13 Hz, and we apply a bandpass filter over this frequency range for all data processing. In Figure 2, UD denotes the vertical component, NS the north-south direction horizontal component, and EW the east-west direction horizontal component. In Figure 2a, the P wave arrives
Near-surface S-wave velocities in Japan a)
surface
UD2
borehole
UD1
borehole
NS1
3
RETRIEVAL OF THE WAVEFIELD BETWEEN RECEIVERS
surface
EW 2
borehole
EW 1 10
20
30
40
50
Time (s)
b) UD2
surface
UD1
borehole
NS2
surface
NS1
borehole
EW 2
surface
EW 1
borehole
0
shear-wave velocity between these sensors as determined in this study. A bias in the velocity estimation due to non-vertical propagation depends on the deviation of cos θ from its value for vertical incidence, cos 0◦ = 1.
surface
NS2
0
317
10
20
30
40
50
Frequency (Hz) Figure 2. An earthquake recorded at all channels of station NIGH13 (latitude 37.0514◦ N and longitude 138.3997◦ E). This earthquake occurred at 14:59:19.56, 27 October 2004. The epicenter is at latitude 37.2204◦ N and longitude 138.5608◦ E and the depth is 11.13 km. The magnitude of this earthquake is MJM A 4.2. UD represents the vertical component, NS the north-south direction horizontal component, EW the east-west direction horizontal component, 1 the borehole seismograph, and 2 the surface seismograph. (a) The bandpass-filtered (1-13 Hz) time series. (b) The power spectra of the unfiltered records.
at around 7 s, and the shear wave arrives at around 14 s. All the used events are at a depth greater than 10 km. Because of this large depth compared to the depth of the boreholes and the low velocity in the near surface, the waves that travel between the sensors at each station propagate in the near-vertical direction as plane waves. We compute the angle of the incoming wave at the borehole receiver by using one-dimensional ray tracing to confirm that the wave propagating between the borehole and surface seismometers, propagates in the near-vertical direction. We use the velocity model of Nakajima et al. (2001) to determine the ray parameter p of the ray between each earthquake and the borehole sensor. The angle of incidence θ of the wave propagating from the borehole p receiver to the surface receiver is given by cos θ = 1 − p2 v 2 , where v is the average
We apply seismic interferometry to the recorded earthquake data of each station for retrieving the wavefield where the borehole receiver behaves as a virtual source. Several algorithms have been used in seismic interferometry to obtain the wavefield. These include cross correlation (e.g., Claerbout, 1968; Bakulin and Calvert, 2004), deconvolution (e.g., Trampert et al.), 1993; Snieder and S ¸ afak, 2006), cross coherence (e.g., Aki, 1957; Prieto et al., 2009), and multidimensional deconvolution (e.g., Wapenaar et al., 2008; Minato et al., 2011). We introduce the cross-correlation and deconvolution algorithms. We denote the wavefield, excited at source location s that strikes the borehole receiver at location rb by u(rb , s, ω) = S(rb , s, ω), where S(rb , s, ω) is the incoming wavefield that includes the source signature of the earthquake and the effect of propagation such as attenuation and scattering, in the frequency domain. The corresponding wavefield recorded at the surface receiver at location rs is given by u(rs , s, ω) = 2G(rs , rb , ω)S(rb , s, ω),
(1)
where the factor 2 is due to the presence of the free surface at rs . Because the wavefield striking the borehole receiver is close to a vertically propagating plane wave, G(rs , rb , ω) is the plane wave Green’s function that accounts for the propagation from the borehole seismometer to the surface seismometer. The cross-correlation approach to retrieve the wavefield in one dimension is given by Wapenaar et al. (2010) |S(rb , s, ω)|2 G(rs , rb , ω) =
2jω u(rs , s, ω)u∗ (rb , s, ω), ρc (2)
where ρ is the mass density of the medium, c the wave propagation velocity, j the imaginary unit, and ∗ the complex conjugate. The regularized deconvolution, which is similar to cross correlation, is given by
G(rs , rb , ω) =
u(rs , s, ω) u(rs , s, ω)u∗ (rb , s, ω) ≈ , u(rb , s, ω) |u(rb , s, ω)|2 + (3)
where is a regularization parameter (Mehta et al., 2007b,a). The deconvolution is potentially unstable due to the spectral devision, and we avoid divergence by adding a positive constant to the denominator (equation 3). Note that the deconvolution eliminates the
318
N. Nakata & R. Snieder
Time (s)
0.14
0.16
0.18
Figure 3. Quadratic interpolation. Using the three largest amplitude points (crosses), we interpolate the highest amplitude point (circle) by estimating the quadratic curve through the three highest amplitude points.
empirically that this is the smallest regularization parameter to obtain stable wavefields. We apply a bandpass filter from 1-13 Hz after applying deconvolution interferometry. In this paper, we average in three ways to interpret the wavefields. The first method is the annual stack, where we average the deconvolved waveforms over the earthquakes recorded in each year. In the second averaging method, we average the deconvolved waves over all earthquakes recorded in each month over the 11 years (January 2000–2010, February 2000–2010, ..., and December 2000–2010). We call this average the monthly stack. In the third method, which we use for analyzing the influence of major earthquakes, we average over three months after a significant earthquake.
4.1 imprint of waveform S(rb , s, ω), which is incident on the borehole receiver. We derive the features of crosscorrelation and deconvolution interferometry in Appendix A.
4
DATA PROCESSING
We use 111,934 earthquake-station pairs that are recorded between 2000 and 2010. The magnitude range is confined between 1.9 and 8.2. The cosine of the angle of incidence cos θ of the wave propagating between the receivers at each station is greater than 0.975, even for the events that are the furthest away. The bias introduced by non-vertical propagation thus is less than 2.5 %, and for most measurements it is much smaller. First, we check the data quality and drop some seismograms by a visual inspection using the signal-to-noise ratio as a criterion. Additionally, we discard stations with a borehole seismometer at a depth greater than 525 m because we focus this study on the near surface. We remove the DC component of the data by subtracting the average of each seismogram. For aligning the directions of the borehole receiver to the exact north-south and east-west directions, we rotate the borehole receiver using the rotation angle provided by NIED (Shiomi et al., 2003). Because the sampling interval is not small compared to the travel time of P waves between the borehole and surface seismometers, we focus on the shear wave and only analyze the horizontal components. We first apply deconvolution interferometry to the motion in the north-south direction of each surfaceborehole pair. In this study, for reasons explained below, deconvolution interferometry gives more consistent estimates of the Green’s function than does crosscorrelation interferometry. We choose in equation 3 to be 1% of the average power spectrum of the borehole receiver in the frequency range 1-13 Hz because we find
Estimating the shear-wave velocity
Before we apply annual stacking or monthly stacking, we resample the data from 0.005-s interval to 0.01-s interval if the data that are stacked include both 0.005-s and 0.01-s sampling-interval data. After stacking, we estimate the arrival time by seeking the three adjacent samples with the largest values and apply quadratic interpolation to find the time at which the deconvolved data have a maximum amplitude (Figure 3). This time is the travel time for a shear wave that propagates between the borehole and surface sensors. We use this travel time to compute the shear-wave velocity of the region between the two receivers.
4.2
Computing the average and standard deviation of the velocity of the annual or monthly stacks
To interpret time-lapse variations in the velocities, we need to compute the average and standard deviation of the velocities within a region. Let us denote the estimated velocity by vi (m, y), where vi is the shear-wave velocity at station i, in month m, and year y. This velocity is already averaged over each month. Each station has a different velocity. In order to quantify the timelapse variations of the velocity, we subtract the average value of each station before calculating temporal variation in the annual or monthly average: ∆vi (m, y) = vi (m, y) − v i ,
(4)
where v i is an average velocity of station i over all months and years. Then we compute either the annual or monthly average of the velocity variation over stations ∆v, and we also compute the standard deviation of this quantity.
Near-surface S-wave velocities in Japan
0.1
0.1
0.2
0.2
Tim e (s)
b) 0.0
Tim e (s)
a) 0.0
0.3
0.3
0.4
0.4
0.5
0.5
0.6
2000
0.6
2010
319
2000
2010
Figure 4. Annual-stacked wavefields by using (a) cross correlation and (b) deconvolution interferometry at station NIGH13. The surface and borehole receiver orientation directions are north-south. Epicenter locations are illustrated in Figure 5.
b) 48 oN
48 oN
100
75
44 oN
44 oN
2005: 19 40 oN
40 oN
2000: 5
2002: 9
36 oN
2004: 110
32 oN
o o 128 oE 132oE 136oE 140oE 144 E 148 E
25
2008: 10 2009: 14
2003: 5 32 oN
2006: 5 2007: 32
2001: 12 36 oN
50
Depth (km )
a)
2010: 10
0
3 5 7
o Magnitude o 128 oE 132oE 136oE 140oE 144 E 148 E
Figure 5. Epicenters used in Figure 4 from 2000 to 2004 (a) and from 2005 to 2010 (b) at station NIGH13. At the right in each panel, we show the number of earthquakes we use to obtain the waveforms in Figure 4 in each year. The size of each circle refers the magnitude of each earthquake and the color indicates the depth. The white triangle illustrates the location of station NIGH13. Because of the proximity of events, many circles overlap.
4.3
Comparison between cross-correlation and deconvolution interferometry
We compare the cross-correlation and deconvolution approaches using the annual-stacked wavefields (Figure 4). We show the locations of the epicenters of the used earthquakes in Figure 5. The annual stacks of the waveforms obtained by cross correlation are shown in Figure 4a. These waveforms are not repeatable from year to year and often do not show a pronounced peak at the
arrival time of the shear wave at around t = 0.15 s. We attribute the variability in these waveforms to variations in the power spectrum |S(rb , s, ω)|2 of the waves incident at the borehole receiver (equation A1). In contrast, the annual stacks of the waveforms obtained by deconvolution shown in Figure 4b are highly repeatable and show a consistent peak at the arrival time of the shear wave. The consistency of these waveforms is due to the deconvolution that eliminates the imprint of the incident wave S(rb , s, ω) (equation A2). Consistent with
320
N. Nakata & R. Snieder
earlier studies (Trampert et al., 1993; Snieder and S ¸ afak, 2006), we use deconvolution to extract the waves that propagate between the seismometers at each KiK-net station.
0.0
4.4
0.2
v(φ) = v0 + v1 cos 2φ + v2 sin 2φ.
(5)
In this expression, p v0 is the isotropic component of the velocity, and v12 + v22 the anisotropy. We assume the splitting time to be much smaller than the period of the wavefield. Because the wavefields of each polarization data are symmetric by 180 degrees, the anisotropy depends on polarization through a dependence of 2φ.
RETRIEVED NEAR-SURFACE SHEAR-WAVE VELOCITIES IN JAPAN
Using deconvolution interferometry at each station, we obtain the wavefield that corresponds to a plane wave propagating in the near-vertical direction (cos θ > 0.975) between the borehole receiver and surface receiver at each station. In this section, we show the wavefields of the annual stack, monthly stack, and shear-wave splitting.
Time (s)
Shear-wave splitting
We investigate shear-wave splitting by measuring the shear-wave velocity as a function of the polarization. We rotate both surface and borehole receivers from 0 to 350 degrees using a 10-degree interval. The north-south direction is denoted by 0 degrees, and the east-west direction by 90 degrees. Because a rotation over 180 degrees does not change the polarization, the 0- to 170-degree wavefields are the same as the 180- to 350-degree data. We apply deconvolution interferometry to the rotated wavefields, located at the surface and borehole receivers with the same polarizations, for determining the velocity of each polarization. Because the velocity for each polarization is related to the velocities of the fastest and slowest shear waves (Appendix B), we can estimate shear-wave splitting from the velocity difference. We cross-correlate the deconvolved wavefield for every used polarization (from 0 to 350 degrees in 10-degree intervals) with the deconvolved wavefield obtained from the motion in the north-south direction. This allows us to quantify the polarization dependence of the shear-wave velocity. Similar to the process described in section 4.1, we compute annual stacks of cross-correlated wavefields and pick the peak amplitudes of stacked wavefields by using quadratic interpolation. We can separate the velocity v(φ) as a function of polarization direction φ into the isotropic and anisotropic terms using a Fourier series expansion (Appendix C):
5
0.1
0.3 0.4 0.5 0.6
2000
2010
Figure 6. Annual-stacked wavefields (curves) with the interpolated largest amplitude (circles) at station NIGH13. The horizontal line at around 0.15 s is the shear-wave arrival time determined from logging data. From left to right, we show annual stacks from 2000 to 2010. The source and receiver polarization directions are the north-south direction.
5.1
Annual and monthly stacks
Figure 6 shows the annual-stacked wavefields at station NIGH13 represented by the large black circle in Figure 1. At this station, the sampling interval is 0.005 s until 2007 and is 0.01 s after 2008. In Figure 6, the deconvolved wavefields have good repeatability and a pronounced peak amplitude. After we apply quadratic interpolation (section 4.1), the determined arrival times (the black circles in Figure 6) correlate well with the travel time which is obtained from logging data (the horizontal line in Figure 6). The logging data is measured using a logging tool and by Vertical Seismic Profiling (VSP). The seismic source of VSP is a verticalcomponent vibrator. For finite offsets this source generates shear waves. We determine the average velocity from the logging data by computing the depth average of the slowness, because this quantity accounts for the vertical travel time. Because of the quadratic interpolation, the measured travel times show variations smaller than the original sampling time. After determining the arrival times of all stations, we compute the shear-wave velocities by using the known depth of the boreholes. Applying triangle-based cubic interpolation (Lawson, 1984) between stations, we create the shear-wave velocity map of Japan in each year (Figure 7). To reduce the uncertainty of the velocity estimation we use only the stations which give deconvolved waves with an arrival time greater than 0.1 s.
Near-surface S-wave velocities in Japan ◦
2 48 N 1 ◦
8E
◦E
2 13
◦E
6 13
◦E
0 14
◦E
◦
14
321
4 E 1 48
Log 44◦N
2000 2000
2005 36◦N
1500
2008 2010
1000
32◦N
Velocity (m /s)
40◦N
2002
500
0
Figure 7. Shear-wave velocities obtained from logging data (Log) and estimated by annual-stacked seismic interferometry using earthquake data at the north-south polarization (excerpted 2000, 2002, 2005, 2008, and 2010). The blue dots on a map represent the station locations which we use to make the map. We interpolate velocities between stations by triangle-based cubic interpolation (Lawson, 1984). The longitude and latitude belong to the map in the upper-left. The number of right-upper side of each map shows the year of data.
Thus, we obtain the near-surface shear-wave velocities throughout Japan by applying seismic interferometry to KiK-net data. The shear-wave velocity obtained from logging data is shown in the top left in Figure 7. Note that the velocities measured in different years are similar. In Figure 8, we crossplot the velocities estimated by interferometry in 2008 and obtained from logging data. The data are concentrated along the black line, which indicates the degree of correlation between the shearwave velocity obtained from logging data and from seismic interferometry. We also analyze seasonal changes and show the monthly-stacked wavefields at station NIGH13 in Figure 9. The monthly stacked wavefields also have good repeatability between different months.
5.2
Shear-wave splitting
In Figure 10a, we show the wavefields of the shear-wave splitting analysis at station NIGH13 in 2010 that are
obtained by the sequence of deconvolution and cross correlation described in section 4.4. Each trace is plotted at the angle that is equal to the shear-wave polarization used to compute that trace. The thick solid line in Figure 10a shows the interpolated maximum amplitude time of each waveform. The dashed circle shows the arrival time for the wave polarized in the northsouth direction. For the polarizations where the thick solid line is outside of the dashed circle, the shear-wave velocity is slower. The fast and slow shear polarization directions in Figure 10a are 22 degrees and −71 degrees clockwise from the north-south direction, respectively. The angle between these fast and slow directions is 93 degrees, which is close to 90 degrees as predicted by theory (Crampin (1985) and Appendix C). The 3-degree discrepancy could be caused by data noise or discretization errors. At station NIGH13 in 2010, the fast polarization shear-wave velocity vf ast is 638 m/s, the slow velocity vslow is 593 m/s, and the anisotropy parameter (vf ast − vslow )/vf ast is 7% (see Figure 10b). The dif-
322
N. Nakata & R. Snieder 0.0
1500
1000
0.2
Time (s)
Velocity (m/s) in 2008
0.1
500
0.3 0.4 0.5
0 0
500
1000
Velocity (m/s) from log
1500
Figure 8. Crossplot of velocities computed from logging data and by seismic interferometry in 2008. The black line indicates equal velocities.
ference of the arrival times between the fast- and slowpolarization velocity wavefields is much smaller than the period of a wavefield when the borehole depth is less than 525 m.
6
6.1
INTERPRETATION OF SHEAR-WAVE VELOCITIES AND SHEAR-WAVE SPLITTING Influence of major earthquakes
The near-surface shear-wave velocity in Japan is similar between years (see Figure 7), which means the nearsurface structure is basically stable. In this section, we focus on a small region. We use ∆v (calculated by the method presented in section 4.2) and the fast shear-wave polarizations shown in Figure 11 to analyze the influence of major earthquakes in the Niigata prefecture (the dark shaded area in Figure 1). Three significant earthquakes, shown by the dashed arrows in each panel, occurred during the 11 years. Figure 11a shows the velocity variation ∆v for the isotropic component v0 computed by equation 5 compiled over periods one year before and three months after the major earthquakes. We use all stations in the Niigata prefecture and compute the average over the stations. Each box depicts the time range (horizontal extent) and the error in the average velocity over that time interval (vertical length). The error in the velocity is given by the standard deviation of measurements from different earthquakes in each time interval (section 4.2). In Figure 11a, these average velocities show significant
0.6
Jan
Dec
Figure 9. Monthly-stacked wavefields (curves) with the interpolated largest amplitude (circles) at station NIGH13. The horizontal line at around 0.15 s is the shear-wave arrival time obtained from logging data. From left to right, we depict monthly stacks from January to December. Each trace is stacked over the 11 years (January 2000–2010, February 2000–2010, ..., and December 2000–2010). The source and receiver polarization directions are the north-south direction.
velocity reduction after the major earthquakes. The average isotropic velocity of all stations in the region from 2000 to 2010 is 662 m/s, and the relative velocity change of each earthquake is around 3-4%. Similar velocity variations caused by major earthquakes were reported earlier; for example, Sawazaki et al. (2009) analyze the variations caused by the 2000 Western-Tottori Earthquake, Yamada et al. (2010) analyze the variations caused by the 2008 Iwate-Miyagi Nairiku earthquake using KiKnet stations, while Nakata and Snieder (2011) observe a velocity reduction of about 5% after the 2011 TohokuOki earthquake. To increase the temporal resolution of the velocity change, we compute velocity changes averaged over periods one year before and three months after the major earthquakes (Figure 11a) because Sawazaki et al. (2009) found that the velocity reduction is sustained over a period of at least three months after an earthquake. The stations on soft-rock sites have a greater velocity reduction than those on hard-rock sites. (We define soft- and hard-rock sites from the estimated shearwave velocity; hard-rock sites have a shear-wave velocity greater than 600 m/s, while soft-rock sites have a shearwave velocity less than 600 m/s.) For the used eventstation pairs, the velocity reduction does not change measurably with the distance from the epicenter. This is an indication that the velocity reduction depends
Near-surface S-wave velocities in Japan
a)
323
b)
N
W
E
Velocity (m /s)
650
625
600 N
E Polarization direction
S
S
Figure 10. (a) Cross correlograms along every 10-degree polarization direction in 2010 at station NIGH13. Each trace is plotted at an angle equal to the polarization direction used to construct that trace. The dashed circle indicates the peak-amplitude time of the north-south direction, and the thick solid line represents the peak-amplitude time for each polarization direction. (b) Shear-wave velocities computed from the thick solid line in panel (a). Black circles represent the quadratic-interpolated fast and slow polarization shear-wave velocities.
b)
20 15
M6.8, Jul. 16, 2007
10
∆ v (m /s)
5
18 33
0 −5 −10
12 15
−15
29 262
−20 M5.3, Jan. 4, 2001 M6.8, Oct. 23, 2004 −25 2000
2002
2004
2006
2008
2010
Fast shearpolarization (degree)
a)
180 160 140 120 100 80 60 40
29 262 M6.8 Oct. 23, 2004
12 15 20 M5.3 Jan. 4, 2001 0 2000 2002
2004
18 33 M6.8 Jul. 16, 2007 2006 2008 2010
Figure 11. (a) The isotropic component of the shear-wave velocity (Fourier coefficient v0 ) averaged over one year before and three months after three major earthquakes. (b) The direction of the fast shear-wave polarization averaged over same intervals. All data are computed for stations in the Niigata prefecture (the dark shaded area in Figure 1). In each panel, the label of the year is placed in the middle of each year. The dashed arrows in panel (a) and the vertical dashed lines in panel (b) indicate the times of the three major earthquakes shown with the black crosses in Figure 1. The numbers at the left and right of dashed arrows (a) and lines (b) are the number of earthquakes we use for determining velocities and polarizations before and after the major earthquakes, respectively. In panel (a), each velocity is the velocity variation ∆v in equation 4. The horizontal extent of each box depicts the time interval used for averaging (one year before and three months after the major earthquakes). The vertical extent of each box represents the standard deviation of the velocity in the area computed by the method in section 4.2. The horizontal line in each box indicates average velocity ∆v in each time interval, and the vertical line the center of each time interval. In panel (b), we use only the stations with significant anisotropy ((vf ast − vslow )/vf ast ≥ 1%). We select 11 stations from 19 stations and depict the fast shear-wave polarization directions with different symbols or lines. Each symbol is placed at the center of each time period.
mostly on the local geology. The velocity reduction can be due either to the opening and closing of existing fractures, to the creation of new fractures, or to the change in the shear modulus caused by changes in the pore fluid
pressure because of shaking-induced compaction (Das, 1993, Figure 4.24). The relative velocity reduction is smaller than the reduction found by Wu et al. (2009) because of the aver-
N. Nakata & R. Snieder
∆ v (m/s)
a)
b)
25
25
20
20
15
15
10
10
5
5
∆ v (m/s)
324
0 −5
0 −5
−10
−10
−15
−15
−20
−20
−25
Jan
Apr
Jul
c)
−25 50
Oct
100
150 200 Precipitation (mm)
250
300
25 20 15
∆ v (m/s)
10 5 0 −5 −10 −15 −20 −25 50
100
150 200 Precipitation (mm)
250
300
Figure 12. Seasonal dependence of shear-wave velocity. (a) Variation of the average monthly-velocities stacked over the period 2000 through 2010 in southern Japan (the light gray area in Figure 1). We use the stations with the 15% slowest velocities in the area. The horizontal extent of each box shows time interval used for averaging, and the vertical extent the standard deviations of all receivers in the time interval computed by the method in section 4.2. The horizontal line in each box indicates average velocity ∆v in each time interval, and the vertical line the center of each time interval. (b) Crossplot between monthly precipitation (provide by JMA) and the average velocity ∆v with error bars. (c) Crossplot between monthly precipitation and the average velocity ∆v with error bars using the stations with the 85% fastest velocities in the area.
aging over stations and over earthquakes recorded over a period of three months. Wu et al. (2009) use a single station located on a soft-rock site and do not average over several months. Wegler et al. (2009) estimate the velocity reduction in deeper parts of the subsurface, and the velocity reduction they find is small (0.3-0.5 %). From this we infer that the velocity reduction due to a major earthquake is most pronounced in near surface, especially for soft-rock sites. We also obtain the polarization directions of the fast shear waves before and after the major earthquakes by averaging over the same time intervals as used in Figure 11a (Figure 11b). The direction of the fast shearwave polarization does not show a significant change after the earthquakes, hence it seems to be unaffected by these earthquakes. The average standard deviation of the polarization direction of all stations in the Niigata
prefecture between 2000 and 2010 is 15 degrees, which represents the accuracy of the fast shear-wave velocity polarization direction. 6.2
Influence of precipitation
We compute the monthly-averaged shear-wave velocities of the north-south polarization (Figure 12a) to investigate a possible seasonal velocity variation related to precipitation. We use only the data in southern Japan (the light shaded area in Figure 1) because that region has a more pronounced seasonal precipitation cycle than northern Japan. Figure 12a illustrates a significant velocity difference between spring/summer and fall/winter. We calculate the average velocities over the stations with the 15% slowest shear-wave velocities in the area because these stations are located at soft-rock
Near-surface S-wave velocities in Japan a)
325
40 oN
Fast shear polarization 39oN
Plate motion
38 oN
37oN
36 oN
b)
136oE
138oE
140oE
o 142 E
40 oN
Fast shear polarization o
39 N
Plate motion
38 oN Figure 14. Crossplot between the direction of the plate motion and the fast shear-wave polarization directions. We use the stations which have faster than 600 m/s shear-wave velocity. Red indicates there are many points. The north-south direction is 0 degrees, and the east-west direction 90 degrees.
37oN
36 oN
136oE
138oE
140oE
o 142 E
Figure 13. (a) Fast shear-wave polarization directions (black lines) and the direction of the plate motion (gray lines) estimated from GPS data (Sagiya et al., 2000) at the stations with significant anisotropy ((vf ast − vslow )/vf ast ≥ 1%). (b) Extracted stations from panel (a) with shear-wave velocity faster than 600 m/s.
sites and are therefore influenced more by precipitation than the station at hard-rock sites. We compare the monthly-averaged velocities with the monthly average of precipitation (observed by the Japan Meteorological Agency (JMA)) computed from precipitation records over 30 years (Figure 12b). Note the negative correlation between the shear-wave velocity and precipitation (i.e., when it rains, the velocity decreases), which is consistent with the findings of Sens-Sch¨ onfelder and Wegler (2006). For comparison, for the stations with the 85% fastest shear-wave velocities in the area (Figure 12c) the shear-wave velocity does not vary with precipitation. The cause of the velocity reduction is the decreased effective stress of the soil due to the infiltration of water that increases the pore pressure (Das, 1993, Section 4.19; Chapman and Godin, 2001; Snieder and van den Beukel, 2004). We assume that for soft-rock sites most of the velocity change is caused by the effective stress change because Snieder and van den Beukel (2004) show that the relative density change with pore pressure is
much smaller than the relative change in the shear modulus. 6.3
Shear-wave splitting and the direction of the plate motion
Using shear-wave splitting analysis, we determine fast shear-wave polarization directions of every station (illustrated by the black arrows in Figure 13a). These directions are averaged over all years from 2000 to 2010 because the temporal changes in the direction are small (see Figure 11b). We plot the directions of all stations which have an anisotropy parameter (vf ast − vslow )/vf ast ≥ 1% because the uncertainty in the direction of the fast shear polarization is large when the anisotropy is small. In Figure 13a, we also plot the direction of the plate motion at each station (the gray arrows), estimated from GPS data (Sagiya et al., 2000). Each arrow is normalized to the same length. In Figure 13b, we plot only the stations which have an anisotropy parameter larger than 1% and a northsouth polarization shear-wave velocity faster than 600 m/s; these stations are located on hard-rock sites. The average absolute angle between the directions of fast shear polarization and the plate motion in these stations is 16 degrees, and this average angle of the stations which have a shear-wave velocity less than 600 m/s is 36 degrees. Therefore, the fast shear-wave polarization on hard-rock sites correlates more strongly with the direc-
326
N. Nakata & R. Snieder
tion of the plate motion than the polarization on softrock sites. The 16-degree angle is close to the 15-degree standard deviation angle of each station computed in section 6.1. The near-surface polarization in the western part of Figure 13b correlates well with observations of the shear-wave polarization at greater depth (Okada et al., 1995; Nakajima and Hasegawa, 2004; Nakajima et al., 2006), but this agreement does not hold in the regions further east. We present a crossplot of the directions of the fast shear-wave polarization and the plate motion for stations all over Japan (Figure 14), where we only used the stations which have an anisotropy parameter greater than 1% and a shear-wave velocity faster than 600 m/s. The red area in Figure 14 indicates that for most stations the direction of the plate motion is between 90 and 140 degrees, and that this direction correlates with the polarization direction of the fast shear wave. The near-surface stress directions on hard-rock sites is presumably related to the plate motion because the stress field related to the plate motion changes the properties of fractures. Note that the used shear waves sample the shallow subsurface (down to about several hundreds of meters). It is remarkable that the shear-wave velocities in the near surface at hard-rock sites correlate with tectonic process (plate motion) that extends several tens of kilometers into the subsurface.
7
CONCLUSION
We obtain annual and monthly averaged near-surface shear-wave velocities throughout Japan by applying seismic interferometry to KiK-net data. Deconvolution interferometry yields more repeatable and higher resolution wavefields than does cross-correlation interferometry. Because picked arrival times in waveforms are generally stable over time and consistent with logging data, the near-surface has a stable shear-wave velocity. After three strong earthquakes in the Niigata prefecture, however, the shear-wave velocity is reduced. By computing the monthly-stacked velocity, we observe a velocity variation on stations placed on soft rock that has a negative correlation with precipitation. We also observe shear-wave splitting. The fast shear-wave polarization direction on a hard-rock site correlates with the direction of the plate motion. Because the shearwave velocity is related to ground soil strength, these velocities are useful for civil engineering, site characterization, and disaster prevention.
8
ACKNOWLEDGMENTS
We are grateful to NIED for providing us with the KiKnet data and to JMA for precipitation data. We thank Diane Witters for her professional help in preparing this
manuscript. We are grateful to Kaoru Sawazaki for suggestions, corrections, and discussions.
REFERENCES Aki, K., 1957, Space and time spectra of stationary stochastic waves, with special reference to microtremors.: Bull. of the Earthquake Res. Inst. Univ. Tokyo, 35, 415–456. Alford, R. M., 1986, Shear data in the presence of azimuthal anisotropy: Dilley, Texas: SEG Expanded Abstracts, 5, 476–479. Bakulin, A., and R. Calvert, 2004, Virtual source: new method for imaging and 4D below complex overburden: SEG Expanded Abstracts, 23, 112–115. ——–, 2006, The virtual source method: Theory and case study: Geophysics, 71, S1139–S1150. Bakulin, A., and A. Mateeva, 2008, Estimating interval shear-wave splitting from multicomponent virtual shear check shots: Geophysics, 73, A39–A43. Brenguier, F., M. Campillo, C. Hadziioannou, N. M. Shapiro, R. M. Nadeau, and E. Larose, 2008, Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations: Science, 321, 1478–1481. Chapman, D., and O. Godin, 2001, Dispersion of interface waves in sediments with power-law shear speed profiles. II. Experimental observations and seismoacoustic inversions: J. Acoust. Soc. Am., 110, 1908– 1916. Claerbout, J. F., 1968, Synthesis of a layered medium from its acoustic transmission response: Geophysics, 33, 264–269. Crampin, S., 1985, Evaluation of anisotropy by shearwave splitting: Geophysics, 50, 142–152. Das, B., 1993, Principles of soil dynamics: PWS-Kent Publishing Co. Draganov, D., X. Campman, J. Thorbecke, A. Verdel, and K. Wapenaar, 2009, Reflection images from ambient seismic noise: Geophysics, 74, A63–A67. Draganov, D., K. Wapenaar, W. Mulder, J. Singer, and A. Verdel, 2007, Retrieval of reflections from seismic background-noise measurements: Geophys. Res. Lett., 34, L04305. Hohl, D., and A. Mateeva, 2006, Passive seismic reflectivity imaging with ocean-bottom cable data: SEG Expanded Abstracts, 25, 1560–1564. Larose, E., L. Margerin, A. Derode, B. van Tiggelen, M. Campillo, N. Shapiro, A. Paul, L. stehly, and M. Ta, 2006, Correlation of random wavefields: An interdisciplinary review: Geophysics, 71, S111–S121. Lawson, C. L., 1984, C1 surface interpolation for scattered data on a sphere: Rocky Mountain J. Math., 14, 177–202. Li, Y.-G., J. E. Vidale, K. Aki, F. Xu, and T. Burdette, 1998, Evidence of shallow fault zone strengthening
Near-surface S-wave velocities in Japan after the 1992 M7.5 Landers, California, earthquake: Science, 279, 217–219. Lin, F.-C., M. H. Ritzwoller, and R. Snieder, 2009, Eikonal tomography: surface wave tomography by phase front tracking across a regional broad-band seismic array: Geophys. J. Int., 177, 1091–1110. Lobkis, O. I., and R. L. Weaver, 2001, On the emergence of the Green’s function in the correlations of a diffuse field: J. Acoust. Soc. Am., 110, 3011–3017. Ma, S., G. A. Prieto, and G. C. Beroza, 2008, Testing community velocity models for southern California using the ambient seismic field: Bull. Seismol. Soc. Am., 98, 2694–2714. Mehta, K., J. L. Sheiman, R. Snieder, and R. Calvert, 2008, Strengthening the virtual-source method for time-lapse monitoring: Geophysics, 73, S73–S80. Mehta, K., R. Snieder, and V. Graizer, 2007a, Downhole receiver function: a case study: Bull. Seismol. Soc. Am., 97, 1396–1403. ——–, 2007b, Extraction of near-surface properties for a lossy layered medium using the propagator matrix: Geophys. J. Int., 169, 271–280. Minato, S., T. Matsuoka, T. Tsuji, D. Draganov, J. Hunziker, and K. Wapenaar, 2011, Seismic interferometry using multidimensional deconvolution and crosscorrelation for crosswell seismic reflection data without borehole sources: Geophysics, 76, SA19– SA34. Miyazawa, M., R. Snieder, and A. Venkataraman, 2008, Application of seismic interferometry to exptract P- and S-wave propagation and observation of shear-wave splitting from noise data at Cold Lake, Alberta, Canada: Geophysics, 73, D35–D40. Nakajima, J., and A. Hasegawa, 2004, Shear-wave polarization anisotropy and subduction-induced flow in the mantle wedge of northeastern Japan: Earth Planet. Sci. Lett., 225, 365–377. Nakajima, J., T. Matsuzawa, A. Hasegawa, and D. Zhao, 2001, Three-dimensional structure of Vp , Vs , and Vp /Vs beneath northeastern Japan: Implications for arc magmatism and fluids: J. Geophys. Res., 106, 21843–21857. Nakajima, J., J. Shimizu, S. Hori, and A. Hasegawa, 2006, Shear-wave splitting beneath the southwestern Kurile arc and northeastern Japan arc: A new insight into mantle return flow: Geophys. Res. Lett., 33, L05305. Nakata, N., and R. Snieder, 2011, Near-surface weakening in Japan after the 2011 Tohoku-Oki earthquake: Geophys. Res. Lett., 38, L17302. Nakata, N., R. Snieder, T. Tsuji, K. Larner, and T. Matsuoka, 2011, Shear-wave imaging from traffic noise using seismic interferometry by cross-coherence: Geophysics, 76, SA97–SA106. Okada, T., T. Matsuzawa, and A. Hasegawa, 1995, Shear-wave polarization anisotropy beneath the north-eastern part of Honshu, Japan: Geophys. J.
327
Int., 123, 781–797. Okada, Y., K. Kasahara, S. Hori, K. Obara, S. Sekiguchi, H. Fujiwara, and A. Yamamoto, 2004, Recent progress of seismic observation networks in Japan — Hi-net, F-net, K-NET and KiK-net —: Earth Planets Space, 56, 15–28. Prieto, G. A., J. F. Lawrence, and G. C. Beroza, 2009, Anelastic earth structure from the coherency of the ambient seismic field: J. Geophys. Res., 114, B07303. Roux, P., and M. Fink, 2003, Green’s function estimation using secondary sources in a shallow wave environment: J. Acoust. Soc. Am., 113, 1406–1416. Ruigrok, E., X. Campman, D. Draganov, and K. Wapenaar, 2010, High-resolution lithospheric imaging with seismic interferometry: Geophys. J. Int., 183, 339– 357. Sagiya, T., S. Miyazaki, and T. Tada, 2000, Continuous GPS array and present-day crustal deformation of Japan: Pure appl. geophy., 157, 2303–2322. Sawazaki, K., H. Sato, H. Nakahara, and T. Nishimura, 2009, Time-lapse changes of seismic velocity in the shallow ground caused by strong ground motion shock of the 2000 western-Tottori earthquake, Japan, as revealed from coda deconvolution analysis: Bull. Seismol. Soc. Am., 99, 352–366. Schaff, D. P., and G. C. Beroza, 2004, Coseismic and postseismic velocity changes measured by repeating earthquakes: J. Geophys. Res., 109, B10302. Schuster, G. T., J. Yu, J. Sheng, and J. Rickett, 2004, Interferometric/daylight seismic imaging: Geophys. J. Int., 157, 838–852. Sens-Sch¨ onfelder, C., and U. Wegler, 2006, Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia: Geophys. Res. Lett., 33, L21302. Shiomi, K., K. Obara, S. Aoi, and K. Kasahara, 2003, Estimation on the azimuth of the Hi-net and KiKnet borehole seismometers (in japanese): Zishin2, 56, 99–110. Snieder, R., and E. S ¸ afak, 2006, Extracting the building response using seismic interferometry: Theory and application to the Millikan library in Pasadena, California: Bull. Seismol. Soc. Am., 96, 586–598. Snieder, R., and A. van den Beukel, 2004, The liquefaction cycle and the role of drainage in liquefaction: Granular Matter, 6, 1–9. Snieder, R., K. Wapenaar, and K. Larner, 2006, Spurious multiples in seismic interferometry of primaries: Geophysics, 71, SI111–SI124. Stehly, L., M. Campillo, B. Froment, and R. L. Weaver, 2008, Reconstruction Green’s function by correlation of the coda of the correlation (C3 ) of ambient seismic noise: J. Geophys. Res., 113, B11306. Thomsen, L., 1988, Reflection seismology over azimuthally anisotropic media: Geophysics, 53, 304– 313. Tonegawa, T., K. Nishida, T. Watanabe, and K. Sh-
328
N. Nakata & R. Snieder
iomi, 2009, Seismic interferometry of teleseismic Swave coda retrieval of body waves: an application to the Philippine Sea slab underneath the Japanese Islands: Geophys. J. Int., 178, 1574–1586. Trampert, J., M. Cara, and M. Frogneux, 1993, SH propagator matrix and Qs estimates from boreholeand surface- recorded earthquake data: Geophys. J. Int., 112, 290–299. Vasconcelos, I., and R. Snieder, 2008, Interferometry by deconvolution: Part 2 - Theory for elastic waves and application to drill-bit seismic imaging: Geophysics, 73, S129–S141. Vidale, J. E., and Y.-G. Li, 2003, Damage to the shallow Landers fault from the nearby Hector Mine earthquake: Nature, 421, 524–526. Wapenaar, K., 2004, Retrieving the elastodynamic Green’s function of an arbitrary inhomogeneous medium by cross correlation: Phys. Rev. Lett., 93, 254–301. Wapenaar, K., D. Draganov, R. Snieder, X. Campman, and A. Verdel, 2010, Tutorial on seismic interferometry: Part 1 - Basic principles and applications: Geophysics, 75, 75A195–75A209. Wapenaar, K., and J. Fokkema, 2006, Green’s function representations for seismic interferometry: Geophysics, 71, SI33–SI46. Wapenaar, K., J. van der Neut, and E. Ruigrok, 2008, Passive seismic interferometry by multidimensional deconvolution: Geophysics, 73, A51–A56. Wegler, U., H. Nakahara, C. Sens-Sch¨ onfelder, M. Korn, and K. Shiomi, 2009, Sudden drop of seismic velocity after the 2004 Mw 6.6 mid-Niigata earthquake, Japan, observed with passive image interferometry: J. Geophys. Res., 114, B06305. Wegler, U., and C. Sens-Sch¨ onfelder, 2007, Fault zone monitoring with passive image interferometry: Geophys. J. Int., 168, 1029–1033. Wu, C., Z. Peng, and D. Assimaki, 2009, Temporal changes in site response associated with the strong ground motion of the 2004 MW 6.6 Mid-Niigata earthquake sequences in Japan: Bull. Seismol. Soc. Am., 99, 3487–3495. Yamada, M., J. Mori, and S. Ohmi, 2010, Temporal changes of subsurface velocities during strong shaking as seen from seismic interferometry: J. Geophys. Res., 115, B03302.
APPENDIX A: 1D SEISMIC INTERFEROMETRY We explain in this appendix why deconvolution interferometry is suitable for this study. Figure A1 illustrates the model of interferometry using a KiK-net station and an earthquake. The incoming wavefield S(rb , s, ω), propagating from source s to receiver rb , is given by G(rb , s, ω)W (s, ω), where G is the Green’s function including any unknown complex effect of wave propaga-
z=0
rs
z=z b
rb
s Figure A1. Geometry of an earthquake and a KiK-net station, where rs is the surface receiver (black triangle), rb the borehole receiver (white triangle), and s the epicenter of earthquake (gray star).
tion such as scattering and attenuation, and W the source signature in the frequency domain. Assuming that the subsurface is homogeneous between the receivers, the received wavefield at surface receiver rs is 2S(rb , s, ω)ejkzb e−γzb , where γ is the attenuation coefficient and k the wave number. Because of the free surface, the amplitude of the wavefield at the surface is multiplied by a factor 2. We assume that there are no multiples between the two receivers. The reflected wavefield from the surface at the borehole receiver rb is S(rb , s, ω)e2jkzb e−2γzb , and the total wavefield at rb is S(rb , s, ω) + S(rb , s, ω)e2jkzb e−2γzb . Applying crosscorrelation interferometry to these wavefields yields
u(rs , s, ω)u∗ (rb , s, ω) = 2S(rb , s, ω)ejkzb e−γzb h i S ∗ (rb , s, ω) + S ∗ (rb , s, ω)e−2jkzb e−2γzb ≈ 2|S(rb , s, ω)|2 ejkzb e−γzb
(A1)
when we consider only the first arrival. Using deconvolution interferometry, we obtain
Near-surface S-wave velocities in Japan [h]
in the horizontal direction. In our study, the virtual source in the borehole has the polarization of the incident wave. The two horizontal components of the virtual source therefore are not independent, so that Alford rotation cannot be applied. The angle of the fast and slow shear-wave polarization directions is 90 degrees because we assume the incoming wave is a plane wave (Crampin, 1985). A waveˆ which is a unit vector, can be field with polarization p, expressed in the polarization of the fast and slow shear waves:
pˆ s
pˆ
pˆ s sin φ
pˆ = pˆf cos φ + pˆs sin φ,
φ pˆ f cos φ
pˆ f
Figure B1. Projection of fast and slow velocity directions, ˆf is the fast polarization direction, p ˆs the slow polarwhere p ˆ an arbitrary direction, and φ the angle ization direction, p ˆ p ˆf , between the fast direction and arbitrary direction. p, ˆs are unit vectors. Dashed arrows show the projection, and p which is shown in equation B1.
jkzb −γzb
e 2S(rb , s, ω)e u(rs , s, ω) = u(rb , s, ω) S(rb , s, ω) + S(rb , s, ω)e2jkzb e−2γzb ≈ 2ejkzb e−γzb ,
329
(A2)
where we also only retain the first arrival. The planewave Green’s function excited at rb and received at rs is equal to 2ejkzb e−γzb in the frequency domain. Let us compare equations A1, A2, and the Green’s function. The wavefield retrieved by cross-correlation interferometry is complicated because equation A1 includes the power spectrum |S(rb , s, ω)|2 . This term is different for different earthquakes. In contrast, deconvolution interferometry eliminates the incoming wave S(rb , s, ω), and thus provides a more accurate estimate of the Green’s function. When we stack deconvolved wavefields over earthquakes, the accuracy of this estimate is improved. Because the deconvolved waves do not depend on the power spectrum of the incident wave |S(rb , s, ω)|2 , the deconvolved wavefields are more reproducible than those obtained from cross correlation.
APPENDIX B: SHEAR-WAVE SPLITTING Usually, shear-wave splitting is analyzed with Alford rotation (Alford, 1986; Thomsen, 1988). This procedure is based on the use of two independent orthogonal sources
(B1)
where φ is the polarization angle in the arbitrary wavefield relative to the direction of the fast shear-wave polarization, and pˆf and pˆs are the unit vectors of the fast and slow velocity wavefields, respectively (see Figure B1). The incoming wavefield ub at the borehole receiver is
ub = S(t)pˆ = S(t)pˆf cos φ + S(t)pˆs sin φ,
(B2)
where S(t) is the incoming wavefield, and t the time. The wavefield at surface receiver us is given by „ « „ « zb zb us = S t − pˆf cos φ + S t − pˆs sin φ, (B3) vf vs where vf is the fast velocity, vs the slow velocity, and zb the distance between the top and bottom receivers. The component of us along pˆ is
us (φ) =(pˆ · us ) „ « zb =S t − (pˆ · pˆf ) cos φ vf „ « zb +S t− (pˆ · pˆs ) sin φ vs « „ « „ zb zb cos2 φ + S t − =S t − sin2 φ. vf vs
(B4)
We express vf and vs in the average velocity v0 and the difference δ: 1 1 = − δ, vf v0
1 1 = + δ, vs v0
(B5)
and assume the splitting time zb δ is much smaller than the period. We insert expression B5 into equation B4 and expand using first-order Taylor expansion in δ:
330
N. Nakata & R. Snieder „
«
zb us (φ) =S t − + zb δ cos2 φ v0 „ « zb − zb δ sin2 φ +S t− v0 „ « zb =S t − (sin2 φ + cos2 φ) v0 „ « zb 0 +S t− zb δ(cos2 φ − sin2 φ) v0 „ « „ « zb zb =S t − + S0 t − zb δ cos 2φ v0 v0 « „ zb + zb δ cos 2φ , (B6) =S t − v0 where S 0 is the time derivative of S. Thus, using Taylor expansion, we obtain the velocity for a shear wave with the polarization of equation B1: v(φ) =
v0 , 1 − v0 δ cos 2φ
(B7)
or to first order in v0 δ: v(φ) = v0 (1 + v0 δ cos 2φ).
(B8)
APPENDIX C: FOURIER COEFFICIENTS We describe the meaning of v0 , v1 , and v2 in equation 5. Expression B8 gives the velocity for a polarization φ relative to the polarization of the fast shear wave. When the azimuth of the fast shear-wave polarization is given by ψ, the angle φ in equation B8 must be changed into φ → φ − ψ. Denoting v02 δ by V , these changes turn equation B8 into v(φ) = v0 + V cos 2(φ − ψ).
(C1)
This can also be written as v(φ) = v0 + v1 cos 2φ + v2 sin 2φ;
(C2)
with
Z π Z 1 1 π v(φ)dφ, v1 = v(φ) cos 2φdφ, 2π −π π −π Z π 1 v2 = v(φ) sin 2φdφ, π −π „ « q 1 v2 V = v12 + v22 , ψ = arctan . 2 v1 v0 =
In expression C2, v0 , v1 , and v2 are the Fourier coefficients of the velocity v(φ). Because v0 does not depend on φ, v0 represents the isotropic velocity. Expression C1 shows that v(ψ) and v(ψ + π/2) are the fastest and slowest velocities, respectively. Therefore, V is the anisotropic velocity and the angle between the fast and slow polarization directions is 90 degrees,
which corresponds to shear-wave splitting of a plane wave (Crampin, 1985).
Estimating near-surface shear-wave velocities in Japan by applying seismic interferometry to KiK-net data Nori Nakata & Roel Snieder Center for Wave Phenomena, Geophysics Department, Colorado School of Mines, Golden, Colorado 80401
ABSTRACT
We estimate shear-wave velocities in the shallow subsurface throughout Japan by applying seismic interferometry to the data recorded with KiK-net, a strongmotion network in Japan. Each KiK-net station has two receivers; one receiver on the surface and the other in a borehole. By using seismic interferometry, we extract the shear wave that propagates between these two receivers. Applying this method to earthquake-recorded data at all KiK-net stations from 2000 to 2010 and measuring the arrival time of these shear waves, we analyze monthly and annual averages of the near-surface shear-wave velocity all over Japan. Shear-wave velocities estimated by seismic interferometry agree well with the velocities obtained from logging data. The estimated shear-wave velocities of each year are stable. For the Niigata region, we observe a velocity reduction caused by major earthquakes. For stations on soft rock, the measured shear-wave velocity varies with the seasons, and we show negative correlation between the shear-wave velocities and precipitation. We also analyze shear-wave splitting by rotating the horizontal components of the surface sensors and borehole sensors and measuring the dependence on the shear-wave polarization. This allows us to estimate the polarization with the fast shear-wave velocity throughout Japan. For the data recorded at the stations built on hard-rock sites, the fast shearwave polarization directions correlate with the direction of the plate motion. Key words: seismic interferometry, shear-wave splitting, shear-wave velocity, precipitation, near surface, KiK-net
1
INTRODUCTION
Seismic interferometry is a powerful tool to obtain the Green’s function that describes wave propagation between two receivers (e.g., Aki, 1957; Claerbout, 1968; Lobkis and Weaver, 2001; Roux and Fink, 2003; Schuster et al., 2004; Wapenaar et al., 2004; Bakulin and Calvert, 2006; Snieder et al., 2006; Wapenaar and Fokkema, 2006). Seismic interferometry is applied to ambient noise (e.g., Hohl and Mateeva, 2006; Draganov et al., 2007, 2009; Brenguier et al., 2008; Stehly et al., 2008; Lin et al., 2009), traffic noise (e.g., Nakata et al., 2011), production noise (e.g., Miyazawa et al., 2008; Vasconcelos and Snieder, 2008), earthquake data (e.g., Larose et al., 2006, Sens-Sch¨ onfelder and Wegler, 2006; Snieder and S ¸ afak, 2006; Ma et al., 2008; Ruigrok et al., 2010), and active sources (e.g., Bakulin and Calvert, 2004; Mehta et al., 2008). In Japan, large seismometer networks, such as Hinet, F-net, K-NET, and KiK-net (Okada et al., 2004),
are deployed. By using these networks for seismic interferometry, Tonegawa et al. (2009) extract the deep subsurface structure of the Philippine Sea slab. These data have also been used to observe time-lapse changes in small regions (Wegler and Sens-Sch¨ onfelder, 2007; Sawazaki et al., 2009; Wegler et al., 2009; Yamada et al., 2010). Each KiK-net station has two receivers, one on the ground surface and the other at the bottom of a borehole. One can estimate the body-wave velocity between two receivers by using seismic interferometry (Trampert et al., 1993; Snieder and S ¸ afak, 2006; Mehta et al., 2007b; Miyazawa et al., 2008). By applying seismic interferometry to KiK-net data, we analyze near-surface velocities throughout Japan. Because KiK-net has recorded strong-motion seismograms continuously since the end of 1990s, the data are available for time-lapse measurements. Measuring time-lapse changes of the shallow subsurface is important for civil engineering and for estimating the site response to earthquakes. Previous studies ex-
316
N. Nakata & R. Snieder
a) 48 oN
b)
o 39 N
Jul. 16, 2007 M6.8
o
44 N o 38 N
o
Oct. 23, 2004 M6.8
NIGH13
40 oN o 37 N
36 N
Jan. 4, 2001 M5.3
32 oN
o 36 N o 137 E o o 128oE 132oE 136oE 140oE 144 E 148 E
o 138 E
o 139 E
o 140 E
Figure 1. KiK-net stations (December, 2010). The black dots on the map represent the locations of the stations. The dark gray shows the area analyzed in section 6.1. The light gray illustrates the area where we apply the analysis for seasonal change (section 6.2). The rectangle area in panel (a) is magnified in panel (b). The large black circle indicates station NIGH13, which we use for examples of analysis in Figures 2, 4, 5, 6, 9, and 10. The black crosses depict the epicenters of three significant earthquakes that occurred in the vicinity.
tracted time-lapse changes caused by earthquakes (Li et al., 1998; Vidale and Li, 2003; Schaff and Beroza, 2004; Wegler and Sens-Sch¨ onfelder, 2007; Brenguier et al., 2008). Interferometry applied to a single KiKnet station has also been used to measure time-lapse change due to earthquakes (Sawazaki et al., 2009; Yamada et al., 2010). Interferometric studies have shown changes in the shear-wave velocity caused by precipitation (Sens-Sch¨ onfelder and Wegler, 2006) and have measured shear-wave splitting (Bakulin and Mateeva, 2008; Miyazawa et al., 2008). We study the annual and monthly averages of the shear-wave velocity and the fast shear-wave polarization directions for stations all over Japan, and the temporal change in shear-wave velocity in the Niigata prefecture for three major earthquakes. This paper presents data processing of KiK-net data with seismic interferometry. We first introduce the properties of KiK-net. Next, we show the data analysis method. Then we present near-surface shear-wave velocities in every part of Japan. Finally, we interpret these velocities to study time-lapse changes, which are related to major earthquakes and precipitation, and present measurements of shear-wave splitting.
2
KIK-NET DATA
About 700 KiK-net stations are distributed in Japan (Figure 1). The stations are operated by the National
Research Institute for Earth Science and Disaster Prevention (NIED). Each station has a borehole and two seismographs which record strong motion at the bottom and top of the borehole. Each seismograph has three components: one vertical component and two horizontal components. Although the two horizontal components of the surface seismograph are oriented in the northsouth and east-west directions, respectively, the horizontal components of the borehole seismograph are not always aligned with the north-south and east-west directions because of technical limitations. Therefore, we rotate the directions of the borehole seismograph northsouth and east-west directions before data processing. The depth of about 25% of the boreholes is 100 m, and the other boreholes are at greater depth. Since our target is the near surface, we use the stations with a depth less than 525 m, which accounts for 94% of the stations. The sampling interval is either 0.005 or 0.01 s, depending on the station and the recording date. We show example records of an earthquake in Figure 2. Figure 2a illustrates bandpass-filtered time series, and Figure 2b the power spectra of the unfiltered data. As shown in Figure 2b, most energy is confined to 1-13 Hz, and we apply a bandpass filter over this frequency range for all data processing. In Figure 2, UD denotes the vertical component, NS the north-south direction horizontal component, and EW the east-west direction horizontal component. In Figure 2a, the P wave arrives
Near-surface S-wave velocities in Japan a)
surface
UD2
borehole
UD1
borehole
NS1
3
RETRIEVAL OF THE WAVEFIELD BETWEEN RECEIVERS
surface
EW 2
borehole
EW 1 10
20
30
40
50
Time (s)
b) UD2
surface
UD1
borehole
NS2
surface
NS1
borehole
EW 2
surface
EW 1
borehole
0
shear-wave velocity between these sensors as determined in this study. A bias in the velocity estimation due to non-vertical propagation depends on the deviation of cos θ from its value for vertical incidence, cos 0◦ = 1.
surface
NS2
0
317
10
20
30
40
50
Frequency (Hz) Figure 2. An earthquake recorded at all channels of station NIGH13 (latitude 37.0514◦ N and longitude 138.3997◦ E). This earthquake occurred at 14:59:19.56, 27 October 2004. The epicenter is at latitude 37.2204◦ N and longitude 138.5608◦ E and the depth is 11.13 km. The magnitude of this earthquake is MJM A 4.2. UD represents the vertical component, NS the north-south direction horizontal component, EW the east-west direction horizontal component, 1 the borehole seismograph, and 2 the surface seismograph. (a) The bandpass-filtered (1-13 Hz) time series. (b) The power spectra of the unfiltered records.
at around 7 s, and the shear wave arrives at around 14 s. All the used events are at a depth greater than 10 km. Because of this large depth compared to the depth of the boreholes and the low velocity in the near surface, the waves that travel between the sensors at each station propagate in the near-vertical direction as plane waves. We compute the angle of the incoming wave at the borehole receiver by using one-dimensional ray tracing to confirm that the wave propagating between the borehole and surface seismometers, propagates in the near-vertical direction. We use the velocity model of Nakajima et al. (2001) to determine the ray parameter p of the ray between each earthquake and the borehole sensor. The angle of incidence θ of the wave propagating from the borehole p receiver to the surface receiver is given by cos θ = 1 − p2 v 2 , where v is the average
We apply seismic interferometry to the recorded earthquake data of each station for retrieving the wavefield where the borehole receiver behaves as a virtual source. Several algorithms have been used in seismic interferometry to obtain the wavefield. These include cross correlation (e.g., Claerbout, 1968; Bakulin and Calvert, 2004), deconvolution (e.g., Trampert et al.), 1993; Snieder and S ¸ afak, 2006), cross coherence (e.g., Aki, 1957; Prieto et al., 2009), and multidimensional deconvolution (e.g., Wapenaar et al., 2008; Minato et al., 2011). We introduce the cross-correlation and deconvolution algorithms. We denote the wavefield, excited at source location s that strikes the borehole receiver at location rb by u(rb , s, ω) = S(rb , s, ω), where S(rb , s, ω) is the incoming wavefield that includes the source signature of the earthquake and the effect of propagation such as attenuation and scattering, in the frequency domain. The corresponding wavefield recorded at the surface receiver at location rs is given by u(rs , s, ω) = 2G(rs , rb , ω)S(rb , s, ω),
(1)
where the factor 2 is due to the presence of the free surface at rs . Because the wavefield striking the borehole receiver is close to a vertically propagating plane wave, G(rs , rb , ω) is the plane wave Green’s function that accounts for the propagation from the borehole seismometer to the surface seismometer. The cross-correlation approach to retrieve the wavefield in one dimension is given by Wapenaar et al. (2010) |S(rb , s, ω)|2 G(rs , rb , ω) =
2jω u(rs , s, ω)u∗ (rb , s, ω), ρc (2)
where ρ is the mass density of the medium, c the wave propagation velocity, j the imaginary unit, and ∗ the complex conjugate. The regularized deconvolution, which is similar to cross correlation, is given by
G(rs , rb , ω) =
u(rs , s, ω) u(rs , s, ω)u∗ (rb , s, ω) ≈ , u(rb , s, ω) |u(rb , s, ω)|2 + (3)
where is a regularization parameter (Mehta et al., 2007b,a). The deconvolution is potentially unstable due to the spectral devision, and we avoid divergence by adding a positive constant to the denominator (equation 3). Note that the deconvolution eliminates the
318
N. Nakata & R. Snieder
Time (s)
0.14
0.16
0.18
Figure 3. Quadratic interpolation. Using the three largest amplitude points (crosses), we interpolate the highest amplitude point (circle) by estimating the quadratic curve through the three highest amplitude points.
empirically that this is the smallest regularization parameter to obtain stable wavefields. We apply a bandpass filter from 1-13 Hz after applying deconvolution interferometry. In this paper, we average in three ways to interpret the wavefields. The first method is the annual stack, where we average the deconvolved waveforms over the earthquakes recorded in each year. In the second averaging method, we average the deconvolved waves over all earthquakes recorded in each month over the 11 years (January 2000–2010, February 2000–2010, ..., and December 2000–2010). We call this average the monthly stack. In the third method, which we use for analyzing the influence of major earthquakes, we average over three months after a significant earthquake.
4.1 imprint of waveform S(rb , s, ω), which is incident on the borehole receiver. We derive the features of crosscorrelation and deconvolution interferometry in Appendix A.
4
DATA PROCESSING
We use 111,934 earthquake-station pairs that are recorded between 2000 and 2010. The magnitude range is confined between 1.9 and 8.2. The cosine of the angle of incidence cos θ of the wave propagating between the receivers at each station is greater than 0.975, even for the events that are the furthest away. The bias introduced by non-vertical propagation thus is less than 2.5 %, and for most measurements it is much smaller. First, we check the data quality and drop some seismograms by a visual inspection using the signal-to-noise ratio as a criterion. Additionally, we discard stations with a borehole seismometer at a depth greater than 525 m because we focus this study on the near surface. We remove the DC component of the data by subtracting the average of each seismogram. For aligning the directions of the borehole receiver to the exact north-south and east-west directions, we rotate the borehole receiver using the rotation angle provided by NIED (Shiomi et al., 2003). Because the sampling interval is not small compared to the travel time of P waves between the borehole and surface seismometers, we focus on the shear wave and only analyze the horizontal components. We first apply deconvolution interferometry to the motion in the north-south direction of each surfaceborehole pair. In this study, for reasons explained below, deconvolution interferometry gives more consistent estimates of the Green’s function than does crosscorrelation interferometry. We choose in equation 3 to be 1% of the average power spectrum of the borehole receiver in the frequency range 1-13 Hz because we find
Estimating the shear-wave velocity
Before we apply annual stacking or monthly stacking, we resample the data from 0.005-s interval to 0.01-s interval if the data that are stacked include both 0.005-s and 0.01-s sampling-interval data. After stacking, we estimate the arrival time by seeking the three adjacent samples with the largest values and apply quadratic interpolation to find the time at which the deconvolved data have a maximum amplitude (Figure 3). This time is the travel time for a shear wave that propagates between the borehole and surface sensors. We use this travel time to compute the shear-wave velocity of the region between the two receivers.
4.2
Computing the average and standard deviation of the velocity of the annual or monthly stacks
To interpret time-lapse variations in the velocities, we need to compute the average and standard deviation of the velocities within a region. Let us denote the estimated velocity by vi (m, y), where vi is the shear-wave velocity at station i, in month m, and year y. This velocity is already averaged over each month. Each station has a different velocity. In order to quantify the timelapse variations of the velocity, we subtract the average value of each station before calculating temporal variation in the annual or monthly average: ∆vi (m, y) = vi (m, y) − v i ,
(4)
where v i is an average velocity of station i over all months and years. Then we compute either the annual or monthly average of the velocity variation over stations ∆v, and we also compute the standard deviation of this quantity.
Near-surface S-wave velocities in Japan
0.1
0.1
0.2
0.2
Tim e (s)
b) 0.0
Tim e (s)
a) 0.0
0.3
0.3
0.4
0.4
0.5
0.5
0.6
2000
0.6
2010
319
2000
2010
Figure 4. Annual-stacked wavefields by using (a) cross correlation and (b) deconvolution interferometry at station NIGH13. The surface and borehole receiver orientation directions are north-south. Epicenter locations are illustrated in Figure 5.
b) 48 oN
48 oN
100
75
44 oN
44 oN
2005: 19 40 oN
40 oN
2000: 5
2002: 9
36 oN
2004: 110
32 oN
o o 128 oE 132oE 136oE 140oE 144 E 148 E
25
2008: 10 2009: 14
2003: 5 32 oN
2006: 5 2007: 32
2001: 12 36 oN
50
Depth (km )
a)
2010: 10
0
3 5 7
o Magnitude o 128 oE 132oE 136oE 140oE 144 E 148 E
Figure 5. Epicenters used in Figure 4 from 2000 to 2004 (a) and from 2005 to 2010 (b) at station NIGH13. At the right in each panel, we show the number of earthquakes we use to obtain the waveforms in Figure 4 in each year. The size of each circle refers the magnitude of each earthquake and the color indicates the depth. The white triangle illustrates the location of station NIGH13. Because of the proximity of events, many circles overlap.
4.3
Comparison between cross-correlation and deconvolution interferometry
We compare the cross-correlation and deconvolution approaches using the annual-stacked wavefields (Figure 4). We show the locations of the epicenters of the used earthquakes in Figure 5. The annual stacks of the waveforms obtained by cross correlation are shown in Figure 4a. These waveforms are not repeatable from year to year and often do not show a pronounced peak at the
arrival time of the shear wave at around t = 0.15 s. We attribute the variability in these waveforms to variations in the power spectrum |S(rb , s, ω)|2 of the waves incident at the borehole receiver (equation A1). In contrast, the annual stacks of the waveforms obtained by deconvolution shown in Figure 4b are highly repeatable and show a consistent peak at the arrival time of the shear wave. The consistency of these waveforms is due to the deconvolution that eliminates the imprint of the incident wave S(rb , s, ω) (equation A2). Consistent with
320
N. Nakata & R. Snieder
earlier studies (Trampert et al., 1993; Snieder and S ¸ afak, 2006), we use deconvolution to extract the waves that propagate between the seismometers at each KiK-net station.
0.0
4.4
0.2
v(φ) = v0 + v1 cos 2φ + v2 sin 2φ.
(5)
In this expression, p v0 is the isotropic component of the velocity, and v12 + v22 the anisotropy. We assume the splitting time to be much smaller than the period of the wavefield. Because the wavefields of each polarization data are symmetric by 180 degrees, the anisotropy depends on polarization through a dependence of 2φ.
RETRIEVED NEAR-SURFACE SHEAR-WAVE VELOCITIES IN JAPAN
Using deconvolution interferometry at each station, we obtain the wavefield that corresponds to a plane wave propagating in the near-vertical direction (cos θ > 0.975) between the borehole receiver and surface receiver at each station. In this section, we show the wavefields of the annual stack, monthly stack, and shear-wave splitting.
Time (s)
Shear-wave splitting
We investigate shear-wave splitting by measuring the shear-wave velocity as a function of the polarization. We rotate both surface and borehole receivers from 0 to 350 degrees using a 10-degree interval. The north-south direction is denoted by 0 degrees, and the east-west direction by 90 degrees. Because a rotation over 180 degrees does not change the polarization, the 0- to 170-degree wavefields are the same as the 180- to 350-degree data. We apply deconvolution interferometry to the rotated wavefields, located at the surface and borehole receivers with the same polarizations, for determining the velocity of each polarization. Because the velocity for each polarization is related to the velocities of the fastest and slowest shear waves (Appendix B), we can estimate shear-wave splitting from the velocity difference. We cross-correlate the deconvolved wavefield for every used polarization (from 0 to 350 degrees in 10-degree intervals) with the deconvolved wavefield obtained from the motion in the north-south direction. This allows us to quantify the polarization dependence of the shear-wave velocity. Similar to the process described in section 4.1, we compute annual stacks of cross-correlated wavefields and pick the peak amplitudes of stacked wavefields by using quadratic interpolation. We can separate the velocity v(φ) as a function of polarization direction φ into the isotropic and anisotropic terms using a Fourier series expansion (Appendix C):
5
0.1
0.3 0.4 0.5 0.6
2000
2010
Figure 6. Annual-stacked wavefields (curves) with the interpolated largest amplitude (circles) at station NIGH13. The horizontal line at around 0.15 s is the shear-wave arrival time determined from logging data. From left to right, we show annual stacks from 2000 to 2010. The source and receiver polarization directions are the north-south direction.
5.1
Annual and monthly stacks
Figure 6 shows the annual-stacked wavefields at station NIGH13 represented by the large black circle in Figure 1. At this station, the sampling interval is 0.005 s until 2007 and is 0.01 s after 2008. In Figure 6, the deconvolved wavefields have good repeatability and a pronounced peak amplitude. After we apply quadratic interpolation (section 4.1), the determined arrival times (the black circles in Figure 6) correlate well with the travel time which is obtained from logging data (the horizontal line in Figure 6). The logging data is measured using a logging tool and by Vertical Seismic Profiling (VSP). The seismic source of VSP is a verticalcomponent vibrator. For finite offsets this source generates shear waves. We determine the average velocity from the logging data by computing the depth average of the slowness, because this quantity accounts for the vertical travel time. Because of the quadratic interpolation, the measured travel times show variations smaller than the original sampling time. After determining the arrival times of all stations, we compute the shear-wave velocities by using the known depth of the boreholes. Applying triangle-based cubic interpolation (Lawson, 1984) between stations, we create the shear-wave velocity map of Japan in each year (Figure 7). To reduce the uncertainty of the velocity estimation we use only the stations which give deconvolved waves with an arrival time greater than 0.1 s.
Near-surface S-wave velocities in Japan ◦
2 48 N 1 ◦
8E
◦E
2 13
◦E
6 13
◦E
0 14
◦E
◦
14
321
4 E 1 48
Log 44◦N
2000 2000
2005 36◦N
1500
2008 2010
1000
32◦N
Velocity (m /s)
40◦N
2002
500
0
Figure 7. Shear-wave velocities obtained from logging data (Log) and estimated by annual-stacked seismic interferometry using earthquake data at the north-south polarization (excerpted 2000, 2002, 2005, 2008, and 2010). The blue dots on a map represent the station locations which we use to make the map. We interpolate velocities between stations by triangle-based cubic interpolation (Lawson, 1984). The longitude and latitude belong to the map in the upper-left. The number of right-upper side of each map shows the year of data.
Thus, we obtain the near-surface shear-wave velocities throughout Japan by applying seismic interferometry to KiK-net data. The shear-wave velocity obtained from logging data is shown in the top left in Figure 7. Note that the velocities measured in different years are similar. In Figure 8, we crossplot the velocities estimated by interferometry in 2008 and obtained from logging data. The data are concentrated along the black line, which indicates the degree of correlation between the shearwave velocity obtained from logging data and from seismic interferometry. We also analyze seasonal changes and show the monthly-stacked wavefields at station NIGH13 in Figure 9. The monthly stacked wavefields also have good repeatability between different months.
5.2
Shear-wave splitting
In Figure 10a, we show the wavefields of the shear-wave splitting analysis at station NIGH13 in 2010 that are
obtained by the sequence of deconvolution and cross correlation described in section 4.4. Each trace is plotted at the angle that is equal to the shear-wave polarization used to compute that trace. The thick solid line in Figure 10a shows the interpolated maximum amplitude time of each waveform. The dashed circle shows the arrival time for the wave polarized in the northsouth direction. For the polarizations where the thick solid line is outside of the dashed circle, the shear-wave velocity is slower. The fast and slow shear polarization directions in Figure 10a are 22 degrees and −71 degrees clockwise from the north-south direction, respectively. The angle between these fast and slow directions is 93 degrees, which is close to 90 degrees as predicted by theory (Crampin (1985) and Appendix C). The 3-degree discrepancy could be caused by data noise or discretization errors. At station NIGH13 in 2010, the fast polarization shear-wave velocity vf ast is 638 m/s, the slow velocity vslow is 593 m/s, and the anisotropy parameter (vf ast − vslow )/vf ast is 7% (see Figure 10b). The dif-
322
N. Nakata & R. Snieder 0.0
1500
1000
0.2
Time (s)
Velocity (m/s) in 2008
0.1
500
0.3 0.4 0.5
0 0
500
1000
Velocity (m/s) from log
1500
Figure 8. Crossplot of velocities computed from logging data and by seismic interferometry in 2008. The black line indicates equal velocities.
ference of the arrival times between the fast- and slowpolarization velocity wavefields is much smaller than the period of a wavefield when the borehole depth is less than 525 m.
6
6.1
INTERPRETATION OF SHEAR-WAVE VELOCITIES AND SHEAR-WAVE SPLITTING Influence of major earthquakes
The near-surface shear-wave velocity in Japan is similar between years (see Figure 7), which means the nearsurface structure is basically stable. In this section, we focus on a small region. We use ∆v (calculated by the method presented in section 4.2) and the fast shear-wave polarizations shown in Figure 11 to analyze the influence of major earthquakes in the Niigata prefecture (the dark shaded area in Figure 1). Three significant earthquakes, shown by the dashed arrows in each panel, occurred during the 11 years. Figure 11a shows the velocity variation ∆v for the isotropic component v0 computed by equation 5 compiled over periods one year before and three months after the major earthquakes. We use all stations in the Niigata prefecture and compute the average over the stations. Each box depicts the time range (horizontal extent) and the error in the average velocity over that time interval (vertical length). The error in the velocity is given by the standard deviation of measurements from different earthquakes in each time interval (section 4.2). In Figure 11a, these average velocities show significant
0.6
Jan
Dec
Figure 9. Monthly-stacked wavefields (curves) with the interpolated largest amplitude (circles) at station NIGH13. The horizontal line at around 0.15 s is the shear-wave arrival time obtained from logging data. From left to right, we depict monthly stacks from January to December. Each trace is stacked over the 11 years (January 2000–2010, February 2000–2010, ..., and December 2000–2010). The source and receiver polarization directions are the north-south direction.
velocity reduction after the major earthquakes. The average isotropic velocity of all stations in the region from 2000 to 2010 is 662 m/s, and the relative velocity change of each earthquake is around 3-4%. Similar velocity variations caused by major earthquakes were reported earlier; for example, Sawazaki et al. (2009) analyze the variations caused by the 2000 Western-Tottori Earthquake, Yamada et al. (2010) analyze the variations caused by the 2008 Iwate-Miyagi Nairiku earthquake using KiKnet stations, while Nakata and Snieder (2011) observe a velocity reduction of about 5% after the 2011 TohokuOki earthquake. To increase the temporal resolution of the velocity change, we compute velocity changes averaged over periods one year before and three months after the major earthquakes (Figure 11a) because Sawazaki et al. (2009) found that the velocity reduction is sustained over a period of at least three months after an earthquake. The stations on soft-rock sites have a greater velocity reduction than those on hard-rock sites. (We define soft- and hard-rock sites from the estimated shearwave velocity; hard-rock sites have a shear-wave velocity greater than 600 m/s, while soft-rock sites have a shearwave velocity less than 600 m/s.) For the used eventstation pairs, the velocity reduction does not change measurably with the distance from the epicenter. This is an indication that the velocity reduction depends
Near-surface S-wave velocities in Japan
a)
323
b)
N
W
E
Velocity (m /s)
650
625
600 N
E Polarization direction
S
S
Figure 10. (a) Cross correlograms along every 10-degree polarization direction in 2010 at station NIGH13. Each trace is plotted at an angle equal to the polarization direction used to construct that trace. The dashed circle indicates the peak-amplitude time of the north-south direction, and the thick solid line represents the peak-amplitude time for each polarization direction. (b) Shear-wave velocities computed from the thick solid line in panel (a). Black circles represent the quadratic-interpolated fast and slow polarization shear-wave velocities.
b)
20 15
M6.8, Jul. 16, 2007
10
∆ v (m /s)
5
18 33
0 −5 −10
12 15
−15
29 262
−20 M5.3, Jan. 4, 2001 M6.8, Oct. 23, 2004 −25 2000
2002
2004
2006
2008
2010
Fast shearpolarization (degree)
a)
180 160 140 120 100 80 60 40
29 262 M6.8 Oct. 23, 2004
12 15 20 M5.3 Jan. 4, 2001 0 2000 2002
2004
18 33 M6.8 Jul. 16, 2007 2006 2008 2010
Figure 11. (a) The isotropic component of the shear-wave velocity (Fourier coefficient v0 ) averaged over one year before and three months after three major earthquakes. (b) The direction of the fast shear-wave polarization averaged over same intervals. All data are computed for stations in the Niigata prefecture (the dark shaded area in Figure 1). In each panel, the label of the year is placed in the middle of each year. The dashed arrows in panel (a) and the vertical dashed lines in panel (b) indicate the times of the three major earthquakes shown with the black crosses in Figure 1. The numbers at the left and right of dashed arrows (a) and lines (b) are the number of earthquakes we use for determining velocities and polarizations before and after the major earthquakes, respectively. In panel (a), each velocity is the velocity variation ∆v in equation 4. The horizontal extent of each box depicts the time interval used for averaging (one year before and three months after the major earthquakes). The vertical extent of each box represents the standard deviation of the velocity in the area computed by the method in section 4.2. The horizontal line in each box indicates average velocity ∆v in each time interval, and the vertical line the center of each time interval. In panel (b), we use only the stations with significant anisotropy ((vf ast − vslow )/vf ast ≥ 1%). We select 11 stations from 19 stations and depict the fast shear-wave polarization directions with different symbols or lines. Each symbol is placed at the center of each time period.
mostly on the local geology. The velocity reduction can be due either to the opening and closing of existing fractures, to the creation of new fractures, or to the change in the shear modulus caused by changes in the pore fluid
pressure because of shaking-induced compaction (Das, 1993, Figure 4.24). The relative velocity reduction is smaller than the reduction found by Wu et al. (2009) because of the aver-
N. Nakata & R. Snieder
∆ v (m/s)
a)
b)
25
25
20
20
15
15
10
10
5
5
∆ v (m/s)
324
0 −5
0 −5
−10
−10
−15
−15
−20
−20
−25
Jan
Apr
Jul
c)
−25 50
Oct
100
150 200 Precipitation (mm)
250
300
25 20 15
∆ v (m/s)
10 5 0 −5 −10 −15 −20 −25 50
100
150 200 Precipitation (mm)
250
300
Figure 12. Seasonal dependence of shear-wave velocity. (a) Variation of the average monthly-velocities stacked over the period 2000 through 2010 in southern Japan (the light gray area in Figure 1). We use the stations with the 15% slowest velocities in the area. The horizontal extent of each box shows time interval used for averaging, and the vertical extent the standard deviations of all receivers in the time interval computed by the method in section 4.2. The horizontal line in each box indicates average velocity ∆v in each time interval, and the vertical line the center of each time interval. (b) Crossplot between monthly precipitation (provide by JMA) and the average velocity ∆v with error bars. (c) Crossplot between monthly precipitation and the average velocity ∆v with error bars using the stations with the 85% fastest velocities in the area.
aging over stations and over earthquakes recorded over a period of three months. Wu et al. (2009) use a single station located on a soft-rock site and do not average over several months. Wegler et al. (2009) estimate the velocity reduction in deeper parts of the subsurface, and the velocity reduction they find is small (0.3-0.5 %). From this we infer that the velocity reduction due to a major earthquake is most pronounced in near surface, especially for soft-rock sites. We also obtain the polarization directions of the fast shear waves before and after the major earthquakes by averaging over the same time intervals as used in Figure 11a (Figure 11b). The direction of the fast shearwave polarization does not show a significant change after the earthquakes, hence it seems to be unaffected by these earthquakes. The average standard deviation of the polarization direction of all stations in the Niigata
prefecture between 2000 and 2010 is 15 degrees, which represents the accuracy of the fast shear-wave velocity polarization direction. 6.2
Influence of precipitation
We compute the monthly-averaged shear-wave velocities of the north-south polarization (Figure 12a) to investigate a possible seasonal velocity variation related to precipitation. We use only the data in southern Japan (the light shaded area in Figure 1) because that region has a more pronounced seasonal precipitation cycle than northern Japan. Figure 12a illustrates a significant velocity difference between spring/summer and fall/winter. We calculate the average velocities over the stations with the 15% slowest shear-wave velocities in the area because these stations are located at soft-rock
Near-surface S-wave velocities in Japan a)
325
40 oN
Fast shear polarization 39oN
Plate motion
38 oN
37oN
36 oN
b)
136oE
138oE
140oE
o 142 E
40 oN
Fast shear polarization o
39 N
Plate motion
38 oN Figure 14. Crossplot between the direction of the plate motion and the fast shear-wave polarization directions. We use the stations which have faster than 600 m/s shear-wave velocity. Red indicates there are many points. The north-south direction is 0 degrees, and the east-west direction 90 degrees.
37oN
36 oN
136oE
138oE
140oE
o 142 E
Figure 13. (a) Fast shear-wave polarization directions (black lines) and the direction of the plate motion (gray lines) estimated from GPS data (Sagiya et al., 2000) at the stations with significant anisotropy ((vf ast − vslow )/vf ast ≥ 1%). (b) Extracted stations from panel (a) with shear-wave velocity faster than 600 m/s.
sites and are therefore influenced more by precipitation than the station at hard-rock sites. We compare the monthly-averaged velocities with the monthly average of precipitation (observed by the Japan Meteorological Agency (JMA)) computed from precipitation records over 30 years (Figure 12b). Note the negative correlation between the shear-wave velocity and precipitation (i.e., when it rains, the velocity decreases), which is consistent with the findings of Sens-Sch¨ onfelder and Wegler (2006). For comparison, for the stations with the 85% fastest shear-wave velocities in the area (Figure 12c) the shear-wave velocity does not vary with precipitation. The cause of the velocity reduction is the decreased effective stress of the soil due to the infiltration of water that increases the pore pressure (Das, 1993, Section 4.19; Chapman and Godin, 2001; Snieder and van den Beukel, 2004). We assume that for soft-rock sites most of the velocity change is caused by the effective stress change because Snieder and van den Beukel (2004) show that the relative density change with pore pressure is
much smaller than the relative change in the shear modulus. 6.3
Shear-wave splitting and the direction of the plate motion
Using shear-wave splitting analysis, we determine fast shear-wave polarization directions of every station (illustrated by the black arrows in Figure 13a). These directions are averaged over all years from 2000 to 2010 because the temporal changes in the direction are small (see Figure 11b). We plot the directions of all stations which have an anisotropy parameter (vf ast − vslow )/vf ast ≥ 1% because the uncertainty in the direction of the fast shear polarization is large when the anisotropy is small. In Figure 13a, we also plot the direction of the plate motion at each station (the gray arrows), estimated from GPS data (Sagiya et al., 2000). Each arrow is normalized to the same length. In Figure 13b, we plot only the stations which have an anisotropy parameter larger than 1% and a northsouth polarization shear-wave velocity faster than 600 m/s; these stations are located on hard-rock sites. The average absolute angle between the directions of fast shear polarization and the plate motion in these stations is 16 degrees, and this average angle of the stations which have a shear-wave velocity less than 600 m/s is 36 degrees. Therefore, the fast shear-wave polarization on hard-rock sites correlates more strongly with the direc-
326
N. Nakata & R. Snieder
tion of the plate motion than the polarization on softrock sites. The 16-degree angle is close to the 15-degree standard deviation angle of each station computed in section 6.1. The near-surface polarization in the western part of Figure 13b correlates well with observations of the shear-wave polarization at greater depth (Okada et al., 1995; Nakajima and Hasegawa, 2004; Nakajima et al., 2006), but this agreement does not hold in the regions further east. We present a crossplot of the directions of the fast shear-wave polarization and the plate motion for stations all over Japan (Figure 14), where we only used the stations which have an anisotropy parameter greater than 1% and a shear-wave velocity faster than 600 m/s. The red area in Figure 14 indicates that for most stations the direction of the plate motion is between 90 and 140 degrees, and that this direction correlates with the polarization direction of the fast shear wave. The near-surface stress directions on hard-rock sites is presumably related to the plate motion because the stress field related to the plate motion changes the properties of fractures. Note that the used shear waves sample the shallow subsurface (down to about several hundreds of meters). It is remarkable that the shear-wave velocities in the near surface at hard-rock sites correlate with tectonic process (plate motion) that extends several tens of kilometers into the subsurface.
7
CONCLUSION
We obtain annual and monthly averaged near-surface shear-wave velocities throughout Japan by applying seismic interferometry to KiK-net data. Deconvolution interferometry yields more repeatable and higher resolution wavefields than does cross-correlation interferometry. Because picked arrival times in waveforms are generally stable over time and consistent with logging data, the near-surface has a stable shear-wave velocity. After three strong earthquakes in the Niigata prefecture, however, the shear-wave velocity is reduced. By computing the monthly-stacked velocity, we observe a velocity variation on stations placed on soft rock that has a negative correlation with precipitation. We also observe shear-wave splitting. The fast shear-wave polarization direction on a hard-rock site correlates with the direction of the plate motion. Because the shearwave velocity is related to ground soil strength, these velocities are useful for civil engineering, site characterization, and disaster prevention.
8
ACKNOWLEDGMENTS
We are grateful to NIED for providing us with the KiKnet data and to JMA for precipitation data. We thank Diane Witters for her professional help in preparing this
manuscript. We are grateful to Kaoru Sawazaki for suggestions, corrections, and discussions.
REFERENCES Aki, K., 1957, Space and time spectra of stationary stochastic waves, with special reference to microtremors.: Bull. of the Earthquake Res. Inst. Univ. Tokyo, 35, 415–456. Alford, R. M., 1986, Shear data in the presence of azimuthal anisotropy: Dilley, Texas: SEG Expanded Abstracts, 5, 476–479. Bakulin, A., and R. Calvert, 2004, Virtual source: new method for imaging and 4D below complex overburden: SEG Expanded Abstracts, 23, 112–115. ——–, 2006, The virtual source method: Theory and case study: Geophysics, 71, S1139–S1150. Bakulin, A., and A. Mateeva, 2008, Estimating interval shear-wave splitting from multicomponent virtual shear check shots: Geophysics, 73, A39–A43. Brenguier, F., M. Campillo, C. Hadziioannou, N. M. Shapiro, R. M. Nadeau, and E. Larose, 2008, Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations: Science, 321, 1478–1481. Chapman, D., and O. Godin, 2001, Dispersion of interface waves in sediments with power-law shear speed profiles. II. Experimental observations and seismoacoustic inversions: J. Acoust. Soc. Am., 110, 1908– 1916. Claerbout, J. F., 1968, Synthesis of a layered medium from its acoustic transmission response: Geophysics, 33, 264–269. Crampin, S., 1985, Evaluation of anisotropy by shearwave splitting: Geophysics, 50, 142–152. Das, B., 1993, Principles of soil dynamics: PWS-Kent Publishing Co. Draganov, D., X. Campman, J. Thorbecke, A. Verdel, and K. Wapenaar, 2009, Reflection images from ambient seismic noise: Geophysics, 74, A63–A67. Draganov, D., K. Wapenaar, W. Mulder, J. Singer, and A. Verdel, 2007, Retrieval of reflections from seismic background-noise measurements: Geophys. Res. Lett., 34, L04305. Hohl, D., and A. Mateeva, 2006, Passive seismic reflectivity imaging with ocean-bottom cable data: SEG Expanded Abstracts, 25, 1560–1564. Larose, E., L. Margerin, A. Derode, B. van Tiggelen, M. Campillo, N. Shapiro, A. Paul, L. stehly, and M. Ta, 2006, Correlation of random wavefields: An interdisciplinary review: Geophysics, 71, S111–S121. Lawson, C. L., 1984, C1 surface interpolation for scattered data on a sphere: Rocky Mountain J. Math., 14, 177–202. Li, Y.-G., J. E. Vidale, K. Aki, F. Xu, and T. Burdette, 1998, Evidence of shallow fault zone strengthening
Near-surface S-wave velocities in Japan after the 1992 M7.5 Landers, California, earthquake: Science, 279, 217–219. Lin, F.-C., M. H. Ritzwoller, and R. Snieder, 2009, Eikonal tomography: surface wave tomography by phase front tracking across a regional broad-band seismic array: Geophys. J. Int., 177, 1091–1110. Lobkis, O. I., and R. L. Weaver, 2001, On the emergence of the Green’s function in the correlations of a diffuse field: J. Acoust. Soc. Am., 110, 3011–3017. Ma, S., G. A. Prieto, and G. C. Beroza, 2008, Testing community velocity models for southern California using the ambient seismic field: Bull. Seismol. Soc. Am., 98, 2694–2714. Mehta, K., J. L. Sheiman, R. Snieder, and R. Calvert, 2008, Strengthening the virtual-source method for time-lapse monitoring: Geophysics, 73, S73–S80. Mehta, K., R. Snieder, and V. Graizer, 2007a, Downhole receiver function: a case study: Bull. Seismol. Soc. Am., 97, 1396–1403. ——–, 2007b, Extraction of near-surface properties for a lossy layered medium using the propagator matrix: Geophys. J. Int., 169, 271–280. Minato, S., T. Matsuoka, T. Tsuji, D. Draganov, J. Hunziker, and K. Wapenaar, 2011, Seismic interferometry using multidimensional deconvolution and crosscorrelation for crosswell seismic reflection data without borehole sources: Geophysics, 76, SA19– SA34. Miyazawa, M., R. Snieder, and A. Venkataraman, 2008, Application of seismic interferometry to exptract P- and S-wave propagation and observation of shear-wave splitting from noise data at Cold Lake, Alberta, Canada: Geophysics, 73, D35–D40. Nakajima, J., and A. Hasegawa, 2004, Shear-wave polarization anisotropy and subduction-induced flow in the mantle wedge of northeastern Japan: Earth Planet. Sci. Lett., 225, 365–377. Nakajima, J., T. Matsuzawa, A. Hasegawa, and D. Zhao, 2001, Three-dimensional structure of Vp , Vs , and Vp /Vs beneath northeastern Japan: Implications for arc magmatism and fluids: J. Geophys. Res., 106, 21843–21857. Nakajima, J., J. Shimizu, S. Hori, and A. Hasegawa, 2006, Shear-wave splitting beneath the southwestern Kurile arc and northeastern Japan arc: A new insight into mantle return flow: Geophys. Res. Lett., 33, L05305. Nakata, N., and R. Snieder, 2011, Near-surface weakening in Japan after the 2011 Tohoku-Oki earthquake: Geophys. Res. Lett., 38, L17302. Nakata, N., R. Snieder, T. Tsuji, K. Larner, and T. Matsuoka, 2011, Shear-wave imaging from traffic noise using seismic interferometry by cross-coherence: Geophysics, 76, SA97–SA106. Okada, T., T. Matsuzawa, and A. Hasegawa, 1995, Shear-wave polarization anisotropy beneath the north-eastern part of Honshu, Japan: Geophys. J.
327
Int., 123, 781–797. Okada, Y., K. Kasahara, S. Hori, K. Obara, S. Sekiguchi, H. Fujiwara, and A. Yamamoto, 2004, Recent progress of seismic observation networks in Japan — Hi-net, F-net, K-NET and KiK-net —: Earth Planets Space, 56, 15–28. Prieto, G. A., J. F. Lawrence, and G. C. Beroza, 2009, Anelastic earth structure from the coherency of the ambient seismic field: J. Geophys. Res., 114, B07303. Roux, P., and M. Fink, 2003, Green’s function estimation using secondary sources in a shallow wave environment: J. Acoust. Soc. Am., 113, 1406–1416. Ruigrok, E., X. Campman, D. Draganov, and K. Wapenaar, 2010, High-resolution lithospheric imaging with seismic interferometry: Geophys. J. Int., 183, 339– 357. Sagiya, T., S. Miyazaki, and T. Tada, 2000, Continuous GPS array and present-day crustal deformation of Japan: Pure appl. geophy., 157, 2303–2322. Sawazaki, K., H. Sato, H. Nakahara, and T. Nishimura, 2009, Time-lapse changes of seismic velocity in the shallow ground caused by strong ground motion shock of the 2000 western-Tottori earthquake, Japan, as revealed from coda deconvolution analysis: Bull. Seismol. Soc. Am., 99, 352–366. Schaff, D. P., and G. C. Beroza, 2004, Coseismic and postseismic velocity changes measured by repeating earthquakes: J. Geophys. Res., 109, B10302. Schuster, G. T., J. Yu, J. Sheng, and J. Rickett, 2004, Interferometric/daylight seismic imaging: Geophys. J. Int., 157, 838–852. Sens-Sch¨ onfelder, C., and U. Wegler, 2006, Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia: Geophys. Res. Lett., 33, L21302. Shiomi, K., K. Obara, S. Aoi, and K. Kasahara, 2003, Estimation on the azimuth of the Hi-net and KiKnet borehole seismometers (in japanese): Zishin2, 56, 99–110. Snieder, R., and E. S ¸ afak, 2006, Extracting the building response using seismic interferometry: Theory and application to the Millikan library in Pasadena, California: Bull. Seismol. Soc. Am., 96, 586–598. Snieder, R., and A. van den Beukel, 2004, The liquefaction cycle and the role of drainage in liquefaction: Granular Matter, 6, 1–9. Snieder, R., K. Wapenaar, and K. Larner, 2006, Spurious multiples in seismic interferometry of primaries: Geophysics, 71, SI111–SI124. Stehly, L., M. Campillo, B. Froment, and R. L. Weaver, 2008, Reconstruction Green’s function by correlation of the coda of the correlation (C3 ) of ambient seismic noise: J. Geophys. Res., 113, B11306. Thomsen, L., 1988, Reflection seismology over azimuthally anisotropic media: Geophysics, 53, 304– 313. Tonegawa, T., K. Nishida, T. Watanabe, and K. Sh-
328
N. Nakata & R. Snieder
iomi, 2009, Seismic interferometry of teleseismic Swave coda retrieval of body waves: an application to the Philippine Sea slab underneath the Japanese Islands: Geophys. J. Int., 178, 1574–1586. Trampert, J., M. Cara, and M. Frogneux, 1993, SH propagator matrix and Qs estimates from boreholeand surface- recorded earthquake data: Geophys. J. Int., 112, 290–299. Vasconcelos, I., and R. Snieder, 2008, Interferometry by deconvolution: Part 2 - Theory for elastic waves and application to drill-bit seismic imaging: Geophysics, 73, S129–S141. Vidale, J. E., and Y.-G. Li, 2003, Damage to the shallow Landers fault from the nearby Hector Mine earthquake: Nature, 421, 524–526. Wapenaar, K., 2004, Retrieving the elastodynamic Green’s function of an arbitrary inhomogeneous medium by cross correlation: Phys. Rev. Lett., 93, 254–301. Wapenaar, K., D. Draganov, R. Snieder, X. Campman, and A. Verdel, 2010, Tutorial on seismic interferometry: Part 1 - Basic principles and applications: Geophysics, 75, 75A195–75A209. Wapenaar, K., and J. Fokkema, 2006, Green’s function representations for seismic interferometry: Geophysics, 71, SI33–SI46. Wapenaar, K., J. van der Neut, and E. Ruigrok, 2008, Passive seismic interferometry by multidimensional deconvolution: Geophysics, 73, A51–A56. Wegler, U., H. Nakahara, C. Sens-Sch¨ onfelder, M. Korn, and K. Shiomi, 2009, Sudden drop of seismic velocity after the 2004 Mw 6.6 mid-Niigata earthquake, Japan, observed with passive image interferometry: J. Geophys. Res., 114, B06305. Wegler, U., and C. Sens-Sch¨ onfelder, 2007, Fault zone monitoring with passive image interferometry: Geophys. J. Int., 168, 1029–1033. Wu, C., Z. Peng, and D. Assimaki, 2009, Temporal changes in site response associated with the strong ground motion of the 2004 MW 6.6 Mid-Niigata earthquake sequences in Japan: Bull. Seismol. Soc. Am., 99, 3487–3495. Yamada, M., J. Mori, and S. Ohmi, 2010, Temporal changes of subsurface velocities during strong shaking as seen from seismic interferometry: J. Geophys. Res., 115, B03302.
APPENDIX A: 1D SEISMIC INTERFEROMETRY We explain in this appendix why deconvolution interferometry is suitable for this study. Figure A1 illustrates the model of interferometry using a KiK-net station and an earthquake. The incoming wavefield S(rb , s, ω), propagating from source s to receiver rb , is given by G(rb , s, ω)W (s, ω), where G is the Green’s function including any unknown complex effect of wave propaga-
z=0
rs
z=z b
rb
s Figure A1. Geometry of an earthquake and a KiK-net station, where rs is the surface receiver (black triangle), rb the borehole receiver (white triangle), and s the epicenter of earthquake (gray star).
tion such as scattering and attenuation, and W the source signature in the frequency domain. Assuming that the subsurface is homogeneous between the receivers, the received wavefield at surface receiver rs is 2S(rb , s, ω)ejkzb e−γzb , where γ is the attenuation coefficient and k the wave number. Because of the free surface, the amplitude of the wavefield at the surface is multiplied by a factor 2. We assume that there are no multiples between the two receivers. The reflected wavefield from the surface at the borehole receiver rb is S(rb , s, ω)e2jkzb e−2γzb , and the total wavefield at rb is S(rb , s, ω) + S(rb , s, ω)e2jkzb e−2γzb . Applying crosscorrelation interferometry to these wavefields yields
u(rs , s, ω)u∗ (rb , s, ω) = 2S(rb , s, ω)ejkzb e−γzb h i S ∗ (rb , s, ω) + S ∗ (rb , s, ω)e−2jkzb e−2γzb ≈ 2|S(rb , s, ω)|2 ejkzb e−γzb
(A1)
when we consider only the first arrival. Using deconvolution interferometry, we obtain
Near-surface S-wave velocities in Japan [h]
in the horizontal direction. In our study, the virtual source in the borehole has the polarization of the incident wave. The two horizontal components of the virtual source therefore are not independent, so that Alford rotation cannot be applied. The angle of the fast and slow shear-wave polarization directions is 90 degrees because we assume the incoming wave is a plane wave (Crampin, 1985). A waveˆ which is a unit vector, can be field with polarization p, expressed in the polarization of the fast and slow shear waves:
pˆ s
pˆ
pˆ s sin φ
pˆ = pˆf cos φ + pˆs sin φ,
φ pˆ f cos φ
pˆ f
Figure B1. Projection of fast and slow velocity directions, ˆf is the fast polarization direction, p ˆs the slow polarwhere p ˆ an arbitrary direction, and φ the angle ization direction, p ˆ p ˆf , between the fast direction and arbitrary direction. p, ˆs are unit vectors. Dashed arrows show the projection, and p which is shown in equation B1.
jkzb −γzb
e 2S(rb , s, ω)e u(rs , s, ω) = u(rb , s, ω) S(rb , s, ω) + S(rb , s, ω)e2jkzb e−2γzb ≈ 2ejkzb e−γzb ,
329
(A2)
where we also only retain the first arrival. The planewave Green’s function excited at rb and received at rs is equal to 2ejkzb e−γzb in the frequency domain. Let us compare equations A1, A2, and the Green’s function. The wavefield retrieved by cross-correlation interferometry is complicated because equation A1 includes the power spectrum |S(rb , s, ω)|2 . This term is different for different earthquakes. In contrast, deconvolution interferometry eliminates the incoming wave S(rb , s, ω), and thus provides a more accurate estimate of the Green’s function. When we stack deconvolved wavefields over earthquakes, the accuracy of this estimate is improved. Because the deconvolved waves do not depend on the power spectrum of the incident wave |S(rb , s, ω)|2 , the deconvolved wavefields are more reproducible than those obtained from cross correlation.
APPENDIX B: SHEAR-WAVE SPLITTING Usually, shear-wave splitting is analyzed with Alford rotation (Alford, 1986; Thomsen, 1988). This procedure is based on the use of two independent orthogonal sources
(B1)
where φ is the polarization angle in the arbitrary wavefield relative to the direction of the fast shear-wave polarization, and pˆf and pˆs are the unit vectors of the fast and slow velocity wavefields, respectively (see Figure B1). The incoming wavefield ub at the borehole receiver is
ub = S(t)pˆ = S(t)pˆf cos φ + S(t)pˆs sin φ,
(B2)
where S(t) is the incoming wavefield, and t the time. The wavefield at surface receiver us is given by „ « „ « zb zb us = S t − pˆf cos φ + S t − pˆs sin φ, (B3) vf vs where vf is the fast velocity, vs the slow velocity, and zb the distance between the top and bottom receivers. The component of us along pˆ is
us (φ) =(pˆ · us ) „ « zb =S t − (pˆ · pˆf ) cos φ vf „ « zb +S t− (pˆ · pˆs ) sin φ vs « „ « „ zb zb cos2 φ + S t − =S t − sin2 φ. vf vs
(B4)
We express vf and vs in the average velocity v0 and the difference δ: 1 1 = − δ, vf v0
1 1 = + δ, vs v0
(B5)
and assume the splitting time zb δ is much smaller than the period. We insert expression B5 into equation B4 and expand using first-order Taylor expansion in δ:
330
N. Nakata & R. Snieder „
«
zb us (φ) =S t − + zb δ cos2 φ v0 „ « zb − zb δ sin2 φ +S t− v0 „ « zb =S t − (sin2 φ + cos2 φ) v0 „ « zb 0 +S t− zb δ(cos2 φ − sin2 φ) v0 „ « „ « zb zb =S t − + S0 t − zb δ cos 2φ v0 v0 « „ zb + zb δ cos 2φ , (B6) =S t − v0 where S 0 is the time derivative of S. Thus, using Taylor expansion, we obtain the velocity for a shear wave with the polarization of equation B1: v(φ) =
v0 , 1 − v0 δ cos 2φ
(B7)
or to first order in v0 δ: v(φ) = v0 (1 + v0 δ cos 2φ).
(B8)
APPENDIX C: FOURIER COEFFICIENTS We describe the meaning of v0 , v1 , and v2 in equation 5. Expression B8 gives the velocity for a polarization φ relative to the polarization of the fast shear wave. When the azimuth of the fast shear-wave polarization is given by ψ, the angle φ in equation B8 must be changed into φ → φ − ψ. Denoting v02 δ by V , these changes turn equation B8 into v(φ) = v0 + V cos 2(φ − ψ).
(C1)
This can also be written as v(φ) = v0 + v1 cos 2φ + v2 sin 2φ;
(C2)
with
Z π Z 1 1 π v(φ)dφ, v1 = v(φ) cos 2φdφ, 2π −π π −π Z π 1 v2 = v(φ) sin 2φdφ, π −π „ « q 1 v2 V = v12 + v22 , ψ = arctan . 2 v1 v0 =
In expression C2, v0 , v1 , and v2 are the Fourier coefficients of the velocity v(φ). Because v0 does not depend on φ, v0 represents the isotropic velocity. Expression C1 shows that v(ψ) and v(ψ + π/2) are the fastest and slowest velocities, respectively. Therefore, V is the anisotropic velocity and the angle between the fast and slow polarization directions is 90 degrees,
which corresponds to shear-wave splitting of a plane wave (Crampin, 1985).
