The misorientation of three‐component seismometers restricts the application of relevant seismic experiments such as ocean‐bottom seismometer (OBS) arrays. Previous orientation determination relied on an assumption that the propagation azimuth of seismic waves follows the great‐circle path (GCP) azimuth. This assumption may yield systematic errors in the estimated orientation, particularly when the ray paths are bent laterally due to velocity heterogeneity in the Earth. Here, we develop a new method for unbiasedly estimating the horizontal orientations of seismic sensors and apply this method to the Blanco transform fault OBS experiment. We first retrieve the orientations relative to the propagation azimuths from the recorded Rayleigh and P waveforms, and then determine the geographic north orientations by calculating the propagation azimuths via an Eikonal‐equation‐based phase‐tracking method that theoretically accounts for the effect of ray bending. Synthetics test validates that the phase‐tracking method can retrieve unbiased propagation azimuths of seismic waves. The final results derived from Rayleigh‐ and P‐wave polarization analyses with the respective phase‐tracked propagation azimuths are more consistent and the orientation errors are smaller, indicating the robustness and accuracy of this method. Comparing the orientations from our phase‐tracking method to those from the GCP assumption, the deviation can reach up to 8° between these two techniques in the study region. Subsequently, when orientations of the synthetics modeled from three‐dimensional elastic waveform simulation are deviated according to the GCP‐predicted orientations, we find nonnegligible bias in the phase and amplitude measurements that could reduce the accuracy and resolution of following inversion, which indicates the significance of our phase‐tracking method in accurate orientation of OBS arrays as well as inland seismic experiments.

Accurate orientation of three‐component seismometers is crucial to various seismic measurements, such as the vertical‐to‐horizontal (Z–H) amplitude ratios (Boore and Toksoz, 1969; Xu et al., 2020), the component‐differential phase delays (Bao and Shen, 2018; Bao et al., 2021), and receiver functions (Langston, 1979; Zhu and Kanamori, 2000; Cheng et al., 2022). The passive‐source experiments of ocean‐bottom seismometers (OBSs) have been rapidly advancing our insight into the dynamics of the oceanic crust and mantle, such as the expansion mechanism of midocean ridges (Forsyth et al., 1998; Zhang et al., 2022), the water cycle at subduction zones (Cai et al., 2018; Zhu et al., 2021), and the slip behavior of transform faults (Froment et al., 2014; Kuna et al., 2019; Shi et al., 2022). However, the unknown horizontal orientations of OBSs attributed to the free‐fall deployment restrict the utilization of the associated waveform‐based phase and amplitude measurements to explore earth structure. Moreover, temporary inland seismic experiments may also suffer the misorientation due to the lack of fiberoptic gyroscopes (Wang et al., 2016; Xu et al., 2018; Zeng et al., 2021).

Previous studies have developed a variety of methods to determine the correct seismometer orientation. Airgun shot arrays are typically effective to orient the horizontal components in active‐source OBS deployments (Anderson et al., 1987; Duennebier et al., 1987), but not common for most passive source broadband OBS experiments due to difficulties of cost and permission. Recent passive broadband OBS experiments have obtained instrument orientations via Rayleigh‐ and P‐wave polarization (hereafter called R‐pol and P‐pol, respectively) from large‐magnitude earthquakes (Stachnik et al., 2012; D’Alessandro et al., 2013; Doran and Laske, 2017; Scholz et al., 2017), ambient noise cross correlations (Zha et al., 2013), or receiver functions (Zheng et al., 2020). These methods estimate the seismic propagation orientation via polarization analyses and then rotate the propagation orientation to the geographic north orientation according to the propagation azimuth. Although similar north orientations were produced by the polarization analysis using different seismic phases, discrepancy remains between the results from the R‐pol and P‐pol analyses (Stachnik et al., 2012; Scholz et al., 2017). As a possible reason, previous studies commonly relied on an assumption that the propagation azimuth of recorded waves follows the great‐circle paths (GCPs). However, this assumption may yield considerable deviation to the estimated orientations when the true ray paths are bent laterally due to the velocity heterogeneity in the Earth, and may also cause the difference between the results from surface waves and body waves which ray paths may be bent to different azimuths. In contrast, by tracking the propagation of a phase within a seismic array based on the Eikonal equation (Aki and Richards, 1980), the true wave propagation azimuths can be retrieved more accurately (Lin et al., 2009), suggesting a new way for the more unbiased determination of OBS orientations when combined with the polarization analyses.

