## Abstract

We analysed background surface waves in seismic ambient noise by cross-correlating continuous records of eight ocean bottom seismometers and nine differential pressure gauges deployed in the northwestern Pacific Ocean by the PLATE project. After estimating the clock delay and instrumental phase responses of differential pressure gauges by using cross-correlation functions, we measured average phase velocities in the area of the array for the fundamental-, first higher- and second higher-mode Rayleigh waves, and the fundamental-mode Love waves at a period range of 3–40 s by waveform fitting. We then measured phase-velocity anomalies of fundamental-mode and first higher-mode Rayleigh waves for each pair of stations at a period range of 5–25 s, and corrected the effect of variation in water-depths. The seismic anomalies imply the presence of strong azimuthal anisotropy beneath the eastern part of array. The direction of maximum velocity is approximately N35°E in the fossil seafloor spreading direction perpendicular to magnetic lineations from the ancient triple junction at this location. The peak-to-peak intensity of shear-wave velocity anisotropy in the mantle is ∼7 per cent.

## INTRODUCTION

We study seismic anisotropy in the oceanic lithosphere beneath old seafloor in the NW Pacific using data collected from a marine seismic deployment, the PLATE Project (Pacific Lithosphere Anisotropy and Thickness Experiment). The region has the oldest identifiable magnetic anomalies in the Pacific, ∼150–160 Ma, and is in an area not affected by subsequent Cretaceous intraplate volcanism (Nakanishi et al. 19891992). We exploit the crystalographic properties of olivine which display lattice preferred orientation (LPO) in response to shear in the uppermost mantle.

Seismic anisotropy beneath oceanic regions has been estimated for understanding deformation in the mantle related to plate tectonics from the 1960s. For example, azimuthal anisotropy of Pn waves obtained by refraction surveys indicates that the mantle flow at the top of the mantle at depths of ∼10–15 km was perpendicular to the ancient mid-ocean ridge (Francis 1969; Raitt et al. 1969).

On the other hand, azimuthal anisotropy of S waves estimated by surface wave tomography at periods of ∼50–200 s indicates that the mantle flow in the oceanic asthenosphere at depths of ∼50–200 km is parallel to current plate motion (Tanimoto and Anderson 1984). For investigating the intermediate depth range, recent studies have analysed surface waves at shorter periods of ∼18–50 s either by measuring group velocities from onland records (Smith et al. 2004) or by measuring phase velocities from records of ocean bottom seismometers (OBSs; Forsyth & Li 2005; Weeraratne et al. 2007).

Another method to analyse surface waves at even shorter periods of ∼1–50 s is to extract background surface wave propagation in ambient noise by cross-correlating continuous records of seismometers (Shapiro and Campillo 2004). As each cross-correlation function (CCF) is mainly sensitive to the structure along the path between the pair of stations (Snieder 2004; Tromp et al. 2010; Nishida 2011), regional tomography is possible to determine S-wave velocity structure of the crust and the uppermost mantle (e.g. Bensen et al. 2008; Nishida et al. 2008) including azimuthal anisotropy (Lin et al. 2010). Harmon et al. (2007) first applied this method to OBS records for measuring phase and group velocities of Rayleigh waves at periods of 2–16 s in the East Pacific Rise (EPR) region. From the phase velocities of first higher-mode Rayleigh waves at periods of 3.5–7 s, they inferred the possibility of azimuthal anisotropy at depths of ∼10–15 km in the EPR region. At shorter periods of 0.5–3 s, Mordret et al. (2013b) analysed fundamental-mode Rayleigh waves, and estimated horizontal distribution of azimuthal anisotropy within the sediment layer in the North Sea.

In this study, we investigate seismic anisotropy in the uppermost mantle at depths of ∼10–50 km by applying the cross-correlation method to records from the PLATE project. The PLATE project was originally designed with 16 OBSs (each with three component seismometers and a differential pressure gauge; DPG) which were deployed in two arrays of eight OBS, which we refer to as the western and eastern arrays (Fig. 1). The two arrays were chosen to compare structures for different azimuths of isochrones, that is, orthogonal directions of seafloor spreading. While this location was specifically chosen for its old seafloor age, the location required deployment at great water depths (∼6000 m) and high pressures that lead to failure of several instruments. The final array of working instruments is indicated in Fig. 1 and was reduced to eight OBSs and nine DPGs, which were not necessarily coincident at any one station.

Figure 1.

Station map of the PLATE project in the lower map. Geographical location of the experiment is shown by the rectangle in the top map. Triangles and crosses show three-component OBSs. Circles show DPGs. The horizontal-component records of OBSs shown by crosses and the DPGs shown by broken circles are not used in this study due to low data quality. White lines show the isochrones based on the magnetic lineations. SR indicates the Shatsky Rise located northeast of the stations.

Figure 1.

Station map of the PLATE project in the lower map. Geographical location of the experiment is shown by the rectangle in the top map. Triangles and crosses show three-component OBSs. Circles show DPGs. The horizontal-component records of OBSs shown by crosses and the DPGs shown by broken circles are not used in this study due to low data quality. White lines show the isochrones based on the magnetic lineations. SR indicates the Shatsky Rise located northeast of the stations.

In this study, we first overcome three problems in our analysis: unknown clock delays, imperfectly known instrumental responses of DPGs, and the overlap between multimode surface waves. Although these problems have not all been addressed in previous studies (e.g. Harmon et al. 2007; Mordret et al. 2013b), the problems may exist in any other records especially collected on the seafloor.

The first problem is the delay of the clock. As the GPS signal is unavailable in the ocean, seafloor instruments are equipped with inner high-quality crystal clocks, whose drifts are typically one second per year or less and can be corrected after the recovery by comparison to GPS signals. If the clock stops during the observation period or if the clock drift is non-linear, however, the delay of the clock needs to be corrected by other methods.

The second problem is the accuracy of the instrumental response when DPGs are used instead of OBSs. The DPG is an instrument that records the difference between seafloor pressure and low-pass filtered pressure, and gives relative pressure change at periods shorter than about 1000 s with some phase delay (Cox et al. 1984). The instrumental responses of DPGs are usually imperfectly known because they depend on the viscosity of the fluid and the diameter (and the cleanness) of the capillary tube (Cox et al. 1984).

The third problem is the overlap between multimode surface waves when the interstation distances are comparable to the wavelength. For example, higher-mode Rayleigh waves are observable and important in oceanic regions at periods shorter than ∼10–20 s, where the fundamental-mode Rayleigh wave is mainly sensitive to P-wave velocity in the ocean and higher-mode Rayleigh waves are sensitive to S-wave velocity in the crust and uppermost mantle (Ewing et al. 1957). Previous studies analysed each mode of Rayleigh wave separately by applying group velocity filters (Harmon et al. 2007; Yao et al. 2011; Takeo et al. 2013). The separation, however, becomes difficult for shorter interstation distances due to overlap between different modes. The separation of Rayleigh and Love waves in horizontal components is also difficult for short interstation distances or at long periods, where the Rayleigh and Love waves appear in the both radial and transverse components (Aki 1957). Although both Rayleigh and Love waves have been analysed in land (Bensen et al. 2008; Nishida et al. 2008) and oceanic regions (Mordret et al. 2013a; Takeo et al. 2013), this effect has not been accounted for.

We therefore first used CCFs for estimating clock delays and the instrumental phase responses of DPGs. We then measured phase velocities of multimode Rayleigh waves and the fundamental-mode Love wave simultaneously by developing a waveform fitting method. The obtained average phase velocities are used for constraining 1-D isotropic structure in the crust and uppermost mantle beneath the array. The phase velocity anomalies of Rayleigh waves are finally used for determining seismic anisotropy in the uppermost mantle.

## DATA

