## 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.* 1989, 1992). 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.

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.

Station | OBS |
---|---|

1 | T-240 |

4 | T-40 |

6 | – |

7 | T-240 |

8 | LDEO |

10 | – |

11 | T-240 |

12 | LDEO |

13 | – |

15 | T-40 |

16 | T-240 |

Station | OBS |
---|---|

1 | T-240 |

4 | T-40 |

6 | – |

7 | T-240 |

8 | 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,

*k*th component record at the

*i*th 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).

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.

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).

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.

### 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 (*c*_{mode}) 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,

*kl*component,

*d*is the distance between the

_{ij}*i*th and the

*j*th stations.

*J*is the

_{n}*n*th Bessel function of the first kind,

*J*

_{0 − 2}(

*x*) =

*J*

_{0}(

*x*) −

*J*

_{2}(

*x*) and

*J*

_{0 + 2}(

*x*) =

*J*

_{0}(

*x*) +

*J*

_{2}(

*x*). We express phase velocities as a summation of B-spline functions,

*g*(

_{m}*f*),

*p*. 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

_{m}*p*, calculated the optimal amplitude of each mode and each component, $$a_{{\rm mode}}^{kl} $$, for each frequency for each assumed combination of

_{m}*p*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.

_{m}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.

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 | 8 | 8 | 5 |

Number of measurements for 1-D inversion | 23 | 13 | 5 | 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 | 8 | 8 | 5 |

Number of measurements for 1-D inversion | 23 | 13 | 5 | 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 (*V _{P}*) in the crust are scaled to

*S*-wave velocity (

*V*) by the scaling relationship based on natural rocks (Christensen and Salisbury 1975) as:

_{S}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

*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).

### 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.9*f*_{0} to 1.1*f*_{0}, where *f*_{0} is the mean frequency. We used PP, PZ and ZZ components. The synthetic cross-spectra is,

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, *R _{SN}*, 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.

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, *R _{SN}*, 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/(2

*R*)]. The error in the phase-velocity anomaly is then Δγ

_{SN}_{1}= 1/(2

*R*ωτ) = 1/(4π

_{SN}*R*) × (λ′/

_{SN}*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}= Δτ′/τ.

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π^{2}*B* × (λ′/*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 (*B*_{1}, *B*_{3}, ⋅⋅⋅) and even orders (*B*_{2}, *B*_{4}, ⋅⋅⋅) 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 (*B*_{2}, *B*_{4}, ⋅⋅⋅) 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 *B*_{2}/*B*_{0} = 18 per cent, *B*_{4}/*B*_{0} = 6 per cent, *B*_{6}/*B*_{0} = 4 per cent and *B*_{8}/*B*_{0} = 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 *A*_{0} + *A*_{2}cos2(ϕ − ϕ_{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, *A*_{2}, 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, 2*A*_{2}, 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).

## 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.

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.* 1989, 1992; 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, 2*A*_{2}, 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).

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.