In this study, we present a joint method of an Eikonal‐equation‐based phase‐tracking and seismic wave polarization analyses, and then apply this method to an OBS array around the Pacific–Juan de Fuca plate boundary, the Blanco transform fault (BTF) OBS Experiment (Nábělek and Braunmiller, 2012). Our method, which differs from the conventional approach in that it considers the impact of ray bending on the estimated orientations, theoretically produces more accurate orientation results. We corroborate the flipping polarity of the OBS data and the determined orientations by comparing the corrected waveforms with the synthetic waveforms from three‐dimensional full wavefield simulation. Based on these synthetics, we also evaluate the improvement of the phase and amplitude measurements with the orientations determined from our method in contrast to the GCP assumption.

Data and R‐pol/P‐pol analyses

The three‐component seismic waveform data used in this study were recorded during the BTF Experiment (Nábělek and Braunmiller, 2012) from September 2012 to October 2013 (Fig. 1a) at 54 OBSs in total, including 30 Woods Hole Oceanographic Institution (WHOI) OBSs and 24 Scripps Institution of Oceanography (SIO) OBSs. First, we remove the instrumental responses and correct for the tilt and compliance noises for the vertical component (Bell et al., 2015; Janiszewski et al., 2019), and verify the polarity of the vertical component based on the Rayleigh‐wave records (Xu et al., 2018; Zhu et al., 2020). See Text S1, available in the supplemental material to this article.

We then calculate the orientation from Rayleigh and P waveforms relative to the propagation azimuth via the R‐pol and P‐pol analyses (Stachnik et al., 2012; D’Alessandro et al., 2013; Doran and Laske, 2017; Scholz et al., 2017; Yang et al., 2020). About 91 Mw>6 earthquakes and 1461 Mw>5 earthquakes are selected for the R‐pol and P‐pol analyses, respectively (Fig. 1b). We discard the low‐quality events according to low correlation coefficients or low linearity (Doran and Laske, 2017; Zhu et al., 2020), the remaining events are defined as the high‐quality events. Stations BS160, BS310, BS360, and BS470 have either no records or extremely high noise level so that we do not determine their orientations. See Text S2 and S3.

The polarization analyses call for a band‐pass filter that balances the low signal‐to‐noise ratio (SNR) at high frequencies and low sensitivity at low frequencies. The WHOI OBSs were designed for a broad frequency band of 0.01–30 Hz (Janiszewski et al., 2019; Yang et al., 2019), whereas the SIO OBSs were short‐period seismometers with less sensitivity to the frequencies lower than 8 Hz. Because teleseismic waves typically have low SNRs at relatively high frequencies, previous studies commonly used band‐pass filters such as 0.02–0.04 Hz, 0.03–0.07 Hz, and 0.05–0.12 Hz, during their determination of OBSs orientations (Stachnik et al., 2012; Scholz et al., 2017; Zhu et al., 2020). Here, we attempt various frequency bands and find that the SIO OBSs have more earthquake signals with a filter of 0.05–0.1 Hz, and thus use it to determine orientations for both WHOI and SIO OBSs.

Phase‐tracking method

We transfer the propagation azimuths to the geographic north orientations by allowing the off‐GCPs for both P and Rayleigh waves due to the heterogeneity of the Earth. The propagation azimuths of Rayleigh waves are determined from the horizontal gradients (at the station locations) of the phase travel‐time surface (T) tracking the Rayleigh‐wave phase travel times within a two‐dimensional seismic array (Lin et al., 2009):
(1)Tk^c.

Here, k^ is the unit wavenumber vector and represents the propagation azimuth, c is the phase speed here. Although equation (1) has been previously applied to structural imaging (commonly called Eikonal tomography) for surface waves (Lin et al., 2009; Lin and Ritzwoller, 2011; Jin and Gaherty, 2015; Bao et al., 2016; Janiszewski et al., 2019), we take its advantage to account for the off‐GCP propagation and estimate the propagation azimuths of Rayleigh waves. This technique, called the phase‐tracking method hereinafter, is different from the previous studies assuming the GCP as the ray path (Stachnik et al., 2012; Doran and Laske, 2017; Scholz et al., 2017). Moreover, we expand its application to measure the propagation azimuths of teleseismic P waves. This treatment is not to yield the full P‐wave vectors, because it is difficult to track the three‐dimensional ray paths of teleseismic body waves based on equation (1) with P‐wave recorded at stations on the surface. In contrast, the horizontal components of the P‐wave vectors are retrievable from the phase‐tracking method given a two‐dimensional seismic array, which essentially provides the horizontal propagation azimuths of P waves needed in our study.