One-year records of eight OBSs and nine DPGs were acquired at 11 stations in the northwestern Pacific Basin by the PLATE project in 2009–2010. Fig. 1 shows the seafloor topography and location of stations. We use five stations (stations 1–8) in the western array, and six stations (stations 10–16) in the eastern array. The interstation distances range from 60 to 600 km. There are three types of OBSs as summarized in Table 1: OBSs with a velocity-flat response at periods shorter than 240 s (Trillium 240; T-240) or shorter than 40 s (Trillium 40; T-40), and OBSs with an acceleration-flat response at periods longer than 1 s developed in Lamont-Doherty Earth Observatory (L-DEO; Webb et al. 2001). If the pressure change is only caused by seismic surface waves, the phase of the pressure is theoretically the same as that of the vertical displacement. We therefore regarded the DPG records as displacement waveforms, and estimated the unknown instrumental phase responses of DPGs from the CCFs between DPG and OBS records in Section 3.2. The instrumental amplitude responses of DPGs were roughly estimated by comparing the amplitudes of teleseismic waveforms of DPGs and OBSs at a period range of 50–100 s. Although the amplitude response may depend on the period of analysis, the amplitude responses estimated at periods of 50–100 s were simply extrapolated to shorter periods. As the sampling rate is 50 Hz for eight stations and 40 Hz for three stations and we are interested in surface wave dispersion at periods greater than 1 s, we first down-sampled all records to 10 Hz.

Table 1.

Types of OBSs.

Station OBS
T-240
T-40
–
T-240
LDEO
10 –
11 T-240
12 LDEO
13 –
15 T-40
16 T-240
Station OBS
T-240
T-40
–
T-240
LDEO
10 –
11 T-240
12 LDEO
13 –
15 T-40
16 T-240

## ANALYSES AND RESULTS

### Calculation of CCFs

We first calculated CCFs between every pair of stations and components by calculating the cross spectra in the frequency domain and by taking the inverse Fourier transform. The method is the same as Takeo et al. (2013). The cross-spectrum is defined as,

(1)
$$S_{ij}^{kl} = \left \langle W_{ij}^{kl} \cdot F_i^k \cdot \left( {F_j^l } \right)^* \right\rangle /\left\langle W_{ij}^{kl} \right\rangle ,$$
where $$F_i^k$$ is a Fourier spectrum of the kth component record at the ith station with a time window of 819.2 s. The instrumental responses were corrected. The weighting term, $$W_{ij}^{kl} = 1/( {| {f_i^k } | \cdot | {f_j^l } |} )$$, is introduced for reducing the contribution of records with high noise level.

Fig. 2 shows the CCFs between DPGs (k = l = P, PP component), DPG and vertical component of OBS (PZ), vertical (ZZ), radial (RR) and transverse (TT) components of OBSs for three bandpass filters: 3–5 s, 5–10 and 10–20 s. The numbers of CCFs are summarized in Table 2. As three stations only have DPG records and two stations only have OBS records, we use the PZ component to increase the number of available pairs. The TT-component CCFs show the propagation of fundamental-mode Love waves (0T mode), whereas the PP-, PZ-, ZZ- and RR- component CCFs show the propagation of multimode Rayleigh waves: fundamental mode (0S mode), first higher mode (1S mode) and second higher mode (2S mode). The Rayleigh waves with a velocity lower than ∼2 km s−1 mainly reflect P-wave velocity in the ocean, ∼1.5 km s−1, whereas the waves with a velocity higher than ∼2 km s−1 mainly reflect S-wave velocity in the solid (Ewing et al. 1957). These waves are called generalized Rayleigh waves (Ewing et al. 1957) or Scholte-Rayleigh waves (Yao et al. 2011). We hereafter refer to them as Rayleigh waves for simplicity. In addition, we call the slow waves ocean-mode Rayleigh waves and the fast waves solid-mode Rayleigh waves. For example, 0S mode transitions from ocean-mode Rayleigh wave to solid-mode Rayleigh wave at periods of ∼10–20 s. At periods of 3–5 s in this water depth, the 1S mode arrives later than the 0S mode due to the sudden decrease of the phase velocity of 1S mode (Ewing et al. 1957; the group velocity governing energy arrival depends on the derivative of the phase velocity).

Figure 2.

CCFs bandpass filtered at 10–40 s, 5–10 and 3–5 s. Dashed pink lines show the arrival of each mode at a typical period. Arrow in PZ component shows unknown phase.

Figure 2.

CCFs bandpass filtered at 10–40 s, 5–10 and 3–5 s. Dashed pink lines show the arrival of each mode at a typical period. Arrow in PZ component shows unknown phase.

Table 2.

Number of CCFs for each component.

Component PP PZ ZZ RR TT
Number of CCFs shown in Fig. 2 36 66 28 10 10
Number of CCFs used for analysis 16 30 12 10 10
Component PP PZ ZZ RR TT
Number of CCFs shown in Fig. 2 36 66 28 10 10
Number of CCFs used for analysis 16 30 12 10 10

The PP and PZ components mainly show the ocean-mode Rayleigh waves, whereas the RR component mainly shows the solid-mode Rayleigh waves. The ZZ component shows both modes simultaneously. The difference between appearances can be understood by considering the eigenfunctions of Rayleigh waves at the seafloor, which are continuous for the vertical displacement but are discontinuous for the horizontal displacement. As a result, the ocean-mode Rayleigh waves appear in both the DPGs and vertical components of OBSs, but the amplitude is smaller in the horizontal components of OBSs.

Fig. 2 further shows an unknown phase marked by a black arrow in the positive lag-time of the PZ component at periods of 5–10 s and interstation distances longer than 300 km. The apparent group velocity at the interstation distance range of 300–700 km is ∼3.5 km s−1, which is close to the group velocity of 1S mode at this period range, 3.6 km s−1. The simplest interpretation is that the signal corresponds to the conversion of 0S mode to 1S mode at seamounts located between two arrays in our study or at the Shatsky Rise at northeast of our arrays (see Fig. 1). As we pointed out already, the ocean-mode Rayleigh waves tend to appear in DPGs, whereas the solid-mode Rayleigh waves tend to appear in the vertical components of OBSs. At periods of 5–10 s, 0S and 1S modes appear in DPGs and vertical components of OBSs, respectively. The converted waves are therefore expected to appear in the CCFs between the DPGs and the vertical components of the OBSs.

### Correction for clock delays and instrumental responses

The accuracy of traveltime anomalies is critical for the estimation of azimuthal anisotropy. For the data set used in this study, the typical travel time is about 30–100 s. The traveltime accuracy of about 0.3 s is then needed to achieve phase-velocity accuracy of 1 per cent and to discuss the azimuthal dependence of phase-velocity anomalies. We therefore corrected the clock delays, instrumental responses and constant delays associated with the recording system in three steps by using CCFs.

In the first step, we corrected the clock delays at station 4 and station 12. At station 4, the total delay of the clock at the end of the project was ∼20 s, an order of magnitude larger than other stations, implying the possibility of significant non-linear drift or a clock jump. At station 12, the drift correction was not possible due to the cessation of recording during the observation period. We therefore estimated the clock delays at these two stations from the temporal changes of CCFs as done by Sens-Schönfelder (2008) for land seismic records without GPS reception and recently by Hannemann et al. (2013) and Gouedard et al. (2014) for oceanic regions. Fig. 3(a) shows the CCFs calculated from 5-d records and the CCFs averaged over the whole year with a bandpass filter of 3–10 s. Although the waveforms seem identical, the lag time of the signal slightly varies with day. Fig. 3(b) shows the time shifts, Δτ, that maximize cross-correlation coefficients between 5-d CCFs and the whole-year CCFs for PP, PZ and ZZ components with a time-window from −300 to 300 s. For pairs of stations in the western array, the time shift varies with day for pairs including station 4, but there is no discernible variation from zero for other pairs of stations (Fig. 3b). For pairs of stations in the eastern array, the time shift varies with day for pairs including station 12, but is zero for other pairs of stations (Fig. 3c). We fit the temporal change by a linear function of time for station 12. On the other hand, the temporal change for station 4 is not linear. We fitted the change by an empirical summation of a linear function and an exponential function with time. A low-order polynomial function works equally well. The CCFs between all the pairs of stations are then re-calculated by assuming the clock delay at the start of the observation period to be zero.