Here, we calculate the gradient of the phase travel‐time surface based on the automated surface‐wave phase‐velocity measuring system (Jin and Gaherty, 2015) for surface waves and then estimate the propagation azimuths. The resulting orientation with respect to the north and the error estimated in the polarization analysis are determined from the mean and the median absolute deviation (Rousseeuw and Croux, 1993; Doran and Laske, 2017), respectively, for the events with high quality (Text S4). According to the resulting orientation, the original horizontal‐component waveforms are rotated to the geographic north and east orientations. We call the waveforms after this rotation the corrected waveforms hereinafter. We do a similar process for P waves except using the same time window as in the P‐pol analysis.

Benchmark and uncertainty estimation

We carry out a test to validate the feasibility and accuracy of the phase‐tracked azimuths. We first calculate synthetic waveforms based on a laterally homogeneous medium, called one‐dimensional synthetics hereinafter, then estimate the phase‐tracked azimuths from the one‐dimensional synthetics at the station locations. The azimuths from phase tracking the one‐dimensional synthetics are theoretically equal to the GCP‐predicted azimuths and thus provide a benchmark for estimating the true azimuths as well as their uncertainties. We employ the Mineos package (Masters et al., 2011) to calculate the one‐dimensional synthetics based on the global model preliminary reference Earth model, PREM (Dziewonski and Anderson, 1981), for the same events and stations as used in the R‐pol analysis. The corresponding earthquake mechanisms are obtained from the Global Centroid Moment Tensor Project (see Data and Resources). Figure 2a–c illustrates the observed waveforms and synthetics for an Mw 6.1 earthquake event (8 November 2012, Vancouver Island event 02:01:50 UTC, 49.231° N, 128.477° W, 13.7 km depth; see Data and Resources). The phase‐tracked azimuths are then compared with GCP‐predicted azimuths. The mean of the difference between the GCP‐predicted and the phase‐tracked azimuths approximates 0° (Fig. 2d), indicating that the phase‐tracking method can retrieve unbiased propagation azimuths of seismic waves in the OBSs array.

Although the uncertainty of the resulting OBS orientations can be partly estimated by the polarization analysis, the phase‐tracking method may introduce three types of errors due to the geometrical asymmetry of OBS arrays (Jin and Gaherty, 2015), the uneven coverage of events’ azimuthal, and the clock drift at each station. The quantities of these error sources rely on the particular array deployment and cannot be clearly divided. Thus, here we assess their overall uncertainty by employing the above one‐dimensional synthetics for the Blanco OBS array and the recorded events same as in our orientation determination. Gardner and Collins (2012) measured the clock drifts for the OBS deployments from 2010 to 2012, suggesting that the mean value of clock drifts is 1.3 ms/day with a median absolute deviation of 1 ms/day. Given that the OBSs get pre‐ and postdeployment clock checks, the maximum clock drift will occur in about half of the deployment cycle (Gardner and Collins, 2012). Based on the earlier criteria, we randomly generate clock drifts in the form of a Gaussian distribution and add them to the one‐dimensional synthetics with a bound of sampled drifts that increases over the deploying time. Then, we calculate the propagation azimuths of the synthetics via the phase‐tracking method. Similarly, we calculate the mean of the difference between the GCP‐predicted and the phase‐tracked azimuths. The uncertainty is determined as two times standard deviation of the mean (95% confidence interval of the mean). We repeat the earlier processes of adding clock drift and calculating uncertainty multiple times. The uncertainty is approximately 1°, suggesting that the clock drift in the OBSs has little effect on the estimated propagation azimuths. Finally, we estimate a total error, Errortotal, by summing up the errors produced by the phase‐tracking method, Errorptm, and the corresponding polarization analyses, Errorpa, based on L2 norm: Errortotal=Errorptm2+Errorpa2.

The angle difference between the phase‐tracked azimuths (k^) and GCP‐predicted azimuths of Rayleigh wave vary with the distribution of the events used for phase tracking. Although these angles are within ±6° for many events (such as an earthquake near New Zealand, Fig. 3a), some events along the east Pacific subduction zone yield ∼40° angle difference (Fig. 3c). In contrast to the significant angle difference (∼40°) for Rayleigh waves, the angle difference for P waves mainly ranges within ±3° (Fig. 3b). Because the phase‐tracking method requires relatively high coherence between the waveforms of nearby stations (Jin and Gaherty, 2015), we are not able to estimate the phase‐tracked azimuth for some event‐station pairs (Table 1 and Fig. 4a,b).

The final orientations and the errors of the BTF OBSs are summarized in Table 1 and Figure 4a–d. Errors from the R‐pol and P‐pol analyses generally vary between 2° and 16° (Table 1 and Fig. 4c,d), yet can be as high as ∼50° for individual stations (e.g., BS520) presumably owing to the high noise on their horizontal components.