Figure 3.

(a) CCFs calculated from 5-d records (black) and the whole-year records (red) as a function of day from the start of observation. (b, c) Time-shifts between 5-d and the whole-year CCFs as a function of the day for pairs of stations in western and eastern arrays, respectively. Red lines show the fitting lines.

Figure 3.

(a) CCFs calculated from 5-d records (black) and the whole-year records (red) as a function of day from the start of observation. (b, c) Time-shifts between 5-d and the whole-year CCFs as a function of the day for pairs of stations in western and eastern arrays, respectively. Red lines show the fitting lines.

The apparent time shift due to inhomogeneous source distribution cannot be neglected especially when the interstation distances are comparable to the wavelength (Yao & van der Hilst 2009). The absence of a temporal change for pairs without station 4 or 12, however, suggests that the effect of seasonal variation of source distribution to the estimation of clock delay is small. If the temporal changes of time shifts for the pairs including station 4 or 12 are caused by seasonal variation of source distribution, the same pattern should appear in the time shifts for the pairs without station 4 or 12.

In the second step, we estimated the response of each DPG. In a previous study, Araki and Sugioka (2009) estimated the response by comparing the DPG record with a record of a quartz pressure gauge. In this study, we estimate the response from CCFs between the DPG record and the vertical-component record of the OBS deployed in the same station. If both DPG and OBS records only contain background surface-waves, the phase of the CCF, Δϕ, reflects the instrumental response of DPG. Fig. 4(a) shows the phase for all pairs of stations. At frequency of 0.1–0.3 Hz, the curves are stable and have slight offsets from zero, which seem to reflect the instrumental responses of the DPGs. On the other hand, the curves are unstable at frequencies higher than 0.3 Hz, where the wavelength is shorter than the water-depth and reverberation of body waves between the seafloor and the sea surface is possible. We therefore fitted the phase curves at frequencies of 0.1–0.3 Hz for obtaining the instrument response of DPGs. Fig. 4 shows the fitting line and the constant phase delay obtained by cross-correlating the DPG and the vertical records of teleseismic Rayleigh waves at the frequency range of 0.016–0.05 Hz (periods of 20–60 s). For the DPG stations without any working OBS (6, 10 and 13), we assume that the response is given by the average of the responses for four stations (7, 11, 15 and 16), which give similar responses (Fig. 4b).

Figure 4.

(a) Phase difference between DPG and vertical component of OBS as a function of frequency. Red solid lines show the linear fitting line, whereas blue dashed lines show the constant value obtained from teleseismic waveforms at periods of 20–60 s. (b) Same as (a) but for the average of four stations (7, 11, 15 and 16), which is used for DPG stations without OBS records (station 6, 10 and 13).

Figure 4.

(a) Phase difference between DPG and vertical component of OBS as a function of frequency. Red solid lines show the linear fitting line, whereas blue dashed lines show the constant value obtained from teleseismic waveforms at periods of 20–60 s. (b) Same as (a) but for the average of four stations (7, 11, 15 and 16), which is used for DPG stations without OBS records (station 6, 10 and 13).

In the third step, we evaluated the effectiveness of the response correction from the time symmetry of CCFs between different stations. We applied bandpass filters with a width of 0.05 Hz for each CCFs, and measured the time shifts that maximize the cross-correlation coefficients between positive and negative lag time of each CCFs. The time shift becomes zero if the CCF is time symmetric and becomes 2Δτ′ if the clock of one station delays by Δτ′ or the phase of one station delays 2πΔτ′/T, where T is period. The half amplitudes of time shifts are non-zero before the correction of DPG responses (Fig. 5a), but become smaller than 0.2 s after the correction of DPG responses even for the DPG station without a corresponding OBS (Fig. 5b). This result indicates that the assumed response for DPG without an OBS is valid with the uncertainty smaller than 0.2 s. Time shifts of about 0.5 s, however, still exist for the pairs including stations 8 and 12 (red lines in Fig. 5b). We concluded that these time shifts are due to recording systems such as the FIR filter because these two stations were operated by different institutions recording at different sampling rates on different recording systems, and corrected the timing of these stations by 0.5 s. Fig. 5(c) shows the time shifts after the DPG responses and the timing of 0.5 s were corrected. The amplitudes of time shifts are almost all smaller than 0.2 s, and indicate that the phase and timing of records are corrected with uncertainties smaller than 0.2 s.

Figure 5.

Half the time-shifts needed to obtain time-symmetric CCFs for each PP, PZ and ZZ component CCFs as a function of frequency (a) after the correction of clock delays with time and before the correction of DPG responses, (b) after the correction of DPG responses and (c) after the corrections of DPG responses and the timing delay of 0.5 s at station 8 and 12. Red lines in (b) show the pairs including station 8 or 12.

Figure 5.

Half the time-shifts needed to obtain time-symmetric CCFs for each PP, PZ and ZZ component CCFs as a function of frequency (a) after the correction of clock delays with time and before the correction of DPG responses, (b) after the correction of DPG responses and (c) after the corrections of DPG responses and the timing delay of 0.5 s at station 8 and 12. Red lines in (b) show the pairs including station 8 or 12.

### Measurement of average phase velocities

Four modes of surface waves can be recognized from CCFs: 0S, 1S, 2S and 0T modes. For obtaining the average phase velocities beneath the array, we measured the phase velocities and amplitudes of all modes simultaneously by assuming a constant phase velocity and amplitude at each mode for all pairs of stations, and by fitting the synthetic CCFs to the observed CCFs in the frequency domain using the notation of spatial autocorrelation (SPAC) method by Aki (1957). Similar multimode methods have already been applied to onland records at a frequency higher than 1 Hz (e.g. Beaty et al. 2002), but have not been applied to broadband records to investigate mantle yet.

The model parameters are phase velocities (cmode) and amplitudes of the four modes in five components at each frequency (f = ω/2π). According to the equation by Aki (1957), the synthetic cross-spectra, $${S\!^{^\prime}}_{ij}^{kl}$$, are given by,

(2)
$${{\it {S\!^{^\prime}}}}_{ij}^{PP} = a_{0S}^{PP} J_0 \left(\frac{\omega d_{ij}}{{c_{0S}}} \right) + a_{1S}^{PP} J_0 \left(\frac{\omega d_{ij}}{c_{1S}}\right),$$

(3)
$${{\it {S\!^{^\prime}}}}_{ij}^{PZ} = a_{0S}^{PZ} J_0 \left( {\frac{{\omega d_{ij} }}{{c_{0S} }}} \right) + a_{1S}^{PZ} J_0 \left( {\frac{{\omega d_{ij} }}{{c_{1S} }}} \right),$$

(4)
$${{\it {S\!^{^\prime}}}}_{ij}^{ZZ} = a_{0S}^{ZZ} J_0 \left( {\frac{{\omega d_{ij} }}{{c_{0S} }}} \right) + a_{1S}^{ZZ} J_0 \left( {\frac{{\omega d_{ij} }}{{c_{1S} }}} \right) + a_{2S}^{ZZ} J_0 \left( {\frac{{\omega d_{ij} }}{{c_{2S} }}} \right),$$

(5)
\begin{eqnarray} {{\it {S\!^{^\prime}}}}_{ij}^{RR} &=& a_{0S}^{RR} J_{0 - 2} \left( {\frac{{\omega d_{ij} }}{{c_{0S} }}} \right) + a_{1S}^{RR} J_{0 - 2} \left( {\frac{{\omega d_{ij} }}{{c_{1S} }}} \right) \nonumber\\ &&+\, a_{2S}^{RR} J_{0 - 2} \left( {\frac{{\omega d_{ij} }}{{c_{2S} }}} \right) + a_{0T}^{TT} J_{0 + 2} \left( {\frac{{\omega d_{ij} }}{{c_{0T} }}} \right),\end{eqnarray}

(6)
\begin{eqnarray} {\it{{S\!^{^\prime}}}}_{ij}^{TT} &=& a_{0S}^{RR} J_{0 + 2} \left( {\frac{{\omega d_{ij} }}{{c_{0S} }}} \right) + a_{1S}^{RR} J_{0 + 2} \left( {\frac{{\omega d_{ij} }}{{c_{1S} }}} \right)\nonumber\\ && +\, a_{2S}^{RR} J_{0 + 2} \left( {\frac{{\omega d_{ij} }}{{c_{2S} }}} \right) + a_{0T}^{TT} J_{0 - 2} \left( {\frac{{\omega d_{ij} }}{{c_{0T} }}} \right),\end{eqnarray}
where $$a_{{\rm mode}}^{kl}$$ is the amplitude of the mode in the kl component, dij is the distance between the ith and the jth stations. Jn is the nth Bessel function of the first kind, J0 − 2(x) = J0(x) − J2(x) and J0 + 2(x) = J0(x) + J2(x). We express phase velocities as a summation of B-spline functions, gm(f),
$$c_{{\rm mode}} \left( f \right) = \mathop \sum \limits_m p_m g_m \left( f \right),$$
and determine the values of coefficients, pm. We estimated the optimal phase velocity that minimizes the mean square misfit between observed and synthetic cross-spectra by the simulated annealing method (Ingber 1989): We assumed 40 000 combinations of pm, calculated the optimal amplitude of each mode and each component, $$a_{{\rm mode}}^{kl}$$, for each frequency for each assumed combination of pm by the least squares method, and searched for the minimum of the misfit. We only used pairs with interstation distance shorter than 250 km for PP, PZ and ZZ components for reducing the number of CCFs and avoiding the unknown phase in the PZ component (shown by an arrow in Fig. 2). The total number of CCFs used for the analysis is summarized in Table 2.

Fig. 6 shows the synthetic and observed CCFs. Phases are matched very well and three modes of Rayleigh waves are fitted even at short distances. The amplitudes are not matched well especially for the ocean-mode Rayleigh waves at a period range of 3–5 s. One possible reason is that the amplitude response of DPGs at this period range is not well determined as already described in Section 2. It does not affect the later estimation of structure and anisotropy in the solid Earth based on phase velocities of solid-mode Rayleigh waves. Fig. 7(a) shows the estimated phase velocities. The phase velocity of 0T mode is almost constant, whereas the phase velocities of 0S and 1S modes suddenly decrease at periods of ∼14 and ∼5 s, respectively. The error of phase-velocity measurement is estimated by the bootstrap method (Efron 1979). We first made 100 aggregates of CCFs by randomly selecting CCFs from all CCFs allowing duplicate selections. The number of CCFs in each aggregate is the same as the number of all CCFs. We then estimated 100 dispersion curves from 100 different aggregates of CCFs, and obtained the error of phase-velocity measurement from one standard deviation of 100 dispersion curves. The numbers of B-spline functions were chosen by the trial and error and summarized in Table 3. If we decrease the number of splines, the fitting between the observed and synthetic waveforms becomes worse. If we increase the number, on the other hand, the uncertainty becomes larger.

Figure 6.

Time-symmetric components of observed CCFs (black) and synthetic CCFs (red). Dashed pink lines show the arrival of each mode at a typical period. Although the gradients of pink lines for RR component are same as those for PP, PZ and ZZ components, they may seem different due to the difference between distance coverage.

Figure 6.

Time-symmetric components of observed CCFs (black) and synthetic CCFs (red). Dashed pink lines show the arrival of each mode at a typical period. Although the gradients of pink lines for RR component are same as those for PP, PZ and ZZ components, they may seem different due to the difference between distance coverage.

Figure 7.

(a) Measured phase velocities of Rayleigh (red circles) and Love (blue triangles) waves, and model phase velocities (black curves). (b) Model group velocities. (c) One-dimensional S-wave and P-wave velocity model shown by solid and dashed lines, respectively. (d) Histograms showing the uncertainty of the top depths of the two crustal layers and the mantle. (e) Same as (d) but for S-wave velocities in the three layers.

Figure 7.

(a) Measured phase velocities of Rayleigh (red circles) and Love (blue triangles) waves, and model phase velocities (black curves). (b) Model group velocities. (c) One-dimensional S-wave and P-wave velocity model shown by solid and dashed lines, respectively. (d) Histograms showing the uncertainty of the top depths of the two crustal layers and the mantle. (e) Same as (d) but for S-wave velocities in the three layers.

Table 3.

Range of analysis for four modes.

Mode 0S 1S 2S 0T
Period range of analysis 3–40 s 3–12 s 3–5s 3–40 s
Number of spline functions for measurement 14
Number of measurements for 1-D inversion 23 13 23
Mode 0S 1S 2S 0T
Period range of analysis 3–40 s 3–12 s 3–5s 3–40 s
Number of spline functions for measurement 14
Number of measurements for 1-D inversion 23 13 23

### Inversion for 1-D isotropic structure

Using average phase velocities, we estimated the 1-D isotropic structure in the crust and uppermost mantle beneath the area of the array by the simulated annealing method (Ingber 1989). The numbers of measurements used for the inversion are summarized in Table 3. The model parameters are the depth of the seafloor (d1), the depths of the bottom of two crustal layers (d2 and d3), S-wave velocity in the two crustal layers (v1 and v2), and the S-wave velocity in the uppermost mantle from the Moho to a depth of 220 km (v3). We assumed monotonic increase of S-wave velocity with depth. With the available period range, there is little resolution of mantle structure deeper than 40 km therefore we assumed a single mantle layer. The density (ρ) and P-wave velocity (VP) in the crust are scaled to S-wave velocity (VS) by the scaling relationship based on natural rocks (Christensen and Salisbury 1975) as:

\begin{eqnarray*} &&V_P = 1.75\,V_S + 0.375\,{\rm km}\;{\rm s}^{ - 1} ,\\ &&\rho = 0.5\;V_S\, {\rm g}\;{\rm cc}^{ - 1} {\rm (km}\;{\rm s}^{ - 1} {\rm )}^{ - 1} {\rm + 1}{\rm .25\,g}\;{\rm cc}^{ - 1} {\rm .}\end{eqnarray*}

These assumptions give the ratio between P and S waves in the crust larger for lower S-wave velocities. The ratio is set to be 1.73 in the mantle according to an experimental study by Anderson and Bass (1984). Other parameters such as the density in the mantle and attenuation coefficients were fixed to the PREM model (Dziewonski and Anderson 1981). We used the Fortran package DISPER80 (Saito 1988) with the physical dispersion (Kanamori and Anderson 1977) to calculate phase velocities of Rayleigh and Love waves from eigenperiods of spheroidal and toroidal modes, respectively.