Generally, the orientations derived from the R‐pol and P‐pol analyses show a relatively good agreement for most stations (Fig. 4b), which is consistent with previous studies (Scholz et al., 2017; Zhu et al., 2020). We only account for the stations for which orientations are calculated from more than 10 high‐quality events (Wang et al., 2016). All the SIO OBSs suffered vertical‐component polarity flipping, besides the station BS250 in which the polarity of one horizontal component also flips.

Advantages of the phase‐tracking method

Applying the phase‐tracked azimuths to the orientation estimation results in unbiased, more robust, and more accurate orientations (Table 1 and Fig. 4). Comparing the orientations calculated from the phase‐tracked azimuths and the GCP‐predicted azimuths, the deviation can reach up to 8° in the study region (Fig. 4e,f). The deviations of individual events can be significantly larger than 8° because the orientation is determined by averaging multiple events with different azimuths. For example, the deviations of the events located in the southeast and northwest Pacific subduction zones can be even over 20° (Fig. 3c,d), within which the phase‐tracked azimuths are typically larger than the GCP‐predicted azimuths for the earthquakes in southeast Pacific (Fig. 3c), yet the GCP‐predicted azimuths are typically larger than the phase‐tracked azimuths for earthquakes in northwest Pacific (Fig. 3d). This suggests that the phase‐tracking method are particularly crucial for both earthquakes and seismic arrays located near subduction zones in which many ray paths are bent by the velocity anomaly of slab particularly with limited events’ azimuthal coverage.

The orientations determined from the R‐pol and P‐pol analyses become more consistent after applying the phase‐tracking azimuths instead of the GCP‐predicted azimuths (Fig. 4g,h). The mean of the orientation difference between the R‐pol and P‐pol analyses approximates to 0° for the phase‐tracking method, but 3° for the GCP prediction (Fig. 4g,h). This consistency demonstrates that the orientations from the phase‐tracking method with either R‐pol or P‐pol analyses are more robust than those from the GCP‐predicted azimuths. However, the phase‐tracking method with the R‐pol analysis may be practically more effective because the P‐pol analysis commonly lacks high‐quality events (Table 1).

The phase‐tracking method does not only yield more robust orientations, but also produces smaller orientation errors, in which the mean value of the error reduced by about 3° (Fig. 4c,d). As Figure 5 illustrates, most events tend to be more aggregated after applying the phase‐tracked azimuths, indicating that the phase‐tracking method reduces the uncertainty in the estimated orientation. The smaller orientation errors imply more accurate orientations.

Comparison with three‐dimensional synthetic waveform

We calculate synthetic waveforms based on a three‐dimensional elastic wavefield forward modeling to corroborate our results, including the validation of the observed polarity flipping of OBS data (Zhu et al., 2020). Here, we select the same earthquake event as in the Benchmark and uncertainty estimation section (Fig. 2a) and compare the three‐component synthetics to the observed waveforms on OBSs. The synthetics are simulated using a three‐dimensional collocated grid finite‐difference algorithm in the curvilinear coordinates (Zhang et al., 2012; Sun et al., 2018, 2021) that takes complex topography, bathymetry, and heterogeneous Earth media into account. The Earth model for the forward modeling is the SPiRal that was derived from both body‐wave travel times and surface‐wave dispersion measurements (Simmons et al., 2021). Although the SPiRal model contains anisotropic parameters, our simulation uses the isotropic shear wavespeed (Voigt average, VS=23VSV2+13VSH2). The horizontal grid spacing used in wave simulation is 1.5 km, and the vertical grid spacing varies from 0.13 km at the surface to 1.3 km at 400 km depth. The model parameters and topography are interpolated into each grid using cubic‐spline interpolation. The simulation takes 160 central processing unit cores and 30 wall‐clock hours.

The high cross‐correlation coefficients between the synthetics and the corrected waveforms on all three components for most OBSs validate the polarity flipping in the SIO OBSs and reinforce the reliability of our orientation results (Fig. 6). The corrected waveforms are determined from the R‐pol analysis with the phase‐tracking method. Similarly, we use the P‐pol analysis with the phase‐tracking method to determine the corrected waveforms for the stations with high‐quality events, which correspondingly present high cross‐correlation coefficients with the synthetics.

Effect of inaccurate orientation on waveforms