Fig. 7(c) shows the obtained structure. The phase and group velocities corresponding to the obtained model are shown in Figs 7(a) and (b), respectively. The model phase velocities match well to the measured phase velocities, fitting within the range of error bars. As average phase velocities of both the Rayleigh and Love waves were satisfied with an isotropic structure, there is no requirement for radial anisotropy given the limited period range measured. The group velocity of 1S mode is lower than that of 0S mode at periods 3–5 s, and verifies our interpretation that the 0S mode arrives faster than 1S mode in Figs 2 and 6. We estimated the model uncertainty by estimating 100 models corresponding to 100 pairs of dispersion curves obtained by the bootstrap method, and show the histograms in Figs 7(d) and (e). As S-wave velocities in the second and third crustal layers (v2 and v3) are close to each other (Fig. 7e), we cannot recognize the boundary between the two layers in Fig. 7(c). The depth of the seafloor and the S-wave velocity in the mantle are well determined, whereas other parameters have larger uncertainty due to trade-off between parameters. Fig. 8 shows sensitivity of phase velocities to velocity structures. The sensitivity of phase velocities (c) at a frequency of f to P- or S-wave velocities (α or β) at a depth of z, Kα(z; f) and Kβ(z; f), is defined by

(7)
$$\frac{\Delta c(f)}{c(f)} = \int \left[ {\frac{{\Delta \alpha \left( z \right)}}{{\alpha \left( z \right)}}K_\alpha \left( {z; f} \right) + \frac{{\Delta \beta \left( z \right)}}{{\beta \left( z \right)}}K_\beta \left( {z; f} \right)} \right]{\rm d}z,$$
where Δc is the perturbation of phase velocity when Δα and Δβ are added to P- and S- wave structures, respectively. The ocean-mode Rayleigh waves mainly have sensitivity to P-wave velocity in the ocean, whereas the solid-mode Rayleigh wave and the fundamental-mode Love wave mainly have sensitivity to S-wave velocity in the solid Earth. The ocean mode of fundamental-mode Rayleigh wave corresponds to the Scholte wave, which propagates near the solid–liquid interface (Yao et al. 2011).

Figure 8.

Sensitivities of phase velocities of surface waves to (a) S-wave velocity and (b) P-wave velocity at periods of 20 s (thin solid lines), 10 s (thick solid lines), 5 s (broken lines) and 3 s (dotted lines). Two horizontal lines in each frame show the depths of the seafloor and the Moho.

Figure 8.

Sensitivities of phase velocities of surface waves to (a) S-wave velocity and (b) P-wave velocity at periods of 20 s (thin solid lines), 10 s (thick solid lines), 5 s (broken lines) and 3 s (dotted lines). Two horizontal lines in each frame show the depths of the seafloor and the Moho.

### Measurement of phase velocity anomalies

We finally measured phase velocity anomalies and amplitude anomalies of 0 S and 1 S modes at periods of 5–25 s by fitting the synthetic CCF to the observed CCF for each pair of stations and each bandpass filter. The frequency range for each filter is from 0.9f0 to 1.1f0, where f0 is the mean frequency. We used PP, PZ and ZZ components. The synthetic cross-spectra is,

(8)
\begin{eqnarray} {\it {S\!\!^{^\sim}}}_{ij}^{kl} &=& \left(1 + \zeta_{ij}^{kl} \right)a_{0S}^{kl} J_0 \left( {\frac{{\omega d_{ij} }}{{c_{0S} \left(1 + \gamma _{ij}^{kl} \right)}}} \right) \nonumber\\ &&+\, \left(1 + {\it{\zeta\!^{^\prime}}}_{ij}^{kl} \right)a_{1S}^{kl} J_0 \left( {\frac{{\omega d_{ij} }}{{c_{1S} \left(1 + {\it{\gamma\!^{^\prime}}}_{ij}^{kl} \right)}}} \right)\end{eqnarray}
where $$\zeta_{ij}^{kl}$$ and $$\gamma_{ij}^{kl}$$ are the amplitude anomaly and phase-velocity anomaly of the 0S mode, whereas $${\it{\zeta\!^{^\prime}}}_{ij}^{kl}$$ and $${\it{\gamma \!^{^\prime}}}_{ij}^{kl}$$ are those of the 1S mode. We obtained synthetic and observed waveforms by applying band-pass filters in the frequency domain, and by calculating inverse Fourier transforms.

In calculating the misfit between the synthetic and observed waveforms, the above envelope function was used as the weighting function, in essence using group velocity windowing to isolate each mode, and estimated optimal values of $$\zeta_{ij}^{kl}$$, $$\gamma_{ij}^{kl}$$, $${\it{\zeta\!^{^\prime}}}_{ij}^{kl}$$ and $${\it{\gamma\!^{^\prime}}}_{ij}^{kl}$$. For several pairs of stations, more than two components exist, and give different phase-velocity anomalies. We choose the value with the highest signal-to-noise ratio, RSN, where the signal amplitude is obtained from the root mean square of CCFs in the time window including the analysing mode, and the noise amplitude is obtained from the time window without the mode.

The phase-velocity anomalies shown in Fig. 9(a) seem to reflect the water-depth variations: Large deviations occur at periods of ∼10–20 s for 0S mode, where the transition from ocean-mode to solid-mode occurs. We estimated the effect of water depth, γd, by modifying our final 1-D model to various water depths, and by calculating phase velocities corresponding to the models. Fig. 9(b) shows the correction value with respect to the water-depth of 5.9 km as a function of period and water-depth, which is the largest at periods of ∼10–15 s. For estimating the correction term for each pairs of stations, we estimated the correction value as a function of longitude and latitude based on the sea topography by ETOPO 2.0, and averaged the values within the area where the distance from the great-circle path between pairs of stations is less than 10 km. Fig. 9(c) shows the phase-velocity anomalies after the correction, γ − γd. The deviation is much smaller than before the correction (Fig. 9a) especially at periods of 10–17 s for 0S mode. We hereafter only consider anomalies in the eastern array considering a smaller number of available measurements and larger uncertainties of each measurement in the western array.

Figure 9.

(a) Phase-velocity anomalies as a function of period. Solid black lines and dashed red lines correspond to pairs in eastern and western arrays, respectively. (b) The perturbations of phase velocities with respective to the water-depth of 5.9 km. (c) Phase-velocity anomalies after the correction of water-depth.

Figure 9.

(a) Phase-velocity anomalies as a function of period. Solid black lines and dashed red lines correspond to pairs in eastern and western arrays, respectively. (b) The perturbations of phase velocities with respective to the water-depth of 5.9 km. (c) Phase-velocity anomalies after the correction of water-depth.

Fig. 10 shows the phase-velocity anomalies with error smaller than 3 per cent as a function of azimuth of the great-circle path between the pair of stations. The error of each phase-velocity anomaly, Δγ, is estimated by considering three factors: signal to noise ratio of each CCF, RSN, the uncertainty of the clock delay or instrumental response, Δτ′ = 0.2s (Fig. 5c), and the azimuthal dependence of source intensity, B(ϕ), where ϕ is the backazimuth. The presence of noise in each CCF causes bias to the phase of signal by about tan − 1[1/(2RSN)]. The error in the phase-velocity anomaly is then Δγ1 = 1/(2RSNωτ) = 1/(4πRSN) × (λ′/d), where τ = d/u is the travel time and u is the model group velocity. The value, λ′ = u/f, is almost equivalent to the wavelength, λ = c/f. The error of the clock and response can be similarly estimated to be Δγ2 = Δτ′/τ.

Figure 10.

(a) Phase-velocity anomalies of the 0S mode in the eastern array as a function of azimuth of one station from another station. (b) Same as (a) but for the 1S mode. (c) Map showing direction of maximum velocity, ϕmax, for the solid-mode Rayleigh waves: 0S mode at periods of 10–21 s (red bars) and 1S mode at periods of 5–9 s (blue bars).

Figure 10.

(a) Phase-velocity anomalies of the 0S mode in the eastern array as a function of azimuth of one station from another station. (b) Same as (a) but for the 1S mode. (c) Map showing direction of maximum velocity, ϕmax, for the solid-mode Rayleigh waves: 0S mode at periods of 10–21 s (red bars) and 1S mode at periods of 5–9 s (blue bars).

The inhomogeneous source distribution of ambient noise may bias the phase-velocity measurements. Weaver et al. (2009) approximated the effect to be Δγ3 = B/8π2B × (λ′/d)2, which becomes $${\rm \Delta }\gamma _3 = \sum \limits_{n = 1}^{\infty} n^2 B_n /8\pi ^2 B \times ( {\lambda}^{\prime} /d )^2$$ for the source distribution of $$B(\phi) = B_0 + \sum \limits_{n = 1}^\infty B_n \cos n( \phi + \phi _n )$$. If the source distribution is B(ϕ) = 1 + 0.1cos 4ϕ, for example, the uncertainty is Δγ3 ∼ 2 × (λ′/d)2 per cent. Cox et al. (1984) showed that the Fourier components of odd orders (B1, B3, ⋅⋅⋅) and even orders (B2, B4, ⋅⋅⋅) only affect time-asymmetric and time-symmetric components of CCFs, respectively. As we only use time-symmetric components of CCFs, the effects of Fourier components of odd orders can be neglected. For the time-symmetric components used in this study, the effect of even orders (B2, B4, ⋅⋅⋅) needs to be considered, but is difficult to estimate in this study due to uncertain amplitude responses of DPGs. We therefore roughly estimated the error based on the source distribution in southern California estimated by Harmon et al. (2010). They estimated the ratio between Fourier components to be B2/B0 = 18 per cent, B4/B0 = 6 per cent, B6/B0 = 4 per cent and B8/B0 = 2 per cent. The mean error of phase-velocity anomalies is then Δγ3 = 2.8(λ′/d)2 per cent. We assume that this value approximates the upper limit for the area of our study. The order of magnitude of the source heterogeneity in the area of our study is expected to be the same or smaller than in southern California because the NW Pacific is surrounded by sources of ambient noise, the ocean and coastlines. The total error is the summation of all errors, Δγ = Δγ1 + Δγ2 + Δγ3. Although the estimated error, Δγ3, depends on the assumption of source distribution, the effect is negligible compared to the large error, Δγ1, due to low signal to ratio of ∼10. Even if the error, Δγ3, is twice or half as large as our estimate, it does not change the discussion below.

The anomalies shown in Fig. 10 tend to be high in the northeast and southwest directions, whereas the anomalies are low in the southeast and northwest directions. Fig. 10(c) shows the azimuth of maximum phase velocity, ϕmax, estimated by fitting A0 + A2cos2(ϕ − ϕmax) to the obtained phase-velocity anomalies for each mode and each frequency range. The minimum number of measurements for the fitting is 10 and 5 for 0S and 1S modes, respectively. The difference between the minimum numbers reflects the difference between available components: PP, PZ and ZZ components for 0S mode, and ZZ component for 1S mode. The uncertainties of ϕmax and the intensity of anisotropy, A2, are estimated by the bootstrap method (Efron 1979) as done for the estimation of the measurement error of phase velocities in Section 3.3. Although the small number of measurements for the 1S mode is not suitable for the bootstrap method, we can roughly evaluate the uncertainty caused by one anomalous measurement if it existed. Fig. 11 shows the estimated values, their uncertainties and the number of phase-velocity anomalies used for the estimation. The uncertainty for the 1S mode is larger than that for the 0S mode especially for the intensity of anisotropy due to smaller number of available phase-velocity anomalies. The peak-to-peak amplitude of anisotropy, 2A2, is ∼2–4 per cent for the solid-mode Rayleigh waves, whereas the values are less than 1 per cent for the ocean-mode Rayleigh waves. The directions of maximum velocity shown in Fig. 10(c) correspond to the values with errors smaller than 30°, and range from N16°E to N47°E, which are almost perpendicular to the direction of the magnetic lineations (N35°E).

Figure 11.

(a) Intensity of azimuthal anisotropy, (b) the azimuth of maximum phase velocity and (c) the number of measurements for each mode and each period band.

Figure 11.

(a) Intensity of azimuthal anisotropy, (b) the azimuth of maximum phase velocity and (c) the number of measurements for each mode and each period band.

## DISCUSSION

### Azimuthal anisotropy and lateral heterogeneity

We interpreted the phase-velocity anomalies by considering azimuthal anisotropy in Section 3.5. The fastest direction is almost perpendicular to magnetic lineations. This kind of anisotropy is frequently observed in oceanic regions, and can be interpreted to be caused by the mantle flow perpendicular to the ancient mid ocean ridge (Raitt et al. 1969; Francis 1969). Here we consider another possible interpretation of the lateral heterogeneity of phase-velocity anomalies. Fig. 12 shows the phase-velocity anomalies as a function of longitude for the 0S mode at a period of 18.8 s and for the 1S mode at a period of 6 s. For both modes, the phase-velocity anomalies are about 3 per cent stronger in the eastern part of eastern array than those in the western part of the eastern array. The variation in phase velocity at these periods becomes about half the variation in S-wave velocity at depths of ∼10–50 km. A phase-velocity anomaly of about 3 per cent, therefore, would indicate an S-wave velocity anomaly of about 6 per cent at depths of 10–50 km in the uppermost mantle. If we assume the percentage of the P-wave velocity anomaly to be the same as the S-wave velocity anomaly, the estimated S velocity anomaly becomes about 5 per cent.

Figure 12.

Phase-velocity anomalies of the fundamental-mode Rayleigh wave in the eastern array as a function of longitude at the centre between two stations.

Figure 12.

Phase-velocity anomalies of the fundamental-mode Rayleigh wave in the eastern array as a function of longitude at the centre between two stations.

Both lateral heterogeneity and azimuthal anisotropy can explain the phase-velocity anomalies, and the identification of the cause of the phase-velocity anomalies from current data is impossible. We, however, consider that the azimuthal anisotropy is the main cause of the phase-velocity anomalies. An S-wave velocity anomaly of 5 per cent corresponds to the change of rock composition or the temperature variation of 500 °K if we use the temperature derivative of olivine crystals by Liu et al. (2005). Such large lateral heterogeneity over a horizontal distance of 50–100 km is not expected in this region of almost uniform age and water depth. In comparison, large azimuthal anisotropy of up to 10 per cent is observed in various oceanic regions by refraction surveys, especially in the region with large shear deformation due to fast spreading (e.g. Song and Kim 2012) including the western array of PLATE project (Oikawa et al. 2010). Pn anisotropy in both the eastern and western arrays has been estimated to be 8–10 per cent by using high frequency phases generated by earthquakes in the western Pacific trenches (Shintaku et al. 2014). The seafloor in the eastern array, where we estimated azimuthal anisotropy in this study, has a high half spreading rate of ∼5 cm yr–1 (Nakanishi et al. 19891992; Müller et al. 2008). The presence of strong azimuthal anisotropy is, therefore, not surprising in this region.

### Intensity of azimuthal anisotropy

The intensity of azimuthal anisotropy is difficult to discuss especially when heterogeneity exists. Here we estimate the intensity by assuming laterally homogeneous structure and by assuming that the direction of maximum velocity is perpendicular to magnetic lineations. Fig. 13(a) shows the peak-to-peak intensity of azimuthal anisotropy, 2A2, obtained by fitting $$A_0 + A_2 {\rm cos}2( {\phi - \phi^{\prime} _{\rm max} } )$$ to the phase-velocity anomalies, where $$\phi^{\prime}_{{\rm max}} = 35^\circ$$. Two curves show the model intensity for the models shown in Fig. 12(b). We assumed an isotropic ocean and crust at depths of 0–13.3 km, and anisotropic mantle at depths of 13.3–220 km. As the azimuthal anisotropy of P-wave velocity slightly affects the azimuthal anisotropy of surface-wave velocity, we further assumed that the intensity of P-wave anisotropy is 1.3 times the intensity S-wave anisotropy based on the ratio between Pn- and Sn- wave anisotropy obtained in the northwest of Shatsky Rise by Shinohara et al. (2008).