We evaluate the impact of inaccurate orientation due to the GCP assumption on the corrected waveforms by simulating the bias of the GCP‐prediction orientations with respect to the phase‐tracked orientations shown in Table 1 upon the three‐dimensional synthetics. Considering that about half of the stations do not yield orientations from the P‐pol analysis due to the lack of high‐quality events (Table 1), we determine the orientation bias from the difference between the R‐pol analysis with the phase‐tracking method and the R‐pol analysis with the GCP assumption. We rotate the orientations of the two horizontal components of the synthetics and calculate the differences in the measured phase and amplitude, respectively, between the rotated and the original synthetics. The phase delays and amplitude ratios determined by cross correlation (Bao and Shen, 2018) between the rotated and the original synthetics, thus denote the effect of inaccurate orientation on the corrected waveforms. Table S1 demonstrates the variations of the phase and amplitude of seismic waves due to the inaccurate orientation.

Inaccurate orientation has a larger effect the surface waves than the P waves (Fig. 7). Although mostly less than 0.01 s, the phase delays at some stations can reach up to 0.05 s (Fig. 7a,b). Given that the radial‐ and vertical‐component phase delays are about 0.5 s in the BTF OBS array after eliminating an 1/4 period difference, an extra error of 0.05 s produced by inaccurate orientation should not be neglected for accurately measuring the component‐differential phase delays (Bao and Shen, 2018; Bao et al., 2021). Meanwhile, the effect of inaccurate orientation on amplitude is more notable (Fig. 7d–f), in which the amplitude variation can be more than 3%. Because the uncertainty of amplitude measurements are usually less than 5% in previous studies (Lin et al., 2012, 2014; Berg et al., 2018), a 3% amplitude variation caused by the inaccurate orientation can produce nonnegligible errors in amplitude measurements such as the ZH ratios. Hence, the phase‐tracking method in the orientation determination of OBS arrays maybe essential for accurate waveform‐phased measurements of phase delays and amplitudes as well as the following inversion.

We present a new workflow for determining unbiased orientations of OBSs by combining the R‐pol/P‐pol analyses and an array‐based phase‐tracking method that accounts for the off‐GCP propagation azimuths, and apply this technique to the BTF OBSs. Our results exhibit more accurate orientations with smaller uncertainties than those from previous techniques assuming the GCP‐predicted propagation azimuths. After applying the phase‐tracking method, the orientations determined from the R‐pol and P‐pol analyses are more consistent, indicating the robustness of combining the phase‐tracking method with either the R‐pol or P‐pol analyses in the orientation estimation of OBS arrays. The results are further corroborated by comparing the three‐component synthetics generated from three‐dimensional elastic full‐wave forward simulation to the observed waveforms. Moreover, we quantitatively assess the effect of biased orientations derived from the GCP‐predicted azimuths on seismic waveforms and find that the biased orientations may produce nonnegligible bias in measured phase travel times and amplitudes. This study demonstrates the promise of the new method in accurate data processing for the ongoing rapidly developed OBS array experiments. Furthermore, temporary inland seismic experiments also suffer from erroneous orientation, suggesting that this workflow can be applied to inland seismic arrays as well.

We used waveform data from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (DMC) MetaData Aggregator available at http://ds.iris.edu/mda (last accessed May 2021). The seafloor topography data were available at https://www.gmrt.org/ (last accessed November 2022). The Global Centroid Moment Tensor project database was available at https://www.globalcmt.org/CMTsearch.html (last accessed November 2022). The detailed information of the Vancouver Island event was available at https://earthquake.usgs.gov/earthquakes/eventpage/usp000jv6k/origin/detail (last accessed November 2022). All figures are plotted by Generic Mapping Tools version 6 (Wessel et al., 2019). The supplemental material includes the introduction concerning the preprocedures for OBSs data, the R‐pol and P‐pol analyses, and the processes for determining orientations and their errors.

The authors acknowledge that there are no conflicts of interest recorded.

The authors wish to thank all persons who contributed to data collection in the Blanco Transform ocean‐bottom seismometer (OBS) experiment (Nábělek and Braunmiller, 2012). The authors wish to thank Václav Kuna and John Nábělek for helpful discussions on the Blanco Transform OBS experiment, Xiaozhou Yang and Hongrui Xu for their help with data analysis, and Le Ba Manh and Ting Yang for discussions on polarity flipping. The authors would like to thank two anonymous reviewers and an associate editor for their constructive comments and suggestions. This work was supported by the National Natural Science Foundation of China (Grant Numbers 41976046, 42174063, 92155307, and 41904046), Guangdong Provincial Key Laboratory of Geophysical High‐resolution Imaging Technology (Grant Number 2022B1212010002), and the Key Special Project for introduced Talents Team of Southern Marine Science and Engineering Guangdong (Guangzhou) (Grant Number GML2019ZD0203).

Supplementary data