Figure 13.

(a) Each point shows the peak-to-peak intensity of azimuthal anisotropy, 2A2, for the phase velocities of Rayleigh waves as a function of period. The fastest azimuth is assumed to be 35°. Horizontal and vertical error bars show the uncertainty and the range of bandpass filter, respectively. Curves show the theoretical values for the model shown in (b). (b) A model for the peak-to-peak intensity of azimuthal anisotropy of S-wave velocity.

Figure 13.

(a) Each point shows the peak-to-peak intensity of azimuthal anisotropy, 2A2, for the phase velocities of Rayleigh waves as a function of period. The fastest azimuth is assumed to be 35°. Horizontal and vertical error bars show the uncertainty and the range of bandpass filter, respectively. Curves show the theoretical values for the model shown in (b). (b) A model for the peak-to-peak intensity of azimuthal anisotropy of S-wave velocity.

We first fit the intensity of phase-velocity anisotropy by a constant model with an intensity of S-wave anisotropy of 7 per cent (solid lines in Fig. 13). The P-wave anisotropy assumed from the ratio between P- and S-wave anisotropy is then 9 per cent, which is slightly larger than the amplitude of the 2ϕ component of Pn anisotropy in the eastern array of about 5–7 per cent, but similar to the total peak-to-peak Pn anisotropy including both 2ϕ and 4ϕ components (Shintaku et al. 2014). The constant model is insufficient to fit the period dependence of intensity in Fig. 13(a). The observed intensity tends to be larger at shorter periods, especially shorter than 6 s for the 1S mode. The discrepancy indicates that the intensity of anisotropy is larger than the constant model at shallow depths of ∼10–20 km, that is, in the crust or in the uppermost mantle. From the observations of anisotropy in other oceanic regions, the intensity of S-wave anisotropy should be less than 30 per cent in the top 500 m of the oceanic crust, less than 2 per cent in the crustal layers below, and less than10 per cent in the mantle (e.g. Barclay & Toomey 2003; Song & Kim 2012). In the range, however, no 1-D model could fit the observed period dependence of the intensity of anisotropy. The discrepancy indicates the possibility of lateral heterogeneity, which cannot be determined from small number of phase-velocity anomalies in this study especially for the 1S mode (Fig. 11). We need to deploy OBSs more densely in laterally homogeneous seafloor for the future discussion of the depth dependence of intensity of anisotropy.

## CONCLUSION

We analysed surface waves in ambient noise by using eight OBSs and nine DPGs deployed southwest of the Shatsky Rise in NW Pacific Ocean. To achieve phase-velocity accuracy of ∼1 per cent, CCFs are used to estimate clock delays, and the instrument phase responses of DPGs. The CCF fitting analysis based on the SPAC method by Aki (1957) allowed us to analyse multimode surface waves including fundamental-, first higher- and second higher-mode Rayleigh waves and fundamental-mode Love waves even for the short interstation distances of ∼60 km. After estimating the effect of water-depth and the estimation error of phase-velocity anomalies, we estimated azimuthal anisotropy of S-wave velocity in the mantle. The direction of maximum velocity at depths of ∼10–50 km is perpendicular to the magnetic lineations, and indicates fossil mineral alignment in the lithosphere that is perpendicular to the ancient ridge spreading centre. Peak-to-peak shear wave azimuthal anisotropy of ∼7 per cent in the mantle is required to fit the average Rayleigh wave anisotropy of solid Earth modes. A larger seismic array with a better station distribution is needed to confirm the depth dependence of anisotropy.

We thank Hitoshi Kawakatsu, Shingo Watada and Takehi Isse in University of Tokyo for discussion and comments about the method and interpretation. This work was supported by Grant-in-Aids for JSPS Fellows (23-8157, 26-847) and Scientific Research (22000003) and by the National Science Foundation under Grant OCE 0648387. Seismometers were supplied by the Ocean Bottom Seismograph Instrument Pool (OBSIP). We thank the members of the SIO and L-DEO instrument centers and the captains and crews of the RV Melville and RV Moana Wave. First author (AT) visited Brown University two months for this study by the support of Institutional Program for Young Researcher Overseas Visits of JSPS. We used GMT (Wessel & Smith 1991) and SAC2000 (Goldstein & Snoke 2005) software during this study.

## REFERENCES

Aki
K.
Space and time spectra of stationary stochastic waves, with special reference to microtremors
Bull. Earthq. Res. Inst.
,
1957
, vol.
35

3
(pg.
415
-
456
)
Anderson
D.L.
Bass
J.D.
Mineralogy and composition of the upper mantle
Geophys. Res. Lett.
,
1984
, vol.
11

7
(pg.
637
-
640
)
Araki
E.
Sugioka
H.
Calibration of deep sea differential pressure gauge
JAMSTEC-R IFREE
,
2009
, vol.
9
(pg.
141
-
148
)
Barclay
A.H.
Toomey
D.R.
Shear wave splitting and crustal anisotropy at the Mid-Atlantic Ridge, 35°N
J. geophys. Res.
,
2003
, vol.
108

B8
pg.
2378

doi:10.1029/2001JB000918
Beaty
K.S.
Schmitt
D.R.
Sacchi
M.
Simulated annealing inversion of multimode Rayleigh wave dispersion curves for geological structure
Geophys. J. Int.
,
2002
, vol.
151

2
(pg.
622
-
631
)
Bensen
G.D.
Ritzwoller
M.H.
Shapiro
N.M.
Broadband ambient noise surface wave tomography across the United States
J. geophys. Res.
,
2008
, vol.
113

B5
(pg.
1
-
21
)
Christensen
N.I.
Salisbury
M.H.
Structure and constitution of the lower oceanic crust
Rev. Geophys.
,
1975
, vol.
13

1
(pg.
57
-
86
)
Cox
C.
Deaton
T.
Webb
S.
A deep-sea differential pressure gauge
J. Atmos. Ocean. Tech.
,
1984
, vol.
1

3
(pg.
237
-
246
)
Dziewonski
A.M.
Anderson
D.L.
Preliminary reference Earth model
Phys. Earth planet. In.
,
1981
, vol.
25

4
(pg.
297
-
356
)
Efron
B.
Bootstrap methods: another look at the jackknife
Ann. Stat.
,
1979
, vol.
7

1
(pg.
1
-
26
)
Ewing
W.M.
Jardetzky
W.S.
Press
F.
Elastic Waves in Layered Media
,
1957
McGraw Hill Book Campany Inc
Francis
T.J.G.
Generation of seismic anisotropy in the upper mantle along the mid-oceanic ridges
Nature
,
1969
, vol.
221

5176
(pg.
162
-
165
)
Forsyth
D.W.
Li
A.
Levander
A.
Nolet
G.
Array analysis of two-dimensional variations in surface wave phase velocity and azimuthal anisotropy in the presence of multipathing interference
Seismic Earth: Array Analysis of Broadband Seismograms
,
2005
American Geophysical Union
(pg.
81
-
97
pp.
Gouedard
P.
Seher
T.
McGuire
J.J.
Collins
J.A.
van der Hilst
R.D.
Correction of ocean-bottom seismometer instrumental clock errors using ambient seismic noise
Bull. seism. Soc. Am.
,
2014
, vol.
104

3
(pg.
1276
-
1288
)
Goldstein
P.
Snoke
A.
SAC availability for the IRIS community
DMS Electr. Newslett.
,
2005
, vol.
7

1

Hannemann
K.
Kruger
F.
Dahm
T.
Measuring of clock drift rates and static time offsets of ocean bottom stations by means of ambient noise
Geophys. J. Int.
,
2013
, vol.
196

2
(pg.
1034
-
1042
)
Harmon
N.
Forsyth
D.
Webb
S.
Using ambient seismic noise to determine short-period phase velocities and shallow shear velocities in young oceanic lithosphere
Bull. seism. Soc. Am.
,
2007
, vol.
97

6
(pg.
2009
-
2023
)
Harmon
N.
Rychert
C.
Gerstoft
P.
Distribution of noise sources for seismic interferometry
Geophys. J. Int.
,
2010
, vol.
183

3
(pg.
1470
-
1484
)
Ingber
L.
Very fast simulated re-annealing
Math. Comp. Modell.
,
1989
, vol.
12

8
(pg.
967
-
973
)
Kanamori
H.
Anderson
D.L.
Importance of physical dispersion in surface wave and free oscillation problems: review
Rev. Geophys.
,
1977
, vol.
15

1
(pg.
105
-
112
)
Lin
F.-C.
Ritzwoller
M.H.
Yang
Y.
Moschetti
M.P.
Fouch
M.J.
Complex and variable crustal and uppermost mantle seismic anisotropy in the western United States
Nat. Geosci.
,
2010
, vol.
4

1
(pg.
55
-
61
)
Liu
W.
Kung
J.
Li
B.
Elasticity of San Carlos olivine to 8 GPa and 1073 K
Geophys. Res. Lett.
,
2005
, vol.
32

16

doi:10.1029/2005GL023453
Mordret
A.
Landes
M.
Shapiro
N.M.
Singh
S.C.
Roux
P.
Barkved
O.I.
Near-surface study at the Valhall oil field from ambient noise surface wave tomography
Geophys. J. Int.
,
2013a
, vol.
193

3
(pg.
1627
-
1643
)
Mordret
A.
Shapiro
N.M.
Singh
S.
Roux
P.
Montagner
J.-P.
Barkved
O.I.
Azimuthal anisotropy at Valhall: the Helmholtz equation approach
Geophys. Res. Lett.
,
2013b
, vol.
40

11
(pg.
2636
-
2641
)
Müller
R.D.
Sdrolias
M.
Gaina
C.
Roest
W.R.
Age, spreading rates, and spreading asymmetry of the world's ocean crust
Geochem. Geophys. Geosyst.
,
2008
, vol.
9

4
(pg.
1
-
19
)
Nakanishi
M.
Tamaki
K.
Kobayashi
K.
Mesozoic magnetic anomaly lineations and seafloor spreading history of the northwestern Pacific
J. geophys. Res.
,
1989
, vol.
94

B11
(pg.
15 437
-
15 462
)
Nakanishi
M.
Tamaki
K.
Kobayashi
K.
Magnetic anomaly lineations from Late Jurassic to Early Cretaceous in the west-central Pacific Ocean
Geophys. J. Int.
,
1992
, vol.
109

3
(pg.
701
-
719
)
Nishida
K.
Two-dimensional sensitivity kernels for cross-correlation functions of background surface waves
C. R. Geosci.
,
2011
, vol.
343

8–9
(pg.
584
-
590
)
Nishida
K.
Kawakatsu
H.
Obara
K.
Three-dimensional crustal S wave velocity structure in Japan using microseismic data recorded by Hi-net tiltmeters
J. geophys. Res.
,
2008
, vol.
113

B10
pg.
B10302

doi:10.1029/2007JB005395
Oikawa
M.
Kaneda
K.
Nishizawa
A.
Seismic structures of the 154–160 Ma oceanic crust and uppermost mantle in the Northwest Pacific Basin
Earth Planets Space
,
2010
, vol.
62

4
(pg.
e13
-
e16
)
Raitt
R.W.
Shor
G.G.
Francis
T.J.G.
Morris
G.B.
Anisotropy of the Pacific upper mantle
J. geophys. Res.
,
1969
, vol.
74

12
(pg.
3095
-
3109
)
Saito
M.
Doornbos
D.J.
DISPER80: a subroutine package for the calculation of seismic normal-mode solutions
Seismological Algorithms
,
1988
(pg.
293
-
319
pp.
Sens-Schönfelder
C.
Synchronizing seismic networks with ambient noise
Geophys. J. Int.
,
2008
, vol.
174

3
(pg.
966
-
970
)
Shapiro
N.M.
Campillo
M.
Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise
Geophys. Res. Lett.
,
2004
, vol.
31

7
pg.
L07614

doi:10.1029/2004GL019491
Shinohara
M.
, et al.  .
Upper mantle and crustal seismic structure beneath the Northwestern Pacific Basin using a seafloor borehole broadband seismometer and ocean bottom seismometers
Phys. Earth planet. Inter.
,
2008
, vol.
170

1–2
(pg.
95
-
106
)
Shintaku
N.
Forsyth
D.W.
Hajewski
C.J.
Weeraratne
D.S.
Pn anisotropy in Mesozoic western Pacific lithosphere
J. geophys. Res.: Solid Earth
,
2014
, vol.
199

4
(pg.
3050
-
3063
)
Smith
D.B.
Ritzwoller
M.H.
Shapiro
N.M.
Stratification of anisotropy in the Pacific upper mantle
J. geophys. Res.
,
2004
, vol.
109

B11
(pg.
1
-
22
)
Snieder
R.
Extracting the Green's function from the correlation of coda waves: a derivation based on stationary phase
Phys. Rev. E
,
2004
, vol.
69
pg.
46610

doi:10.1103/PhysRevE.69046610
Song
T.-R.A.
Kim
Y.
Anisotropic uppermost mantle in young subducted slab underplating Central Mexico
Nat. Geosci.
,
2012
, vol.
5

1
(pg.
55
-
59
)
Takeo
A.
Nishida
K.
Isse
T.
Kawakatsu
H.
Shiobara
H.
Sugioka
H.
Kanazawa
T.
Radially anisotropic structure beneath the Shikoku Basin from broadband surface wave analysis of ocean bottom seismometer records
J. geophys. Res.
,
2013
, vol.
118

6
(pg.
2878
-
2892
)
Tanimoto
T.
Anderson
D.L.
Mapping convection in the mantle
Geophys. Res. Lett.
,
1984
, vol.
11

4
(pg.
287
-
290
)
Tromp
J.
Luo
Y.
Hanasoge
S.
Peter
D.
Noise cross-correlation sensitivity kernels
Geophys. J. Int.
,
2010
, vol.
183

2
(pg.
791
-
819
)
Weaver
R.
Froment
B.
Campillo
M.
On the correlation of non-isotropically distributed ballistic scalar diffuse waves
J. acoust. Soc. Am.
,
2009
, vol.
126

4
(pg.
1817
-
1826
)
Webb
S.C.
Deaton
T.K.
Lemire
J.C.
A broadband ocean-bottom seismometer system based on a 1-Hz natural period geophone
Bull. seism. Soc. Am.
,
2001
, vol.
91

2
(pg.
304
-
312
)
Weeraratne
D.S.
Forsyth
D.W.
Yang
Y.
Webb
S.C.
Rayleigh wave tomography beneath intraplate volcanic ridges in the South Pacific
J. geophys. Res.
,
2007
, vol.
112

B6
(pg.
1
-
18
)
Wessel
P.
Smith
W.H.F.
Free software helps map and display data
EOS, Trans. Am. geophys. Un.
,
1991
, vol.
72

41
(pg.
441
-
441
)
Yao
H.
van der Hilst
R.D.
Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to SE Tibet
Geophys. J. Int.
,
2009
, vol.
179

2
(pg.
1113
-
1132
)
Yao
H.
Gouédard
P.
Collins
J.A.
McGuire
J.J.
van der Hilst
R.D.
Structure of young East Pacific Rise lithosphere from ambient noise correlation analysis of fundamental- and higher-mode Scholte-Rayleigh waves
C. R. Geosci.
,
2011
, vol.
343
(pg.
571
-
583
)