# Measurement of seismometer orientation using the tangential P-wave receiver function based on harmonic decomposition

Measurement of seismometer orientation using the tangential P-wave receiver function based on... Summary Accurate determination of the seismometer orientation is a prerequisite for seismic studies including, but not limited to seismic anisotropy. While borehole seismometers on land produce seismic waveform data somewhat free of human-induced noise, they might have a drawback of an uncertain orientation. This study calculates a harmonic decomposition of teleseismic receiver functions from the P and PP phases and determines the orientation of a seismometer by minimizing a constant term in a harmonic expansion of tangential receiver functions in backazimuth near and at 0 s. This method normalizes the effect of seismic sources and determines the orientation of a seismometer without having to assume for an isotropic medium. Compared to the method of minimizing the amplitudes of a mean of the tangential receiver functions near and at 0 s, the method yields more accurate orientations in cases where the backazimuthal coverage of earthquake sources (even in the case of ocean bottom seismometers) is uneven and incomplete. We apply this method to data from the Korean seismic network (52 broad-band velocity seismometers, 30 of which are borehole sensors) to estimate the sensor orientation in the period of 2005−2016. We also track temporal changes in the sensor orientation through the change in the polarity and the amplitude of the tangential receiver function. Six borehole stations are confirmed to experience a significant orientation change (10°−180°) over the period of 10 yr. We demonstrate the usefulness of our method by estimating the orientation of ocean bottom sensors, which are known to have high noise level during the relatively short deployment period. Body waves, Seismic anisotropy, Seismic instruments, Wave propagation, Crustal structure 1 INTRODUCTION Three components of a seismometer are used in most seismic studies for constraining earthquake sources and structural complexities on Earth. The orientation of horizontal components is critical for various seismic methods such as teleseismic receiver functions, studies of anisotropy, body- and surface-wave polarization and surface wave dispersion. Among these methods, the studies of seismic anisotropy in particular require accurate orientation of the horizontal components for meaningful interpretation, not only on the state of deformation in the crust and mantle in various tectonic settings but their deformation history in the past (e.g. Long & Silver 2009; Long & Becker 2010). The orientation of a seismometer can be misaligned during its installation and maintenance. Errors in the orientation of the horizontal components were previously estimated based on the polarization of body wave (Yoshizawa et al.1999; Schulte-Pelkum et al.2001) and surface wave (Laske 1995; Laske & Masters 1996; Larson 2000; Larson & Ekstrom 2002; Stachnik et al.2012; Zha et al.2013). In some cases, the sensor orientation is misaligned by more than 10° even in quality stations in Global Seismographic Network (GSN) (e.g. Larson & Ekström 2002). However, these methods typically assume that the medium beneath the station is isotropic and that any effect from structure and anisotropy beneath the station on body-wave and surface-wave polarization and arrival angle can be minimized by averaging over a large number of measurements from different backazimuth. Since the backazimuth path coverage is typically incomplete, the estimated orientation can deviate from the true sensor orientation by over 10° (Schulte-Pelkum et al.2001; Wang et al.2016). In this study, we design a new algorithm to scan and detect the change of sensor orientation by utilizing teleseismic receiver functions, which, by construction, remove the effect of source complexities and source-side structure response and isolate structure response beneath the receiver (e.g. Vinnik 1977; Langston 1979). Radial receiver function (R-RF) and tangential receiver function (T-RF) are calculated by deconvolving the radial and tangential components of teleseismic waves, respectively, from the vertical component. Through harmonic decomposition of stacked R-RF and T-RF in backazimuth gathers (e.g. Park & Levin 2016), we determine the sensor orientation angle by minimizing a constant harmonic term of the T-RF at and near 0 s. To demonstrate the utility of our proposed method, we measure the sensor orientation of surface and borehole broad-band seismometers installed at 52 sites in South Korea and compare them against previous estimates (Figs 1 and 2). We also apply this method to relatively lower signal-to-noise ratio (SNR) data from two ocean bottom seismometers (OBSs) in Cascadia Initiative (Toomey et al.2014) to discuss the usefulness of our method when the azimuthal data gap is relatively large during the temporary seismic deployment. Figure 1. View largeDownload slide Map of seismic stations. The seismic networks are operated by Korea Meteorological Administration (KMA; denoted as KS) and Korea Institute of Geoscience and Mineral Resources (KIGAM; denoted as KG). GSN denotes Global Seismograph Network. A station code, SEO&2, denotes stations SEO and SEO2 that are located 0.7 km apart from each other. Of the total 53 station locations plotted, data from the station HDB is excluded in the analysis because of unstable sensor performance (Lee & Sheen 2015). An inset shows the locations of Korean seismic network and ocean bottom sensor network of Cascadia Initiative (Toomey et al.2014). Figure 1. View largeDownload slide Map of seismic stations. The seismic networks are operated by Korea Meteorological Administration (KMA; denoted as KS) and Korea Institute of Geoscience and Mineral Resources (KIGAM; denoted as KG). GSN denotes Global Seismograph Network. A station code, SEO&2, denotes stations SEO and SEO2 that are located 0.7 km apart from each other. Of the total 53 station locations plotted, data from the station HDB is excluded in the analysis because of unstable sensor performance (Lee & Sheen 2015). An inset shows the locations of Korean seismic network and ocean bottom sensor network of Cascadia Initiative (Toomey et al.2014). Figure 2. View largeDownload slide Distribution of teleseismic earthquakes in spatial and temporal scales, recorded from Korean seismic network. (a) The numbers of earthquakes for P and PP phases are 3051 and 1651, respectively, in 2005–2016. The location of the stations (Fig. 1) is indicated by a red rectangle. (b) Cumulative coverage in backazimuth ray coverage for 6 yr, discretized in 72 bins, after November 2007. By incorporating both P and PP phases, the backazimuth coverage increases to 70, 80, 90 and 95 per cent during the station operational period of about 8, 10, 21 and 30 months, respectively. (c) Distribution of the backazimuth and slowness of the earthquakes in (b). Figure 2. View largeDownload slide Distribution of teleseismic earthquakes in spatial and temporal scales, recorded from Korean seismic network. (a) The numbers of earthquakes for P and PP phases are 3051 and 1651, respectively, in 2005–2016. The location of the stations (Fig. 1) is indicated by a red rectangle. (b) Cumulative coverage in backazimuth ray coverage for 6 yr, discretized in 72 bins, after November 2007. By incorporating both P and PP phases, the backazimuth coverage increases to 70, 80, 90 and 95 per cent during the station operational period of about 8, 10, 21 and 30 months, respectively. (c) Distribution of the backazimuth and slowness of the earthquakes in (b). 2 METHOD 2.1 Harmonic decomposition of the receiver function Typically, the amplitude of the R-RFs is backazimuth independent when the subsurface structure is flat-layered isotropic media. However, azimuthally varying arrivals are often observed on both R- and T-RFs, indicating a breakdown of P–SV to SH decoupling in the presence of dipping layer and/or anisotropy. The azimuthal (i.e. horizontal symmetry axis) anisotropy leads to a characteristic 180°-periodic backazimuthal pattern in RFs (Levin & Park 1998), whereas a dipping interface or dipping symmetry axis introduces a 360°-periodic backazimuthal pattern (Cassidy 1992). Previous studies (e.g. Girardin & Farra 1998; Farra & Vinnik 2000; Bianchi et al.2010; Vinnik et al.2012; Schulte-Pelkum & Mahan 2014; Audet 2015; Park & Levin 2016) applied the harmonic decomposition method to extract the periodicity of T-RFs in backazimuth and constrained the velocity structure and the presence of anisotropy at depths. In this study, we utilize the harmonic decomposition method to the direct teleseismic P and PP phases in the context of R-RF and T-RF to estimate the orientation of the horizontal components. 2.2 Estimation of sensor orientation Park & Levin (2016) showed theoretically that the R-RF can be decomposed into constant and sinusoidal harmonic terms, and the T-RF only to the sinusoidal harmonic terms when either an anisotropic layer or a dipping isotropic layer is present beneath the receiver. In this study, the unmodelled sensor misorientation is introduced as a constant harmonic term in the T-RFs. Following the approach of Park & Levin (2016) and their eq. (44), we theoretically lay out how to determine the sensor orientation in this section. The R- and T-RFs can be fitted by regression with five harmonic terms:   \begin{equation*}{\rm{R}}{{\rm{F}}_R} = {H_{R1}}{\rm{\ }} + {H_{R2}}\cos \theta + {H_{R3}}\sin \theta + {H_{R4}}\cos 2\theta + {H_{R5}}\sin 2\theta \end{equation*}and   \begin{eqnarray} {\rm{R}}{{\rm{F}}_T} \!=\! {H_{T1}}{\rm{\ }} \!+\! {H_{T2}}\cos \theta \!+\! {H_{T3}}\sin \theta \!+\! {H_{T4}}\cos 2\theta \!+\! {H_{T5}}\sin 2\theta \end{eqnarray} (1)where the θ is backazimuth. The first terms on the right-hand side of the equations for RFR and RFT are called a constant harmonic term and the others are called harmonic sinusoidal terms. Given N data in backazimuth, the regression is applied to solve the following set of two linear equations:   $${\bf G}{{\bf m}_R} = {{{\bf d}}_R}\quad {\rm {and}}\quad {{\bf G}} {{{\bf m}}_T} = {{{\bf d}}_T}\$$ (2)where dR and dT are data matrices, with dimension of N × M. Here, M indicates the data length of the RFs. The G is a matrix with the size of N × 5 and the kth row of G is [1 cosθk sinθk cos2θk sin2θk], where θk is the backazimuth corresponding to the kth row of dR and dT. In our analysis, N is no more than 72 because both R- and T-RFs are stacked over a 5° interval in backazimuth to improve the SNR. The harmonic terms of the R- and T-RFs are   \begin{eqnarray} {{\bf m}}_R^T &=& \left[ {\begin{array}{*{20}{c}} {{H_{R1}}}&\quad {{H_{R2}}}&\quad {{H_{R3}}}&\quad {{H_{R4}}}&\quad {{H_{R5}}} \end{array}} \right]\quad {\rm{and\ }} \nonumber\\ {{\bf m}}_T^T &=& \left[ {\begin{array}{*{20}{c}} {{H_{T1}}}&\quad {{H_{T2}}}&\quad {{H_{T3}}}&\quad {{H_{T4}}}&\quad {{H_{T5}}} \end{array}} \right]\end{eqnarray} (3)where the superscript T denotes a transpose operation. The harmonic terms of the eq. (3) can be solved by the standard least squares:   $${{{\bf m}}_R} = {\left( {{{{\bf G}}^T}{{\bf G}}} \right)^{ - 1}}\ {{{\bf G}}^T}{{{\bf d}}_R}\quad {\rm{and\ }}\quad {{{\bf m}}_T} = {\left( {{{{\bf G}}^T}{{\bf G}}} \right)^{ - 1}}\ {{{\bf G}}^T}{{{\bf d}}_T}.$$ (4)By introducing an arbitrary angle φ in clockwise direction from the north, we can now represent $${{\bf d}}_R^\prime$$ and $${{\bf d}}_T^\prime$$ as the data matrix of the misoriented sensor using the rotation matrix,   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf d}}_R^\prime}\\ {{{\bf d}}_T^\prime} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {\cos \varphi {{{\bf I}}_N}}&\quad {\sin \varphi {{{\bf I}}_N}}\\ { - \sin \varphi {{{\bf I}}_N}}&\quad {\cos \varphi {{{\bf I}}_N}} \end{array}} \right]\ \left[ {\begin{array}{@{}*{1}{c}@{}} {{{{\bf d}}_R}}\\ {{{{\bf d}}_T}} \end{array}} \right],$$ (5)where 0 and I are zero and identity matrices, respectively, and their subscripts indicate the dimension. Then, the resultant harmonic terms due to the misalignment, $${{\bf m}}_R^\prime$$ and $${{\bf m}}_T^\prime$$, can be expressed as   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf m}}_R^\prime}\\ {{{\bf m}}_T^\prime} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}}&\quad {{{\boldsymbol 0}_{{\bf 5} \times N}}}\\ {{{\boldsymbol 0}_{{\bf 5} \times N}}}&\quad {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}} \end{array}} \right]\ \left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf d}}_R^\prime}\\ {{{\bf d}}_T^\prime} \end{array}} \right].$$ (6)Combining the eqs (5) and (6), we have   \begin{eqnarray}\left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf m}}_R^\prime}\\ {{{\bf m}}_T^\prime} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}}&\quad {{{\boldsymbol 0}_{{\bf 5}\times N}}}\\ {{{\boldsymbol 0}_{{\bf 5}\times N}}}&\quad {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}} \end{array}} \right]\nonumber\\ \times \, \left[ {\begin{array}{@{}*{2}{c}@{}} {\cos \varphi {{{\bf I}}_N}}&\quad {\sin \varphi {{{\bf I}}_N}}\\ { - \sin \varphi {{{\bf I}}_N}}&\quad {\cos \varphi {{{\bf I}}_N}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{{{\bf d}}_R}}\\ {{{{\bf d}}_T}} \end{array}} \right]. \end{eqnarray} (7)By placing the rotation matrix in front, we can further decrease the dimension from 2N × 2N to 10 × 10, and eq. (7) can be written as   \begin{eqnarray*}\ \left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf m}}_R^\prime}\\ {{{\bf m}}_T^\prime} \end{array}} \right] &=& \left[ {\begin{array}{@{}*{2}{c}@{}} {\cos \varphi {{{\bf I}}_5}}&\quad {\sin \varphi {{{\bf I}}_5}}\\ { - \sin \varphi {{{\bf I}}_5}}&\quad {\cos \varphi {{{\bf I}}_5}} \end{array}} \right]\nonumber\\ &&\times \,\left[ {\begin{array}{@{}*{2}{c}@{}} {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ -1}}{{{\bf G}}^T}}&\quad {{{\boldsymbol 0}_{{\bf 5} \times N}}}\\ {{{\boldsymbol 0}_{{\bf 5} \times N}}}&\quad {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{{{\bf d}}_R}}\\ {{{{\bf d}}_T}} \end{array}} \right].\end{eqnarray*}That is,   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf m}}_R^\prime}\\ {{{\bf m}}_T^\prime} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {\cos \varphi {{{\bf I}}_5}}&\quad {\sin \varphi {{{\bf I}}_5}}\\ { - \sin \varphi {{{\bf I}}_5}}&\quad {\cos \varphi {{{\bf I}}_5}} \end{array}} \right]\ \left[ {\begin{array}{@{}*{1}{c}@{}} {{{{\bf m}}_R}}\\ {{{{\bf m}}_T}} \end{array}} \right]$$ (8)where mR and mT are the harmonic terms of RFs from the misoriented sensor, and $${{\bf m}}_R^\prime$$ and $${{\bf m}}_T^\prime$$ are the terms from the orientation correction. The eqs (7) and (8) mean that the regression and rotation are commutative. We note that solving eq. (8) is computationally more cost-effective than solving eq. (7). Using eq. (3), eq. (8) can be written as   \begin{eqnarray*} \ H_{Ri}^\prime &=& {{\bf m}}_R^\prime \left( {i,:} \right) = \cos \varphi \ {{{\bf m}}_R}\left( {i,:} \right) + \sin \varphi {{{\bf m}}_T} \left( {i,:} \right) \nonumber\\ &=& \cos \varphi \ {H_{Ri}} + \sin \varphi {H_{Ti}}\end{eqnarray*}and   \begin{eqnarray} H_{Ti}^\prime &=& {{\bf m}}_T^\prime \left( {i,:} \right) = - \sin \varphi \ {{{\bf m}}_R}\left( {i,:} \right) + \cos \varphi {{{\bf m}}_T}\left( {i,:} \right) \nonumber\\ &=& - \sin \varphi \ {H_{Ri}} + \cos \varphi {H_{Ti}}\end{eqnarray} (9)where the integer i is from 1 to 5. We then define a misfit function f (φ), which is based on the root-mean-square in a window bounded by integers M1 and M2 (1 ≤ M1 < M2 ≤ M), as   \begin{eqnarray} f\left( {\rm{\varphi }} \right)\ &=& \sqrt {\frac{1}{{{M_2} - {M_1}}}\mathop \sum \limits_{j = {M_1}}^{{M_2}} {{\left\{ {H_{Ti = 1}^\prime\left( j \right)} \right\}}^2}} {\rm{\ }} \nonumber\\ &=& \sqrt {\frac{1}{{{M_2} - {M_1}}}\mathop \sum \limits_{j = {M_1}}^{{M_2}} {{\left\{ { - \sin \varphi {{{\bf m}}_R}\left( {1,j} \right) + \cos \varphi {{{\bf m}}_T}\left( {1,j} \right)} \right\}}^2}}\nonumber\\ \end{eqnarray} (10)where φ ranges from 0° to 180°. One can determine the orientation φmin by minimizing the f (φ) through a grid-search scheme (with an increment of 0.01° in this study). The error of φmin can be estimated by bootstrapping 90 per cent random selection from a row of the data matrices dR and dT in eq. (2) without repetition. Since we stack the R- and T-RFs over a 5° interval in backazimuth, the size of randomly selected data matrices in bootstrapping is no more than 64 × M. Ambiguity between the φmin and φmin + 180° can be resolved by choosing an orientation that yields a positive polarity of a constant harmonic term in the R-RF. 3 SYNTHETIC TEST In order to examine the robustness of our proposed method, we first compute synthetic RFs (Levin & Park 1997; Frederiksen & Bostock 2000) and show results from a series of synthetic tests based on several velocity models and ranges of backazimuth data coverage. Here we select the data window of −1.0 s and 1.0 s in the T-RFs in the misfit calculation. We compare these results with those from the method of minimizing the amplitudes of a mean of the T-RFs within the data window. Also, we consider more realistic event distribution (from Korean seismic network and ocean bottom seismic network) to test the effect of slowness variation on the harmonic terms of the RFs. The event distribution from the Korean network ensures nearly complete backazimuth coverage, whereas that from the ocean bottom seismic network does not. Lastly, we examine how the level of noise and range of backazimuthal data coverage interplay in the estimation of the orientation angle. In Sections 3.1–3.3, we discuss synthetic tests performed free of noise, whereas Section 3.4 illustrates synthetic test performed with the addition of noise. 3.1 Synthetic test for six representative models with full backazimuth data coverage We perform synthetic tests from a few velocity models, which produce some notable peaks in the T-RFs near 0 s as a function of backazimuth. In this test, we assume full backazimuthal coverage of earthquakes and slowness range from 0.04 to 0.08 s km−1. Fig. 3 shows six representative layered velocity models and corresponding synthetic T-RFs in backazimuth within the data window. Anisotropy strength, its fast axis direction, and a thickness of the anisotropic layer are chosen arbitrarily in this test. We first examine the effect of a different thickness (0.5 or 3 km) of the topmost anisotropic layer with (1) a horizontal fast axis over an isotropic half space (Figs 3a and d) or with (2) a tilted (45°) fast axis over an isotropic half space (Figs 3g and j). Fig. 3(m) displays a model with two isotropic layers divided by an interface with a 10° dip, whereas Fig. 3(p) includes an additional anisotropic dipping layer on top of the isotropic medium. Also, we test the effect of a sensor misorientation of 1° using the six models. Figure 3. View largeDownload slide Six velocity structure models and corresponding synthetic T-RFs with backazimuth. The velocity model is shown in (a), (d), (g), (j), (m) and (p), and corresponding synthetic T-RFs in (b), (e), (h), (k), (n) and (q), assuming the sensor is properly oriented. The T-RFs from the sensor misoriented by 1° (clockwise rotation) are shown in (c), (f), (i), (l), (o) and (r). The thickness of a green layer (in panels m and p; above the location of a dashed line) and isotropic velocities of P- and S-waves of green and blue layers are taken from those of the upper and mid crusts in Kim et al. (2011). The north is towards the right. The RFs from the models (a, d, g and j) are calculated by the code anirec (Levin & Park 1997) and the RFs from the models (m and p) by the code raysum (Frederiksen & Bostock 2000). The angle of incidence is 25°. We note that the amplitudes of the synthetics are not scaled. The maximum amplitudes within −1 to 1 s are (b) 0.4 per cent, (c) 1.7 per cent, (e) 2.9 per cent, (f) 3.6 per cent, (h) 0.6 per cent, (i) 2.0 per cent, (k) 5.9 per cent, (l) 6.3 per cent, (n) 3.1 per cent, (o) 4.3 per cent, (q) 3.1 per cent, and (r) 4.3 per cent. Figure 3. View largeDownload slide Six velocity structure models and corresponding synthetic T-RFs with backazimuth. The velocity model is shown in (a), (d), (g), (j), (m) and (p), and corresponding synthetic T-RFs in (b), (e), (h), (k), (n) and (q), assuming the sensor is properly oriented. The T-RFs from the sensor misoriented by 1° (clockwise rotation) are shown in (c), (f), (i), (l), (o) and (r). The thickness of a green layer (in panels m and p; above the location of a dashed line) and isotropic velocities of P- and S-waves of green and blue layers are taken from those of the upper and mid crusts in Kim et al. (2011). The north is towards the right. The RFs from the models (a, d, g and j) are calculated by the code anirec (Levin & Park 1997) and the RFs from the models (m and p) by the code raysum (Frederiksen & Bostock 2000). The angle of incidence is 25°. We note that the amplitudes of the synthetics are not scaled. The maximum amplitudes within −1 to 1 s are (b) 0.4 per cent, (c) 1.7 per cent, (e) 2.9 per cent, (f) 3.6 per cent, (h) 0.6 per cent, (i) 2.0 per cent, (k) 5.9 per cent, (l) 6.3 per cent, (n) 3.1 per cent, (o) 4.3 per cent, (q) 3.1 per cent, and (r) 4.3 per cent. Four-lobed patterns with backazimuth are distinctively shown in Figs 3(b) and (e), and two-lobed patterns in Figs 3(h), (k), (n) and (q). In particular, the T-RFs from the models with the thin topmost anisotropic layer (Figs 3a, d, g and j) show distinctive polarity reversal in the vicinity of 0 s, caused by both the coupling between the P, SV and SH at the interface between the thin anisotropic and isotropic layers (Levin & Park 1998). Synthetic T-RFs from the model with an anisotropic topmost layer near 0 s, as expected, display waveform shape similar to the derivative of a Gaussian pulse (Figs 3b, e, h and k; Levin & Park 1998). Even if the sensor is misoriented by only 1°, the mean of the T-RFs from the six models is no longer zero (Figs 3c, f, i, l, o and r). Maximum amplitudes of the T-RFs within the data window in Figs 3(b), (e), (h), (k), (n) and (q) are 0.4, 2.9, 0.6, 5.9, 3.1 and 3.1 per cent (with respect to vertical P), respectively. If the sensor is misoriented by 1°, maximum amplitudes of the T-RFs in Figs 3(c), (f), (i), (l), (o) and (r) become 1.7, 3.6, 2.0, 6.3, 4.3 and 4.3 per cent (with respect to vertical P), respectively. For comparison, robust signals from the transition zone seismic discontinuities are typically observed at 2–5 per cent with respect to vertical P. This indicates that the sensitivity of T-RFs to the sensor orientation is probably on the order of 1° in this instance. 3.2 Synthetic test with non-uniform backazimuth data coverage We next test a case where the backazimuthal distribution of earthquakes is not uniform and incomplete. In this case, we set up a model with a 3 km thick topmost anisotropic layer (Fig. 3j), which produces the two-lobed pattern of T-RFs against backazimuth (Fig. 3k). Fig. 4 shows R- and T-RFs, their means, and their five harmonic terms in the case of complete and poor backazimuthal coverage. In the case of complete backazimuthal data coverage, we observe that the constant harmonic term in the T-RFs (HT1) is identical to a mean value of the T-RFs and they are both zero (Fig. 4a). In the case of incomplete backazimuthal distribution, the mean of the T-RFs is no longer zero near 0 s within the data window and can be biased in the direction where earthquakes are concentrated (Fig. 4b). We note that T-RFs with a single large data gap in backazimuth yields the largest peak in the mean of T-RFs than those with several small data gaps (e.g. Fig. 4b). However, HT1 remains zero regardless of incomplete backazimuthal distribution. See the Supporting Information (Fig. S1) for the case from the model (Fig. 3d), which shows the four-lobed pattern of T-RFs in backazimuth. Figure 4. View largeDownload slide Synthetic test results showing means and harmonic terms of both R- and T-RFs in the case of full (a) and poor (b) backazimuthal coverage of earthquakes. We use the model with the 3 km thick topmost anisotropic layer (Fig. 3j) to compute the R- and T-RFs (left panel), and their means and harmonics (right panel). The RF is sampled with an interval of 5° in backazimuth. Figure 4. View largeDownload slide Synthetic test results showing means and harmonic terms of both R- and T-RFs in the case of full (a) and poor (b) backazimuthal coverage of earthquakes. We use the model with the 3 km thick topmost anisotropic layer (Fig. 3j) to compute the R- and T-RFs (left panel), and their means and harmonics (right panel). The RF is sampled with an interval of 5° in backazimuth. 3.3 Synthetic test with the backazimuth coverage based on real earthquake distribution Using more realistic distribution of the earthquakes, we examine the effect of the variation in the slowness on the harmonic terms of the R- and T-RFs under the noise-free condition (Fig. 5). First, we consider the event distribution from the Korean seismic network, which well exceeds 90 per cent coverage of backazimuth during the period of ∼10 yr (Fig. 2). We compute synthetic R- and T-RFs, which are based on the backazimuth and slowness of the available 2430 earthquakes recorded from the station JJB (with the longest operation period) in South Korea from both P and PP waves. Then, these RFs are stacked with a bin of 5° in backazimuth. Fig. 5(a) shows that the amplitudes of the HT1 are nearly zero. Alternatively, we consider the earthquake distribution for both P and PP arrivals from the OBS J61C from the Cascadia Initiative (Toomey et al.2014), and calculate the harmonics of the RFs. We note that the backazimuth coverage for the OBS during 9 months is far less complete than that for Korean seismic network (Supporting Information Fig. S2). We still observe that the amplitudes of the HT1 are nearly zero from the OBS data (Fig. 5b). Fig. 6 shows the misfit functions f (φ) of the synthetic RFs using the model (Fig. 3j) and the earthquake distributions (e.g. Fig. 2 for the station JJB and Supporting Information Fig. S2 for the OBS J61C). Two minima of the f (φ) based on the synthetics occur at 0° even with the realistic variations in both backazimuth and slowness (Fig. 6). The ambiguity between the minima can be resolved by selecting the φ that makes the HR1 positive. Figure 5. View largeDownload slide Synthetic test results showing means and harmonic terms of both R- and T-RFs in the cases of realistic distributions of the earthquakes from (a) station JJB in South Korea (Fig. 2) and (b) OBS J61C (Cascadia Initiative). We use the model (Fig. 3j) to compute the R- and T-RFs (left panel), and their means and harmonics (right panel). The slowness is calculated with a 1-D velocity model (Kennett et al.1995). The RFs are stacked by a bin of 5° in backazimuth. Figure 5. View largeDownload slide Synthetic test results showing means and harmonic terms of both R- and T-RFs in the cases of realistic distributions of the earthquakes from (a) station JJB in South Korea (Fig. 2) and (b) OBS J61C (Cascadia Initiative). We use the model (Fig. 3j) to compute the R- and T-RFs (left panel), and their means and harmonics (right panel). The slowness is calculated with a 1-D velocity model (Kennett et al.1995). The RFs are stacked by a bin of 5° in backazimuth. Figure 6. View largeDownload slide The function f (φ) based on synthetic R- and T-RFs. The φmin with a positive sign of the mean of the HR1 is 0.0°. The variations in the backazimuth and slowness are from (a) South Korea (station JJB) and (b) Cascadia Initiative (OBS J61C). Figure 6. View largeDownload slide The function f (φ) based on synthetic R- and T-RFs. The φmin with a positive sign of the mean of the HR1 is 0.0°. The variations in the backazimuth and slowness are from (a) South Korea (station JJB) and (b) Cascadia Initiative (OBS J61C). 3.4 Case for various noise levels on synthetic data In this section, we further explore how various levels of noise and backazimuthal data coverage affect the orientation estimates. We generate synthetic seismograms using the model shown in Fig. 3(j) (Fig. 7a), considering various levels of random noise (0–200 per cent) and data gap in backazimuth (0°–300°). The 50 per cent level of noise means that the root-mean-square of the noise is close to a half of the amplitude of the vertical component. The same level of uncorrelated noise is added to all three components. A line in Fig. 7(b) roughly indicates the angle estimate deviated from the true orientation as the noise level and backazimuth coverage range vary. Two green boxes in Fig. 7(b) roughly indicate the ranges of both noise level and backazimuth of earthquakes from the Korean seismic network and the Cascadia Initiative. The upper limit of the noise level is a reciprocal of an average of the SNR in the vertical component. Figure 7. View largeDownload slide Synthetic test results showing how various noise levels and a degree of backazimuthal coverage affect the orientation estimate. (a) An example of three-component synthetic seismograms with random noise of 2.3 per cent at backazimuth of 110° where the tangential component is maximized, based on the model with the 3 km thick topmost anisotropic layer (Fig. 3j). Amplitudes of the seismograms are scaled to the maximum amplitude of the vertical-component seismogram. (b) The angle estimate (deviated from the true orientation) by varying levels of noise and backazimuth coverage range. Red dashed lines roughly indicate the uncertainty range of 0.2°–30°. Empty range in backazimuth is centered on 270° following the set up in Fig. 4(b). Green boxes are approximate ranges in cases of South Korea (SK) and the Cascadia Initiative (CI). Figure 7. View largeDownload slide Synthetic test results showing how various noise levels and a degree of backazimuthal coverage affect the orientation estimate. (a) An example of three-component synthetic seismograms with random noise of 2.3 per cent at backazimuth of 110° where the tangential component is maximized, based on the model with the 3 km thick topmost anisotropic layer (Fig. 3j). Amplitudes of the seismograms are scaled to the maximum amplitude of the vertical-component seismogram. (b) The angle estimate (deviated from the true orientation) by varying levels of noise and backazimuth coverage range. Red dashed lines roughly indicate the uncertainty range of 0.2°–30°. Empty range in backazimuth is centered on 270° following the set up in Fig. 4(b). Green boxes are approximate ranges in cases of South Korea (SK) and the Cascadia Initiative (CI). This synthetic test for the orientation estimate as a function of the noise level and the backazimuthal coverage provides useful information on determining a sensor operational period (Fig. 7b). As the operational period extends, the backazimuth data gap (horizontal axis of the plot) decreases and so is the uncertainty of the estimate, and vice versa (Fig. 7b). While the noise level is site-specific and can be decreased with filtering and applying SNR criterion, a high SNR condition can cause the poor backazimuthal coverage. 3.5 Correlation between radial and tangential harmonic terms Park & Levin (2016) analytically showed that HR2 and HR4 correlate negatively with HT3 and HT5, respectively, and HR3 and HT5 correlate positively with HT2 and HR4, respectively. The polarities of the correlations from our data are opposite to what Park & Levin (2016) presented because we define our tangential component opposite to what they used in the left-handed coordinate system in defining radial, tangential, and vertical direction. The correlations between the radial and tangential harmonic terms are observed except for the HR3, HT2, HR5 and HT4, which have nearly zero amplitudes (Fig. 4). The correlation between HR2 and HT3 is observed even in the case with the variations in slowness and backazimuth (Fig. 5). To quantify a degree of the observed correlation, we estimate a confidence range for non-randomness in the correlation between the harmonics in a time window of −2.0 to 2.0 s (Bendat & Piersol 2000). A degree of freedom is set as the length of the time window (4.0 s) times a corner frequency (2.5 and 1.5 Hz for South Korea and the Cascadia Initiative, respectively). 4 DATA ANALYSIS 4.1 Data acquisition and processing of Korean network data For the RF calculation, we collect teleseismic earthquakes of a magnitude greater than 5.5 recorded at 52 stations (which includes 30 borehole sensors) in South Korea (Fig. 1) from 2005 to 2016 (Fig. 2a). We analyse teleseismic P and PP arrivals to help mitigate uneven distribution of earthquake sources and provide a more filled backazimuthal coverage (Fig. 2a). For the P-wave RFs, the epicentral distance range of 30°–100° is chosen to avoid complex triplicated mantle P waves (less than 30°) and complication from the core-mantle boundary (distances greater than 100°). Similarly, for the PP-wave RFs, the epicentral distance range is chosen at 100°–180°. By including the PP-wave RFs, the backazimuthal gaps of P-RFs in the East Pacific Rise and the Mexico-Peru-Chile subduction zone (Fig. 2a) can be filled. The total number of earthquakes for P and PP phases is 3051 and 1651, respectively. All waveforms (neglecting the SNR) are cut to 30 s before and 180 s after P and PP arrival times before the RF calculation. The R- and T-RFs are calculated in the frequency domain with the water level of 10−2 (Langston 1979). In order to remove the high frequency noise, a Gaussian pulse with a half-width (1σ) of 2.5 Hz is convoluted with the RFs. We observe that the half-width from 1.0 to 4.0 Hz yields stable orientation estimates (Supporting Information Fig. S3). After the RF calculation, we stack both R- and T-RFs over a 5° interval in backazimuth. 4.2 Data acquisition and processing of OBS data We select earthquakes recorded from two broad-band OBSs (J61C and J39C) in the Cascadia Initiative. The OBSs J61C and J39C are deployed in deep water at 2673 m and 2656 m water depth and involve relatively lower noise level compared to those in shallow water (Lodewyk & Sumy 2015). From these two OBSs, the earthquakes with the magnitude greater than 5.5 are collected from September 2013 to June 2014. The ranges of the epicentral distance for the P and PP phases for the RF calculation are the same as those for the Korean seismic network. The numbers of earthquakes for the P and PP phases are 276 and 95, respectively (Supporting Information Fig. S2). All waveforms are cut to 50 s before and 250 s after P and PP arrival times before deconvolution. Prior to the deconvolution procedure, we apply a bandpass filter between 0.1 and 1.5 Hz, the same frequency band used by Audet (2016) for the Cascadia Initiative data and Akuhara & Mochizuki (2015) for OBSs offshore Japan. The lower cut-off frequency of 0.1 Hz is effective to remove an influence of infragravity wave. The water level in the deconvolution is set at 10−2. We stack both R- and T-RFs over a 5° interval in backazimuth. 5 RESULTS 5.1 Orientation estimates for Korean seismic network The orientations of all 52 stations in South Korea are presented in Table 1. The orientation for each station is estimated for a specific period (denoted as validation period), which is defined by a start and an end of the period when the polarity and/or the amplitude of the T-RFs in the data window shifts significantly. Alternatively, the validation period is manually set when the sensor replacement/maintenance record is available (Supporting Information Table S1). Table 1 includes a total of 65 measurements of orientation for different validation periods from two methods, which are (1) the minimization of HT1 in the data window and (2) the minimization of the amplitudes of a mean of the T-RFs. In Table 1 and Supporting Information Table S2, we include results from the method (2) for comparison. In addition, Supporting Information Table S2 shows orientation estimates based on non-stacking approach. The difference between the orientation estimates from the two methods is illustrated in Supporting Information Fig. S4, and it increases as the backazimuthal coverage decreases. However, all values are confined within and less than 1°, except for three stations that show large differences (2.3° of SEHB, 7.6° of INCN00 and 4.6° of KOHB; not shown in Supporting Information Fig. S4). We note that differences in orientations between two cases (with and without stacking) are overall small for Korean network data except for the stations INCN00, KOHB and NPR (Supporting Information Table S2). Table 1. Orientation correction estimates for the stations in South Korea. STA  Orientation  Error  %Backazimuth  Orientation  Validation period  Previously reported data (°)    (°; from harmonic decomposition)  (°; 1σ)  coverage  (°; from mean of T-RF)  From  To    BAR  5.5  0.4  94  5.5  2007–12  2015–01        BGD  359.5  0.6  100  359.6  2005–01  2015–01  359.8a  357.2b    BOSB  31.5  0.6  87  31.3  2012–08  2015–01  33.6c      BUS2  9.3  0.3  99  9.4  2005–01  2015–01        CHC2  3.8  0.6  90  3.9  2012–08  2015–01        CHJ2  10.5  0.4  99  10.6  2005–01  2015–01        CHNB  19.4  0.6  100  19.3  2005–01  2015–01  13.5a  15.3b    DACB  268.4  0.3  94  268.5  2010–08  2015–01  267.6c      DAG2  3.5  0.3  93  3.5  2010–08  2015–01        DGY2  9.6  0.4  94  9.5  2010–08  2015–01        EMSB  338.2  0.6  93  338.2  2012–12  2015–01  337.0c      EURB  104.4  0.9  86  105.1  2012–06  2014–08  106.0c      GAHB  128.5  0.3  94  128.6  2008–10  2015–01  130.7c      GKP1  5.4  0.5  97  5.4  2005–01  2015–01  4.2a  3.4b    GOCB  234.5  0.5  89  234.3  2012–08  2015–01  235.0c      GSU  146.2  0.6  99  146.5  2005–01  2015–01  142.6a  147.7c  140.8b  GWYB  120.5  0.7  90  120  2012–08  2015–01  123.3c      HALB  355.4  1.1  92  355.1  2012–08  2015–01  221.3c      HAMB  121.0  0.5  90  120.8  2012–08  2015–01  125.0c      HAWB  246.7  0.6  89  247.1  2012–08  2015–01  246.2c      HKU  8.2  0.5  96  8.3  2005–01  2015–01  9.1a  5.6b    HSB  144.7  0.8  83  144.5  2005–01  2007–01        HSB  228.5  0.4  96  228.6  2008–01  2015–01  224.3a  226.2c  221.6b  HWCB  3.3  0.3  93  3.4  2010–08  2015–01  2.4c      HWSB  60.9  0.9  96  60.7  2008–01  2012–10  59.0a  66.9c  59.0b  HWSB  –  –  –  –  2014–01  2014–06c  155.6c      IMWB  196.7  0.4  92  197.1  2012–08  2015–01  196.1c      INCN00  335.3  5.0  92  327.7  2009–10  2012–10  356d  359d    INCN00  4.6  1.1  90  4  2013–10  2016–01  358d      INCN10  359.2  0.6  94  359.2  2009–10  2012–10  355d  357d  359d  INCN10  0.1  0.4  90  359.7  2013–10  2016–01  0d      JEO2  2.7  0.4  89  3  2012–08  2015–01        JJB  237.1  0.7  94  236.5  2007–10  2012–10  235.8a  237.6b    JJB  233.8  1.2  82  233.3  2014–01  2016–01        JJU  1.9  0.6  100  1.8  2005–01  2015–01        JRB  194.7  0.3  99  194.7  2006–10  2015–01  193.7a  191.6c  192.5b  JSB  245.0  0.3  99  245.2  2006–10  2015–01  245.5a  245.3c  243.0b  KOHB  256.3  4.7  93  251.7  2009–08  2013–12  190.4c      KOHB  –  –  –  –  2014–04  2014–06c  3.5c      KSA  4.4  0.3  99  4.2  2005–01  2015–01  0.7a  1.2b    MGB  134.7  0.5  96  134.8  2005–03  2009–08  134.2a  133.2b    MGB  10.8  0.9  87  10.8  2010–11  2012–12        MGB  359.2  0.7  73  359.2  2013–06c  2014–06c  356.0c      NAWB  148.6  0.7  90  148.8  2012–08  2015–01  151.2c      NPR  2.9  1.0  100  3.2  2005–01  2015–01  4.3a  4.2b    OKCB  251.1  0.5  92  250.9  2012–08  2015–01  253.5c      OKEB  35.2  1.1  96  34.8  2011–10  2015–01  18.2c      SEHB  264.5  0.5  92  264.2  2008–11  2011–01        SEHB  3.5  1.7  59  1.2  2011–02  2012–07        SEHB  77.2  0.4  93  77.2  2012–08  2015–01  74.9c      SEO  358.1  0.3  100  358.1  2005–01  2015–01        SEO2  0.9  0.4  93  1.2  2010–08  2015–01        SES2  7.3  0.6  89  7.2  2012–08  2015–01        SHHB  326.3  0.4  92  326.2  2011–10  2015–01  325.6c      SMKB  355.0  3.8  63  355.9  2013–09c  2014–03c  2.4c      SMKB  172.6  0.9  44  172.3  2014–04c  2014–06c  165.4c      SND  0.9  0.2  99  0.9  2005–01  2015–01  359.0a  358.4b    SNU  3.8  0.6  92  3.8  2005–01  2007–09  5.5a,†  5.2b,†    SNU  –  –  –  –  2007–11  2008–01        SNU  3.5  0.3  96  3.9  2008–03  2015–01        TJN  1.6  0.4  99  1.5  2005–01  2015–01  0.1a  355.0b    ULJ2  358.0  0.3  94  358.1  2010–08  2015–01        ULL  9.7  0.9  96  9.4  2005–01  2009–12        ULLB  10.4  0.5  99  10.2  2006–10  2015–01  9.4a  20.0c  2.5b  YKB  194.9  0.3  97  194.9  2006–10  2015–01  191.0a  187.0c  195.7b  YNCB  309.0  0.7  90  309.1  2008–10  2010–08        YNCB  21.0  0.5  94  20.8  2011–01  2015–01  28.4c      YSB  13.6  0.4  94  13.8  2009–02  2015–01  13.8c      STA  Orientation  Error  %Backazimuth  Orientation  Validation period  Previously reported data (°)    (°; from harmonic decomposition)  (°; 1σ)  coverage  (°; from mean of T-RF)  From  To    BAR  5.5  0.4  94  5.5  2007–12  2015–01        BGD  359.5  0.6  100  359.6  2005–01  2015–01  359.8a  357.2b    BOSB  31.5  0.6  87  31.3  2012–08  2015–01  33.6c      BUS2  9.3  0.3  99  9.4  2005–01  2015–01        CHC2  3.8  0.6  90  3.9  2012–08  2015–01        CHJ2  10.5  0.4  99  10.6  2005–01  2015–01        CHNB  19.4  0.6  100  19.3  2005–01  2015–01  13.5a  15.3b    DACB  268.4  0.3  94  268.5  2010–08  2015–01  267.6c      DAG2  3.5  0.3  93  3.5  2010–08  2015–01        DGY2  9.6  0.4  94  9.5  2010–08  2015–01        EMSB  338.2  0.6  93  338.2  2012–12  2015–01  337.0c      EURB  104.4  0.9  86  105.1  2012–06  2014–08  106.0c      GAHB  128.5  0.3  94  128.6  2008–10  2015–01  130.7c      GKP1  5.4  0.5  97  5.4  2005–01  2015–01  4.2a  3.4b    GOCB  234.5  0.5  89  234.3  2012–08  2015–01  235.0c      GSU  146.2  0.6  99  146.5  2005–01  2015–01  142.6a  147.7c  140.8b  GWYB  120.5  0.7  90  120  2012–08  2015–01  123.3c      HALB  355.4  1.1  92  355.1  2012–08  2015–01  221.3c      HAMB  121.0  0.5  90  120.8  2012–08  2015–01  125.0c      HAWB  246.7  0.6  89  247.1  2012–08  2015–01  246.2c      HKU  8.2  0.5  96  8.3  2005–01  2015–01  9.1a  5.6b    HSB  144.7  0.8  83  144.5  2005–01  2007–01        HSB  228.5  0.4  96  228.6  2008–01  2015–01  224.3a  226.2c  221.6b  HWCB  3.3  0.3  93  3.4  2010–08  2015–01  2.4c      HWSB  60.9  0.9  96  60.7  2008–01  2012–10  59.0a  66.9c  59.0b  HWSB  –  –  –  –  2014–01  2014–06c  155.6c      IMWB  196.7  0.4  92  197.1  2012–08  2015–01  196.1c      INCN00  335.3  5.0  92  327.7  2009–10  2012–10  356d  359d    INCN00  4.6  1.1  90  4  2013–10  2016–01  358d      INCN10  359.2  0.6  94  359.2  2009–10  2012–10  355d  357d  359d  INCN10  0.1  0.4  90  359.7  2013–10  2016–01  0d      JEO2  2.7  0.4  89  3  2012–08  2015–01        JJB  237.1  0.7  94  236.5  2007–10  2012–10  235.8a  237.6b    JJB  233.8  1.2  82  233.3  2014–01  2016–01        JJU  1.9  0.6  100  1.8  2005–01  2015–01        JRB  194.7  0.3  99  194.7  2006–10  2015–01  193.7a  191.6c  192.5b  JSB  245.0  0.3  99  245.2  2006–10  2015–01  245.5a  245.3c  243.0b  KOHB  256.3  4.7  93  251.7  2009–08  2013–12  190.4c      KOHB  –  –  –  –  2014–04  2014–06c  3.5c      KSA  4.4  0.3  99  4.2  2005–01  2015–01  0.7a  1.2b    MGB  134.7  0.5  96  134.8  2005–03  2009–08  134.2a  133.2b    MGB  10.8  0.9  87  10.8  2010–11  2012–12        MGB  359.2  0.7  73  359.2  2013–06c  2014–06c  356.0c      NAWB  148.6  0.7  90  148.8  2012–08  2015–01  151.2c      NPR  2.9  1.0  100  3.2  2005–01  2015–01  4.3a  4.2b    OKCB  251.1  0.5  92  250.9  2012–08  2015–01  253.5c      OKEB  35.2  1.1  96  34.8  2011–10  2015–01  18.2c      SEHB  264.5  0.5  92  264.2  2008–11  2011–01        SEHB  3.5  1.7  59  1.2  2011–02  2012–07        SEHB  77.2  0.4  93  77.2  2012–08  2015–01  74.9c      SEO  358.1  0.3  100  358.1  2005–01  2015–01        SEO2  0.9  0.4  93  1.2  2010–08  2015–01        SES2  7.3  0.6  89  7.2  2012–08  2015–01        SHHB  326.3  0.4  92  326.2  2011–10  2015–01  325.6c      SMKB  355.0  3.8  63  355.9  2013–09c  2014–03c  2.4c      SMKB  172.6  0.9  44  172.3  2014–04c  2014–06c  165.4c      SND  0.9  0.2  99  0.9  2005–01  2015–01  359.0a  358.4b    SNU  3.8  0.6  92  3.8  2005–01  2007–09  5.5a,†  5.2b,†    SNU  –  –  –  –  2007–11  2008–01        SNU  3.5  0.3  96  3.9  2008–03  2015–01        TJN  1.6  0.4  99  1.5  2005–01  2015–01  0.1a  355.0b    ULJ2  358.0  0.3  94  358.1  2010–08  2015–01        ULL  9.7  0.9  96  9.4  2005–01  2009–12        ULLB  10.4  0.5  99  10.2  2006–10  2015–01  9.4a  20.0c  2.5b  YKB  194.9  0.3  97  194.9  2006–10  2015–01  191.0a  187.0c  195.7b  YNCB  309.0  0.7  90  309.1  2008–10  2010–08        YNCB  21.0  0.5  94  20.8  2011–01  2015–01  28.4c      YSB  13.6  0.4  94  13.8  2009–02  2015–01  13.8c      Notes: The orientation angle ranges from 0° to 360° with clockwise direction, and our values indicate the angles that we need to correct from the misaligned north. We note that we apply this correction scheme to previous estimates (Ekström & Busby 2008; Shin et al.2009; Lee & Rhie 2015; Lee & Sheen 2015). ‘00’ and ‘10’: location code given for the station INCN (http://ds.iris.edu/mda/IU/INCN) by the Data Management Center (DMC) of the Incorporated Research Institutions for Seismology (IRIS). aShin et al. (2009). bLee & Rhie (2015). cLee & Sheen (2015). dEkström & Busby (2008); http://www.ldeo.columbia.edu/∼ekstrom/Projects/WQC/COMB_QC/POL_IU_S_ALL_EPO.html. †Their periods of used data (January 2007 to September 2008) contain the replacement dates of sensors: 23 October 2007 and 18 February 2008 (Supporting Information Table S1). View Large Uncertainty estimates (1σ) of the orientation from the method (1) are shown in Table 1. The 57 out of 65 measurements have an uncertainty of less than 1° (Table 1). Table 1 also includes previous estimates for some stations (Ekström & Busby 2008; Shin et al.2009; Lee & Rhie 2015; Lee & Sheen 2015) for comparisons. We note that estimates from station SMKB (for a validation period of April 2014–June 2014) involve relatively large uncertainty value (Table 1) due to poor backazimuthal coverage (44 per cent). Lee & Sheen (2015) previously reported changes in the orientation for that station, which have not been reported by the seismic network operator. We confirm that significant temporal variations of the polarity and amplitude of the T-RFs likely had stemmed during sensor replacement/maintenance (e.g. Supporting Information Table S1). Of the 52 stations examined, we present results from three stations (HSB, MGB and SEHB) that show substantial changes in both polarity and amplitude of the T-RFs. First, we illustrate such changes by showing T-RFs stacked over an interval of two months without overlap (Supporting Information Figs S5a, c and e for stations HSB, MGB and SEHB, respectively). Second, we show the maximum and minimum amplitudes of the stacked T-RFs in the data window (Supporting Information Figs S5b, d and f for stations HSB, MGB and SEHB, respectively; Supporting Information Fig. S6 for the other stations). Our results show that six borehole stations experience sudden orientation change ranging from 12° (at station MGB) to 178° (at station SMKB) over the period of 10 yr (Table 1). Fig. 8 shows the misfit function f (φ) of stations HSB, MGB and SEHB at different periods. The data from station HSB shows that f (φ) for the period of January 2005–January 2007 is minimized at 144.7° with a plus sign for the mean of HR1 (Fig. 8a). The orientation for the period of January 2008– January 2015 is estimated as 228.5° (48.5° + 180°) because f (φ) is minimized at 48.5° with a minus sign for the mean of HR1 (Fig. 8b). The orientation estimates for station MGB are 134.7°, 10.8° and 359.2° (179.2° + 180°) for three different periods (Figs 8c–e). The estimates for station SEHB are 264.5° (84.5° + 180°) and 77.2° in two different periods (Figs 8f and g). Figure 8. View largeDownload slide The function f (φ) for stations HSB, MGB and SEHB plotted against φ for different periods. The φmin for the different periods are 144.7° and 48.5° for HSB (a and b), 134.7°, 10.8° and 179.2° for MGB (c, d, e), and 84.5° and 77.2° for SEHB (f and g). Figure 8. View largeDownload slide The function f (φ) for stations HSB, MGB and SEHB plotted against φ for different periods. The φmin for the different periods are 144.7° and 48.5° for HSB (a and b), 134.7°, 10.8° and 179.2° for MGB (c, d, e), and 84.5° and 77.2° for SEHB (f and g). We show temporal changes of f (φ) and φmin for stations HSB, MGB and SEHB with a moving window of one year (Fig. 9, background colour and solid black lines, respectively). The orientation angles for the stations and their validation periods (Table 1) are also marked as references (Fig. 9; white ticked lines). We find periods that involve a sudden change in φmin, which are 2007–2008 for HSB (Fig. 9a) and October 2009–April 2010, June 2010–December 2010 and October 2012–February 2013 for MGB (Fig. 9b). These correspond to the time when the orientation had changed, which is consistent with a record of the sensor replacement/maintenance (Supporting Information Table S1). Besides the record, we also identify other periods when either maximum or minimum amplitude of the T-RFs is statically changed (Supporting Information Figs S5b, d and f for stations HSB, MGB and SEHB, respectively). We report that the period of 2011–2012 for station SEHB involves a significant change in φmin (Fig. 9c), which is considered to be the orientation change but not noted by the seismic network operator. In addition, we observe that both f (φ) and φmin slightly deviate from the reference during the periods of August 2005–December 2006, October 2008–August 2009 and June 2013–October 2014 for HSB (Fig. 9a) and during the periods of February 2007–June 2007 and October 2011– October 2012 for MGB (Fig. 9b). Such periods moderately coincide with the temporal fluctuation of backazimuthal coverage of earthquakes (Supporting Information Figs S5a and b). Figure 9. View large Download slide The f (φ) and φmin as a function of time for stations HSB, MGB and SEHB. The f (φ) is calculated with an interval of two months and a moving window of 1 yr. The f (φ) and φmin are represented by colour and a solid black line, respectively. The validation periods for our orientation estimates (Table 1) are indicated by white ticked lines. Blue arrows indicate reported replacement date of sensors (Supporting Information Table S1). Grey background colour means no data. Figure 9. View large Download slide The f (φ) and φmin as a function of time for stations HSB, MGB and SEHB. The f (φ) is calculated with an interval of two months and a moving window of 1 yr. The f (φ) and φmin are represented by colour and a solid black line, respectively. The validation periods for our orientation estimates (Table 1) are indicated by white ticked lines. Blue arrows indicate reported replacement date of sensors (Supporting Information Table S1). Grey background colour means no data. We next show how the harmonic terms of R- and T-RFs, and also means both R- and T-RFs, for the three stations change with the orientation correction based on the method (1) (Fig. 10). The amplitudes of HT1 for the different periods are minimized and essentially nearly zero after the orientation correction (Fig. 10). We now observe strong coherence in waveform between the harmonic terms for different validation periods, and between the means after the correction. This correlation is expected unless there is any temporal change in structure and significant change in the distribution of the earthquakes. We observe that HR2 and HR4 correlate negatively with HT3 and HT5, respectively, and HR3 and HT5 correlate positively with HT2 and HR4, respectively (Fig. 10; Supporting Information Figs S7–S9). We note that the amplitudes of HR4, HR5, HT4, and HT5 are small. Thus, we do not see a clear correlation in waveform between the different harmonics (see Supporting Information Figs S7–S9). Figure 10. View large Download slide Harmonic terms of R- and T-RFs, and their means for the stations HSB, MGB and SEHB prior to (left panel) and after the correction (right panel). Note that the amplitudes of the HT1 are minimized to nearly zero by the orientation correction. See the Supporting Information (Figs S7–S9) for HR4, HR5, HT4 and HT5 for the stations. Figure 10. View large Download slide Harmonic terms of R- and T-RFs, and their means for the stations HSB, MGB and SEHB prior to (left panel) and after the correction (right panel). Note that the amplitudes of the HT1 are minimized to nearly zero by the orientation correction. See the Supporting Information (Figs S7–S9) for HR4, HR5, HT4 and HT5 for the stations. Lastly, we show T-RFs and HT1 for the selected station HSB with and without the orientation correction (Fig. 11). The similar two-lobed patterns of the T-RFs with backazimuth are revealed by the correction (Figs 11c and d). See Supporting Information Fig. S10 for the stations MGB and SEHB, respectively. Figure 11. View large Download slide Stacked T-RFs plotted according to backazimuths for station HSB for each different validation period, prior to (a and b) and after the orientation correction (c and d). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). See Supporting Information Fig. S10 for stations MGB and SEHB. Figure 11. View large Download slide Stacked T-RFs plotted according to backazimuths for station HSB for each different validation period, prior to (a and b) and after the orientation correction (c and d). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). See Supporting Information Fig. S10 for stations MGB and SEHB. 5.2 Orientation estimates for OBS network Our estimates based on the harmonic decomposition are compared against previous estimates based on the minimization of stacked T-RFs at 0.0–5.0 s (Janiszewski & Abers 2015) and surface wave arrival angles (Lodewyk & Sumy 2015; Table 2). Our estimates are similar to the estimates by Lodewyk & Sumy (2015) and the estimates by Janiszewski & Abers (2015). We note that slight deviation between the estimates might stem from the different data selection criteria and processing. Table 2. Orientation correction estimates for the stations from Cascadia Initiative (Toomey et al.2014). STA  Orientation (°; from harmonic decomposition)  Error (°; 1σ)  %Backazimuth coverage  Orientation (°; from mean of T-RF)  Previously reported data (°; 1σ)  J61C  190.2  6.6  79  185.9  200 ± 2a  187 ± 9b  J39C  89.0  4.1  79  97.1  105 ± 1a  92 ± 9b  STA  Orientation (°; from harmonic decomposition)  Error (°; 1σ)  %Backazimuth coverage  Orientation (°; from mean of T-RF)  Previously reported data (°; 1σ)  J61C  190.2  6.6  79  185.9  200 ± 2a  187 ± 9b  J39C  89.0  4.1  79  97.1  105 ± 1a  92 ± 9b  Notes: The orientation angle ranges from 0° to 360° with clockwise direction, and our values indicate the angles that we need to correct from the misaligned north. We note that we apply this correction scheme to previous estimates (Janiszewski & Abers 2015; Lodewyk & Sumy 2015). For the uncertainty estimates, Janiszewski & Abers (2015) used the 95 per cent confidence bounds from the F-test, with degrees of freedom determined from the net filter response of the signal. Our uncertainty estimates for J61C and J39 from the F-test are 3.3° and 2.5°, respectively. The difference in the two estimates for each station may arise from slightly different dataset and filtering. aJaniszewski & Abers (2015). bLodewyk & Sumy (2015). View Large We observe very small positive amplitudes of HR1 near 0 s even after the orientation corrections (Fig. 12). Such small amplitude of HR1 makes the range of f (φ) very narrow (Supporting Information Fig. S11). From the data from J61C, the orientation is estimated as 190.2° (10.2° + 180°) (Supporting Information Fig. S11a). The orientation estimate for J39C is 89.0° (Supporting Information Fig. S11b). Figure 12. View largeDownload slide Harmonic terms of R- and T-RFs, and their means for the OBSs (a) J61C and (b) J39C prior to (first two columns) and after the correction (third and fourth columns). Figure 12. View largeDownload slide Harmonic terms of R- and T-RFs, and their means for the OBSs (a) J61C and (b) J39C prior to (first two columns) and after the correction (third and fourth columns). The correlations between HR2 and HT3 and between HR3 and HT2 of the J61C are observed strong, whereas those between HR4 and HT5 and between HR5 and HT4 are weak (Fig. 12a). The correlations between the harmonic terms for J39C are not all visibly clear because of small amplitudes at and near 0 s (Fig. 12b). Fig. 13 shows the T-RFs before and after the orientation correction for the two OBSs. Figure 13. View largeDownload slide Stacked T-RFs plotted according to backazimuths for the OBSs J61C and J39C, prior to (a and b) and after the orientation correction (c and d). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). Figure 13. View largeDownload slide Stacked T-RFs plotted according to backazimuths for the OBSs J61C and J39C, prior to (a and b) and after the orientation correction (c and d). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). 6 DISCUSSION Most indirect seismic methods in determining the orientation inevitably involve inaccuracy due to uncertain medium seismic properties beneath the receiver. The medium is assumed to be isotropic and homogeneous in the horizontal direction (e.g. Ekström & Busby 2008; Shin et al.2009; Zha et al.2013; Lee & Rhie 2015; Lee & Sheen 2015; Wang et al.2016). If the anisotropy of the medium were strong, the particle motion of the P wave is no longer parallel to its propagation direction. The same applies to the Rayleigh wave. Furthermore, the minimum-time path between the source and receiver might not coincide with the shortest-distance path if the medium is laterally heterogeneous. Although our estimates based on the teleseismic RFs can also be influenced by the uncertainty in the nature of the medium, the degree of its dependency can be much less than that of the studies based on the particle motions. Our method works well in the presence of anisotropy and/or dipping structures beneath the receiver (e.g. Fig. 3). Nevertheless, to make a single measurement of orientation with an uncertainty less than 1°, we require a relatively large dataset with nearly full backazimuthal coverage of teleseismic earthquakes (e.g. Fig. 2). 6.1 Orientation estimates of the Korean seismic network Most of the orientations for stations in South Korea are measured with over 85 per cent of backazimuthal coverage (Table 1). Using both P and PP phases, the 70, 80 and 90 per cent of the backazimuthal coverage are guaranteed by over 8, 10 and 21 months of station operational periods (Figs 2b and 14) or 250, 350 and 750 earthquakes, respectively (Supporting Information Fig. S12a). Even in the case of relatively poor backazimuthal coverage at station MGB (73 per cent), our value (359.2° ± 0.7°; Table 1) still agrees reasonably well with the previous estimate (Lee & Sheen 2015; a difference of 3.2°). It is difficult to statistically compare our estimate with the previous estimate because different types of data (e.g. teleseismic body-wave phase versus regional body-wave phase versus) and different data processing (e.g. passband) are involved in determining the orientation. The measurement for station SEHB includes the smallest backazimuthal coverage (59 per cent), and our estimate for this station is 3.5° ± 1.7°. Figure 14. View largeDownload slide Global distribution of the earthquakes and the period for the 70 per cent backazimuthal coverage. (a) Distribution of the earthquakes in 14 months from November 2007 with magnitude larger than 5.5. (b) The variation of the period for the 70 per cent of the backazimuthal coverage using P and PP (distance ranging from 30° to 180°). Figure 14. View largeDownload slide Global distribution of the earthquakes and the period for the 70 per cent backazimuthal coverage. (a) Distribution of the earthquakes in 14 months from November 2007 with magnitude larger than 5.5. (b) The variation of the period for the 70 per cent of the backazimuthal coverage using P and PP (distance ranging from 30° to 180°). For Korean network data (the case with the nearly full backazimuthal coverage), the backazimuthal coverage at 20 months reaches to about 90 per cent (Fig. 15a; 750 earthquakes, Supporting Information Fig. S12a), and the orientation (φmin) converges to a steady value after 20 months (Fig. 15c). The orientation angle at 20 months slightly deviates from the best estimate by 0.3°, which is smaller than the uncertainty. We note that just over three months the measured orientation deviates from the best angle only by 6° with about 60 per cent backazimuthal coverage. Figure 15. View largeDownload slide Convergence of the orientation (φmin) and corresponding backazimuthal coverage, and correlation (in −2.0 to 2.0 s) between the harmonic terms with cumulative data (a, b, c) for the station CHJ2 in South Korea and (d, e, f) for the OBS J61C of the Cascadia Initiative. The confidence level for non-randomness (Bendat & Piersol 2000) is marked as pink colour. Figure 15. View largeDownload slide Convergence of the orientation (φmin) and corresponding backazimuthal coverage, and correlation (in −2.0 to 2.0 s) between the harmonic terms with cumulative data (a, b, c) for the station CHJ2 in South Korea and (d, e, f) for the OBS J61C of the Cascadia Initiative. The confidence level for non-randomness (Bendat & Piersol 2000) is marked as pink colour. The harmonics HR2 and HT3 and the HR3 and HT2 correlate negatively and positively, respectively (Fig. 15b). The observed correlation between the harmonic terms is quantified in terms of a confidence range for randomness in the correlation (Bendat & Piersol 2000) (Fig. 15b, pink colour). While the correlation value between the HR2 and HT3 reaches to −1 just after one month, the correlation between the HR3 and HT2 slowly converges to 0.6 (Fig. 15b). The observed slow convergence and relatively low correlation (still within the confidence range of 95 per cent) for the HR3 and HT2 are primarily due to small amplitudes of the HR3 and HT2. Similarly, the correlation between the HR4 and HT5 and HR5 and HT4 are scattered within the confidence range during the whole period of 10 yr. Assuming backazimuthal coverage of earthquakes is good, strong correlation (beyond the confidence range of 95 per cent) between the orthogonal harmonic terms may further supports the reliability of our orientation estimates (Figs 10 and 15). The orientation angle is deviated by about 5° from the final estimate during the period of 2005−mid-2006, although a degree of correlations between HR3 and HT2 and between HR5 and HT4 in 2005−mid-2006 is observed high (Fig. 15b). We observe that the harmonic terms from most stations (except for the stations in oceanic islands) in South Korea share similar patterns in waveforms within the data window near 0 s (Fig. 10), and we suspect that this is in part due to the structure. For instance, a negative pulse just before 0 s and a positive pulse just after 0 s are shown in the HR2 (conversely in the HT3) of the stations (Fig. 10). The common features of the harmonic terms might indicate an approximately southward dipping of the fast axis in the anisotropic top layer beneath South Korea, as illustrated in Figs 3(g) and (j). A separate analysis of P coda phases is required for a precise investigation on the crustal seismic structure of Korea. 6.2 Applicability of the method on OBS data One clear difference in the harmonic terms of the RFs between the land-based stations and OBSs is the amplitude of the HR1 within the data window near and at 0 s. Observed small amplitudes of the HR1 (∼0.1 per cent) in the OBS data can be a source of a problem in the orientation correction for the OBS (Fig. 12). Small amplitudes of the R-RFs typically arise from a nearly vertical incidence of the teleseismic P phase due to the presence of the sedimentary layer with low P-wave velocity (Kawakatsu & Abe 2016). The pseudo vertical incidence of P can yield a very shallow range of the f (φ), which makes it difficult to determine φmin in the presence of noise since the orientation estimate depends on the amplitude of the HR1 near 0 s. Despite the high background noise embedded in the OBS data (Webb 1998; Lin et al.2010), the RFs at OBS station J61C show clear correlation between HR2 and HT3 and between HR3 and HT2 (Fig. 12a), as expected in an anisotropic or dipping media (Park & Levin 2016). Our estimates generally agree with previous estimates based on slightly different data processing scheme (Janiszewski & Abers 2015) and different seismic phase (Lodewyk & Sumy 2015). Because of shorter operation period for the OBSs (∼9 months), the backazimuthal coverage reaches to about 80 per cent in 9 months (Fig. 15d; 380 earthquakes in Supporting Information Fig. S12b) and the orientation estimate (φmin) just begins to converge (Fig. 15f). If the noise level is assumed at 80 per cent in the OBS data, about 70 per cent of backazimuthal coverage is required to obtain the orientation estimate with 5° uncertainty (Fig. 7b). For example, this can be achieved approximately within 4−7 months of the deployment in Pacific, Atlantic, Indian Oceans, and North America and 9−14 months in southeastern Asia, South America and Africa (Fig. 14). 7 CONCLUSIONS We calculate teleseismic receiver functions from the P and PP phases using the harmonic decomposition method and determine the orientation of a seismometer by minimizing a constant term in a harmonic expansion of T-RFs in backazimuth near and at 0 s. This approach is applied to data from 52 stations in South Korea and two OBSs from the Cascadia Initiative project. Our method works well in the presence of anisotropy and/or dipping structures beneath the receiver. Also, the method yields a much more robust estimate than those obtained by minimizing the average amplitude of the T-RFs near 0 s. In particular, the method effectively determines the orientation in the case of deficient backazimuthal coverage of earthquakes (i.e. from the OBSs). The station operational period of about one year guarantees ∼80 per cent of the backazimuthal coverage with the usage of both P and PP phases from earthquakes with a magnitude over 5.5 recorded in the seismic network in South Korea. Our analysis shows that the operational period of the station of 20 months (recording about 750 earthquakes) is required for a single measurement of orientation with an uncertainty less than 1°. We keep track of the history of the orientation of a total 52 broad-band velocity seismometers for the period of 2005–2016 by detecting its period when the polarity and the amplitude of the T-RFs are largely changed. Of the 30 borehole stations in South Korea, six stations are confirmed to experience a significant orientation change (12° to ∼178°) during the 10 yr period. Our method can be applied to any broad-band data, regardless of sensor types and orientation, for seismic studies (i.e. receiver functions). ACKNOWLEDGEMENTS Two sets of seismic data from Korea Institute of Geoscience and Mineral Resources (KIGAM) and Korea Meteorological Administration (KMA) are available with permissions at http://quake.kigam.re.kr (last accessed September 2014) and http://necis.kma.go.kr (last accessed April 2015), respectively. The data of INCN (Fig. 1) can be obtained from the IRIS Data Management Center at www.iris.edu (last accessed May 2015). Earthquake catalogue was obtained at US Geological Survey (http://earthquake.usgs.gov/earthquakes/search; last accessed March 2016). Receiver function calculation was performed using processRFmatlab (https://github.com/iwbailey/processRFmatlab; see Bailey et al. (2012) for detail). Most of the figures are plotted by Generic Mapping Tools version 4.5.9 (http://www.soest.hawaii.edu/gmt; Wessel & Smith 1998). YK and HL acknowledge supports from the Korea Meteorological Administration Research and Development Program under Grant KMIPA2015-7020, and from National Research Foundation of Korea Grant funded by the Korean Government (NRF-2014S1A2A2027609). TAS acknowledges the support from the Natural Environment Research Council (NE/P001378/1). Finally, authors thank the editor and two reviewers for helpful comments which improved the manuscript. REFERENCES Akuhara T., Mochizuki K., 2015. Hydrous state of the subducting Philippine Sea plate inferred from receiver function image using onshore and offshore data, J. geophys. Res. , 120( 12), 8461– 8477. Google Scholar CrossRef Search ADS   Audet P., 2015. Layered crustal anisotropy around the San Andreas fault near Parkfield, California, J. geophys. Res. , 120( 5), 3527– 3543. Google Scholar CrossRef Search ADS   Audet P., 2016. Receiver functions using OBS data: promises and limitations from numerical modelling and examples from the Cascadia Initiative, Geophys. J. Int. , 205( 3), 1740– 1755. https://doi.org/10.1093/gji/ggw111 Google Scholar CrossRef Search ADS   Bailey I.W., Miller M.S., Liu K., Levander A., 2012. Vs and density structure beneath the Colorado Plateau constrained by gravity anomalies and joint inversions of receiver function and phase velocity data, J. geophys. Res. , 117, B02313, doi:10.1029/2011JB008522. Bendat J.S., Piersol A.G., 2000. Random Data Analysis and Measurement Procedures , 3rd edn, pp. 108– 111, Wiley. Bianchi I., Park J., Agostinetti N.P., Levin V., 2010. Mapping seismic anisotropy using harmonic decomposition of receiver functions: an application to Northern Apennines, Italy, J. geophys. Res. , 115( B12), doi:10.1029/2009JB007061. https://doi.org/10.1029/2009JB007061 Cassidy J., 1992. Numerical experiments in broadband receiver function analysis, Bull. seism. Soc. Am. , 82, 1453– 1474. Ekström G., Busby R.W., 2008. Measurements of seismometer orientation at USArray transportable array and backbone stations, Seismol. Res. Lett. , 79( 4), 554– 561. https://doi.org/10.1785/gssrl.79.4.554 Google Scholar CrossRef Search ADS   Farra V., Vinnik L., 2000. Upper mantle stratification by P and S receiver functions, Geophys. J. Int. , 141( 3), 699– 712. https://doi.org/10.1046/j.1365-246x.2000.00118.x Google Scholar CrossRef Search ADS   Frederiksen A.W., Bostock M.G., 2000. Modelling teleseismic waves in dipping anisotropic structures, Geophys. J. Int. , 141( 2), 401– 412. https://doi.org/10.1046/j.1365-246x.2000.00090.x Google Scholar CrossRef Search ADS   Girardin N., Farra V., 1998. Azimuthal anisotropy in the upper mantle from observations of P-to-S converted phases: application to southeast Australia, Geophys. J. Int. , 133( 3), 615– 629. https://doi.org/10.1046/j.1365-246X.1998.00525.x Google Scholar CrossRef Search ADS   Janiszewski H.A., Abers G.A., 2015. Imaging the plate interface in the Cascadia seismogenic zone: new constraints from offshore receiver functions, Seismol. Res. Lett. , 86, doi:10.1785/0220150104. Kawakatsu H., Abe Y., 2016. Significance of sediment reverberations on receiver functions of broadband OBS data: Comments on Olugboji et al. (2016) nature of the seismic lithosphere-asthenosphere boundary within normal oceanic mantle from high-resolution receiver functions, Geochem. Geophys. Geosyst. , 17( 8), 3488– 3492. https://doi.org/10.1002/2016GC006418 Google Scholar CrossRef Search ADS   Kennett B.L.N., Engdahl E.R., Buland R., 1995. Constraints on seismic velocities in the Earth from traveltimes, Geophys. J. Int. , 122( 1), 108– 124. https://doi.org/10.1111/j.1365-246X.1995.tb03540.x Google Scholar CrossRef Search ADS   Kim S., Rhie J., Kim G., 2011. Forward waveform modelling procedure for 1-D crustal velocity structure and its application to the southern Korean Peninsula, Geophys. J. Int. , 185( 1), 453– 468. https://doi.org/10.1111/j.1365-246X.2011.04949.x Google Scholar CrossRef Search ADS   Langston C.A., 1979. Structure under Mount Rainier, Washington, Inferred from Teleseismic Body Waves, J. geophys. Res. , 84( B9), 4749–4762. https://doi.org/10.1029/JB084iB09p04749 Larson E.W., 2000. Measuring Refraction and Modeling Velocities of Surface Waves, PhD dissertation , Harvard University, Cambridge, Massachusetts. Larson E.W., Ekström G., 2002. Determining surface wave arrival angle anomalies, J. geophys. Res. , 107( B6), doi:10.1029/2000JB000048. https://doi.org/10.1029/2000JB000048 Laske G., 1995. Global observation of off-great-circle propagation of long-period surface waves, Geophys. J. Int. , 123( 1), 245– 259. https://doi.org/10.1111/j.1365-246X.1995.tb06673.x Google Scholar CrossRef Search ADS   Laske G., Masters G., 1996. Constraints on global phase velocity maps from long-period polarization data, J. geophys. Res. , 101( B7), 16 059–16 075. https://doi.org/10.1029/96JB00526 Google Scholar CrossRef Search ADS   Lee H., Sheen D., 2015. A study on determination of orientation of borehole seismometer, J. Geol. Soc. Korea , 51( 1), 93– 103 (in Korean). https://doi.org/10.14770/jgsk.2015.51.1.93 Google Scholar CrossRef Search ADS   Lee S.-J., Rhie J., 2015. Determining the orientations of broadband stations in South Korea using ambient noise cross-correlation, Geophys. Geophys. Explor. , 18( 2), 85– 90 (in Korean) https://doi.org/10.7582/GGE.2015.18.2.085 Google Scholar CrossRef Search ADS   Levin V., Park J., 1997. P-SH conversions in a flat-layered medium with anisotropy of arbitrary orientation, Geophys. J. Int. , 131( 2), 253– 266. https://doi.org/10.1111/j.1365-246X.1997.tb01220.x Google Scholar CrossRef Search ADS   Levin V., Park J., 1998. P-SH conversions in layered media with hexagonally symmetric anisotropy: a cookbook, in Geodynamics of Lithosphere & Earth's Mantle , pp. 669– 697, Springer, doi:10.1007/s000240050136. Lin C.R., Kuo B.Y., Liang W.T., Chi W.C., Huang Y.C., Collins J., Wang C.Y., 2010. Ambient noise and teleseismic signals recorded by ocean-bottom seismometers offshore eastern Taiwan, Terr. Atmos. Ocean. Sci. , 21( 5), 743– 755. https://doi.org/10.3319/TAO.2009.09.14.01(T) Google Scholar CrossRef Search ADS   Lodewyk J., Sumy D., 2015. Cascadia Amphibious Array Ocean Bottom Seismograph Horizontal Component Orientations, 2013–2014 OBS Deployments, OBSIP Management Office Incorporated Research Institutions for Seismology. Long M.D., Becker T.W., 2010. Mantle dynamics and seismic anisotropy, Earth planet. Sci. Lett. , 297( 3–4), 341– 354. https://doi.org/10.1016/j.epsl.2010.06.036 Google Scholar CrossRef Search ADS   Long M.D., Silver P.G., 2009. Shear wave splitting and mantle anisotropy: measurements, interpretations, and new directions, Surv. Geophys. , 30( 4–5), 407– 461. https://doi.org/10.1007/s10712-009-9075-1 Google Scholar CrossRef Search ADS   Park J., Levin V., 2016. Anisotropic shear zones revealed by backazimuthal harmonics of teleseismic receiver functions, Geophys. J. Int. , 207( 2), 1216– 1243. https://doi.org/10.1093/gji/ggw323 Google Scholar CrossRef Search ADS   Schulte-Pelkum V., Mahan K.H., 2014. A method for mapping crustal deformation and anisotropy with receiver functions and first results from USArray, Earth planet. Sci. Lett. , 402, 221– 233. https://doi.org/10.1016/j.epsl.2014.01.050 Google Scholar CrossRef Search ADS   Schulte-Pelkum V., Masters G., Shearer P.M., 2001. Upper mantle anisotropy from long-period P-polarization, J. geophys. Res. , 106( B10), 21 917–21 934. https://doi.org/10.1029/2001JB000346 Google Scholar CrossRef Search ADS   Shin J.S., Sheen D.-H., Shin I.C., 2009. Orientation correction for borehole seismic stations in South Korea, J. Geol. Soc. Korea , 45, 47– 54 (in Korean). Stachnik J.C., Sheehan A.F., Zietlow D.W., Yang Z., Collins J., Ferris A., 2012. Determination of New Zealand ocean bottom seismometer orientation via Rayleigh-wave polarization, Seismol. Res. Lett. , 83( 4), 704– 713. https://doi.org/10.1785/0220110128 Google Scholar CrossRef Search ADS   Toomey D.R. et al.  , 2014. The Cascadia Initiative: a sea change in seismological studies of subduction zones, Oceanography , 27( 2), 138– 150. https://doi.org/10.5670/oceanog.2014.49 Google Scholar CrossRef Search ADS   Vinnik L., 1977. Detection of waves converted from P to SV in the mantle, Phys. Earth planet. Inter. , 15( 1), 39– 45. https://doi.org/10.1016/0031-9201(77)90008-5 Google Scholar CrossRef Search ADS   Vinnik L., Kiselev S., Weber M., Oreshin S., Makeyeva L., 2012. Frozen and active seismic anisotropy beneath southern Africa, Geophys. Res. Lett. , 39( 8), doi:10.1029/2012GL051326. https://doi.org/10.1029/2012GL051326 Wang X., Chen Q.F., Li J., Wei S., 2016. Seismic sensor misorientation measurement using P-wave particle motion: an application to the NECsaids Array, Seismol. Res. Lett. , 87( 4), 901– 911. https://doi.org/10.1785/0220160005 Google Scholar CrossRef Search ADS   Webb S.C., 1998. Broadband seismology and noise under the ocean, Rev. Geophys. , 36( 1), 105– 142. https://doi.org/10.1029/97RG02287 Google Scholar CrossRef Search ADS   Wessel P., Smith W.H.F., 1998. New, improved version of the generic mapping tools released, EOS, Trans. Am. geophys. Un. , 79, 579. Google Scholar CrossRef Search ADS   Yoshizawa K., Yomogida K., Tsuboi S., 1999. Resolving power of surface wave polarization data for higher-order heterogeneities, Geophys. J. Int. , 138, 205– 220. https://doi.org/10.1046/j.1365-246x.1999.00861.x Google Scholar CrossRef Search ADS   Zha Y., Webb S.C., Menke W., 2013. Determining the orientations of ocean bottom seismometers using ambient noise correlation, Geophys. Res. Lett. , 40, 3585– 3590. https://doi.org/10.1002/grl.50698 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Table S1. Replacement date of velocity sensor of KG network in South Korea. Table S2. Comparison of the orientations based on different methods. Figure S1. Synthetic test results showing means and five harmonic terms of both R- and T-RFs in the case of full (a) and poor (b) backazimuthal coverage of earthquakes. We use the model with the 3 km thick topmost anisotropic layer (Fig. 3d) to compute the R- and T-RFs (left panel), and their means and harmonic terms (right panel). The RF is sampled with an interval of 5° in backazimuth. Figure S2. Distribution of teleseismic earthquakes in spatial and temporal scales, recorded from the OBS J61C (Cascadia Initiative; Toomey et al. 2014). (a) The numbers of earthquakes for P and PP phases are 276 and 95, respectively, in 2013–2014. The location of the station is indicated by a red rectangle. (b) Cumulative coverage in backazimuth ray coverage for nine months, discretized in 72 bins, after September 2013. (c) Distribution of the backazimuth and slowness of the earthquake in (b). Figure S3. Orientation estimations using various half-widths in Gaussian low-pass filter for the station HSB. (a) Period of January 2005−January 2007 and (b) period of January 2008−January 2015. The infinity symbol means no filter. Figure S4. Difference between the orientation estimates from two methods, which are (1) the minimization of HT1 and (2) the minimization of the mean of T-RF near 0 s (Table S2). The stations SEHB, INCN00 and KOHB are excluded in this figure because of large differences (2.3°, 7.6° and 4.6°, respectively) in this figure. Figure S5. Stacked T-RFs and temporal change of their maximum and minimum amplitudes for stations HSB, MGB and SEHB. (a, c, e) Stacked T-RFs (lower panel) and a histogram showing the number of earthquakes in each bin with a two-month window (upper panel) for stations HSB, MGB and SEHB. Blue arrows indicate reported replacement date of sensors (Table S1). (b, d, f) Temporal change of the maximum and minimum amplitudes (black and blue lines, respectively) of the T-RFs within −1.2 to 1.2 s for stations HSB, MGB and SEHB. Coloured background indicates the number of earthquakes in time and backazimuth. The fluctuation of either maximum or minimum amplitude can be partially explained by heterogeneous space–time occurrence of earthquakes. Grey background colour means no data. Figure S6. Temporal change of the maximum and minimum amplitudes (black and blue lines, respectively) of the T-RF within −1.2 to 1.2 s for 51 stations (including the station INCN ‘00’ and ‘10’). See Figs S5b, d and f for stations HSB, MGB and SEHB. Coloured background indicates the number of earthquakes in time and backazimuth. A blue arrow indicates reported replacement date of a sensor (Table S1). The fluctuation of either maximum or minimum amplitude can be partially explained by heterogeneous space–time occurrence of earthquakes. Grey background colour means no data. Figure S7. Harmonic terms of R- and T-RFs, and their means for the station HSB during the two different validation periods. In each panel, first two columns show the results prior to the correction, and the third and fourth columns after the correction. Figure S8. Harmonic terms of R- and T-RFs, and their means for the station MGB during the three different validation periods. In each panel, first two columns show the results prior to the correction, and the third and fourth columns after the correction. Figure S9. Harmonic terms of R- and T-RFs, and their means for the station SEHB during the two different validation periods. In each panel, first two columns show the results prior to the correction, and the third and fourth columns after the correction. Figure S10. Stacked T-RFs plotted according to backazimuths for stations HSB, MGB and SEHB for each different validation period, prior to (a, c, e, g, i, k, m) and after the orientation correction (b, d, f, h, j, l, n). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). Figure S11. The function f (φ) for OBSs (a) J61C and (b) J39C plotted against φ. Note that the range of the function f (φ) is quite limited because of small amplitudes of the T-RFs at and near 0 s for the OBS data. Regardless of such small range in f (φ), the φmin can be determined as shown in the figure. Figure S12. Backazimuthal coverage plotted against the cumulative numbers of earthquakes for (a) the station JJB and (b) OBS J61C. Figure S13. Histogram showing our uncertainties and differences in the orientation estimates shown in Table 1. (a) The error based on the bootstrapping method (the third column in Table 1). (b) Difference between the orientations determined by the minimization of HT1 (the second column) and the mean of T-RFs (the fifth column) (see also Fig. S4). Difference between ours (the second column) and previous estimates, (c) Shin et al. (2009) (the eighth column), (d) Lee & Sheen (2015) (the ninth column), (e) Lee & Rhie (2015) (the tenth column), respectively. The differences shown in (b), (c), (d) and (e) are root-mean-square values, and we exclude a few values which exceed 10°. Dif.—difference; col.—column; stddev—standard deviation. Figure S14. Synthetic test results showing means and harmonic terms of both R- and T-RFs in the cases of (a and b) for Korean seismic network and (c and d) the Cascadia Initiative when the RFs are stacked by 5° of bin in backazimuth (a and c) and not stacked (b and d). Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Measurement of seismometer orientation using the tangential P-wave receiver function based on harmonic decomposition

, Volume 212 (3) – Mar 1, 2018
19 pages

/lp/ou_press/measurement-of-seismometer-orientation-using-the-tangential-p-wave-9PQn0FbWhZ
Publisher
The Royal Astronomical Society
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx515
Publisher site
See Article on Publisher Site

### Abstract

Summary Accurate determination of the seismometer orientation is a prerequisite for seismic studies including, but not limited to seismic anisotropy. While borehole seismometers on land produce seismic waveform data somewhat free of human-induced noise, they might have a drawback of an uncertain orientation. This study calculates a harmonic decomposition of teleseismic receiver functions from the P and PP phases and determines the orientation of a seismometer by minimizing a constant term in a harmonic expansion of tangential receiver functions in backazimuth near and at 0 s. This method normalizes the effect of seismic sources and determines the orientation of a seismometer without having to assume for an isotropic medium. Compared to the method of minimizing the amplitudes of a mean of the tangential receiver functions near and at 0 s, the method yields more accurate orientations in cases where the backazimuthal coverage of earthquake sources (even in the case of ocean bottom seismometers) is uneven and incomplete. We apply this method to data from the Korean seismic network (52 broad-band velocity seismometers, 30 of which are borehole sensors) to estimate the sensor orientation in the period of 2005−2016. We also track temporal changes in the sensor orientation through the change in the polarity and the amplitude of the tangential receiver function. Six borehole stations are confirmed to experience a significant orientation change (10°−180°) over the period of 10 yr. We demonstrate the usefulness of our method by estimating the orientation of ocean bottom sensors, which are known to have high noise level during the relatively short deployment period. Body waves, Seismic anisotropy, Seismic instruments, Wave propagation, Crustal structure 1 INTRODUCTION Three components of a seismometer are used in most seismic studies for constraining earthquake sources and structural complexities on Earth. The orientation of horizontal components is critical for various seismic methods such as teleseismic receiver functions, studies of anisotropy, body- and surface-wave polarization and surface wave dispersion. Among these methods, the studies of seismic anisotropy in particular require accurate orientation of the horizontal components for meaningful interpretation, not only on the state of deformation in the crust and mantle in various tectonic settings but their deformation history in the past (e.g. Long & Silver 2009; Long & Becker 2010). The orientation of a seismometer can be misaligned during its installation and maintenance. Errors in the orientation of the horizontal components were previously estimated based on the polarization of body wave (Yoshizawa et al.1999; Schulte-Pelkum et al.2001) and surface wave (Laske 1995; Laske & Masters 1996; Larson 2000; Larson & Ekstrom 2002; Stachnik et al.2012; Zha et al.2013). In some cases, the sensor orientation is misaligned by more than 10° even in quality stations in Global Seismographic Network (GSN) (e.g. Larson & Ekström 2002). However, these methods typically assume that the medium beneath the station is isotropic and that any effect from structure and anisotropy beneath the station on body-wave and surface-wave polarization and arrival angle can be minimized by averaging over a large number of measurements from different backazimuth. Since the backazimuth path coverage is typically incomplete, the estimated orientation can deviate from the true sensor orientation by over 10° (Schulte-Pelkum et al.2001; Wang et al.2016). In this study, we design a new algorithm to scan and detect the change of sensor orientation by utilizing teleseismic receiver functions, which, by construction, remove the effect of source complexities and source-side structure response and isolate structure response beneath the receiver (e.g. Vinnik 1977; Langston 1979). Radial receiver function (R-RF) and tangential receiver function (T-RF) are calculated by deconvolving the radial and tangential components of teleseismic waves, respectively, from the vertical component. Through harmonic decomposition of stacked R-RF and T-RF in backazimuth gathers (e.g. Park & Levin 2016), we determine the sensor orientation angle by minimizing a constant harmonic term of the T-RF at and near 0 s. To demonstrate the utility of our proposed method, we measure the sensor orientation of surface and borehole broad-band seismometers installed at 52 sites in South Korea and compare them against previous estimates (Figs 1 and 2). We also apply this method to relatively lower signal-to-noise ratio (SNR) data from two ocean bottom seismometers (OBSs) in Cascadia Initiative (Toomey et al.2014) to discuss the usefulness of our method when the azimuthal data gap is relatively large during the temporary seismic deployment. Figure 1. View largeDownload slide Map of seismic stations. The seismic networks are operated by Korea Meteorological Administration (KMA; denoted as KS) and Korea Institute of Geoscience and Mineral Resources (KIGAM; denoted as KG). GSN denotes Global Seismograph Network. A station code, SEO&2, denotes stations SEO and SEO2 that are located 0.7 km apart from each other. Of the total 53 station locations plotted, data from the station HDB is excluded in the analysis because of unstable sensor performance (Lee & Sheen 2015). An inset shows the locations of Korean seismic network and ocean bottom sensor network of Cascadia Initiative (Toomey et al.2014). Figure 1. View largeDownload slide Map of seismic stations. The seismic networks are operated by Korea Meteorological Administration (KMA; denoted as KS) and Korea Institute of Geoscience and Mineral Resources (KIGAM; denoted as KG). GSN denotes Global Seismograph Network. A station code, SEO&2, denotes stations SEO and SEO2 that are located 0.7 km apart from each other. Of the total 53 station locations plotted, data from the station HDB is excluded in the analysis because of unstable sensor performance (Lee & Sheen 2015). An inset shows the locations of Korean seismic network and ocean bottom sensor network of Cascadia Initiative (Toomey et al.2014). Figure 2. View largeDownload slide Distribution of teleseismic earthquakes in spatial and temporal scales, recorded from Korean seismic network. (a) The numbers of earthquakes for P and PP phases are 3051 and 1651, respectively, in 2005–2016. The location of the stations (Fig. 1) is indicated by a red rectangle. (b) Cumulative coverage in backazimuth ray coverage for 6 yr, discretized in 72 bins, after November 2007. By incorporating both P and PP phases, the backazimuth coverage increases to 70, 80, 90 and 95 per cent during the station operational period of about 8, 10, 21 and 30 months, respectively. (c) Distribution of the backazimuth and slowness of the earthquakes in (b). Figure 2. View largeDownload slide Distribution of teleseismic earthquakes in spatial and temporal scales, recorded from Korean seismic network. (a) The numbers of earthquakes for P and PP phases are 3051 and 1651, respectively, in 2005–2016. The location of the stations (Fig. 1) is indicated by a red rectangle. (b) Cumulative coverage in backazimuth ray coverage for 6 yr, discretized in 72 bins, after November 2007. By incorporating both P and PP phases, the backazimuth coverage increases to 70, 80, 90 and 95 per cent during the station operational period of about 8, 10, 21 and 30 months, respectively. (c) Distribution of the backazimuth and slowness of the earthquakes in (b). 2 METHOD 2.1 Harmonic decomposition of the receiver function Typically, the amplitude of the R-RFs is backazimuth independent when the subsurface structure is flat-layered isotropic media. However, azimuthally varying arrivals are often observed on both R- and T-RFs, indicating a breakdown of P–SV to SH decoupling in the presence of dipping layer and/or anisotropy. The azimuthal (i.e. horizontal symmetry axis) anisotropy leads to a characteristic 180°-periodic backazimuthal pattern in RFs (Levin & Park 1998), whereas a dipping interface or dipping symmetry axis introduces a 360°-periodic backazimuthal pattern (Cassidy 1992). Previous studies (e.g. Girardin & Farra 1998; Farra & Vinnik 2000; Bianchi et al.2010; Vinnik et al.2012; Schulte-Pelkum & Mahan 2014; Audet 2015; Park & Levin 2016) applied the harmonic decomposition method to extract the periodicity of T-RFs in backazimuth and constrained the velocity structure and the presence of anisotropy at depths. In this study, we utilize the harmonic decomposition method to the direct teleseismic P and PP phases in the context of R-RF and T-RF to estimate the orientation of the horizontal components. 2.2 Estimation of sensor orientation Park & Levin (2016) showed theoretically that the R-RF can be decomposed into constant and sinusoidal harmonic terms, and the T-RF only to the sinusoidal harmonic terms when either an anisotropic layer or a dipping isotropic layer is present beneath the receiver. In this study, the unmodelled sensor misorientation is introduced as a constant harmonic term in the T-RFs. Following the approach of Park & Levin (2016) and their eq. (44), we theoretically lay out how to determine the sensor orientation in this section. The R- and T-RFs can be fitted by regression with five harmonic terms:   \begin{equation*}{\rm{R}}{{\rm{F}}_R} = {H_{R1}}{\rm{\ }} + {H_{R2}}\cos \theta + {H_{R3}}\sin \theta + {H_{R4}}\cos 2\theta + {H_{R5}}\sin 2\theta \end{equation*}and   \begin{eqnarray} {\rm{R}}{{\rm{F}}_T} \!=\! {H_{T1}}{\rm{\ }} \!+\! {H_{T2}}\cos \theta \!+\! {H_{T3}}\sin \theta \!+\! {H_{T4}}\cos 2\theta \!+\! {H_{T5}}\sin 2\theta \end{eqnarray} (1)where the θ is backazimuth. The first terms on the right-hand side of the equations for RFR and RFT are called a constant harmonic term and the others are called harmonic sinusoidal terms. Given N data in backazimuth, the regression is applied to solve the following set of two linear equations:   $${\bf G}{{\bf m}_R} = {{{\bf d}}_R}\quad {\rm {and}}\quad {{\bf G}} {{{\bf m}}_T} = {{{\bf d}}_T}\$$ (2)where dR and dT are data matrices, with dimension of N × M. Here, M indicates the data length of the RFs. The G is a matrix with the size of N × 5 and the kth row of G is [1 cosθk sinθk cos2θk sin2θk], where θk is the backazimuth corresponding to the kth row of dR and dT. In our analysis, N is no more than 72 because both R- and T-RFs are stacked over a 5° interval in backazimuth to improve the SNR. The harmonic terms of the R- and T-RFs are   \begin{eqnarray} {{\bf m}}_R^T &=& \left[ {\begin{array}{*{20}{c}} {{H_{R1}}}&\quad {{H_{R2}}}&\quad {{H_{R3}}}&\quad {{H_{R4}}}&\quad {{H_{R5}}} \end{array}} \right]\quad {\rm{and\ }} \nonumber\\ {{\bf m}}_T^T &=& \left[ {\begin{array}{*{20}{c}} {{H_{T1}}}&\quad {{H_{T2}}}&\quad {{H_{T3}}}&\quad {{H_{T4}}}&\quad {{H_{T5}}} \end{array}} \right]\end{eqnarray} (3)where the superscript T denotes a transpose operation. The harmonic terms of the eq. (3) can be solved by the standard least squares:   $${{{\bf m}}_R} = {\left( {{{{\bf G}}^T}{{\bf G}}} \right)^{ - 1}}\ {{{\bf G}}^T}{{{\bf d}}_R}\quad {\rm{and\ }}\quad {{{\bf m}}_T} = {\left( {{{{\bf G}}^T}{{\bf G}}} \right)^{ - 1}}\ {{{\bf G}}^T}{{{\bf d}}_T}.$$ (4)By introducing an arbitrary angle φ in clockwise direction from the north, we can now represent $${{\bf d}}_R^\prime$$ and $${{\bf d}}_T^\prime$$ as the data matrix of the misoriented sensor using the rotation matrix,   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf d}}_R^\prime}\\ {{{\bf d}}_T^\prime} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {\cos \varphi {{{\bf I}}_N}}&\quad {\sin \varphi {{{\bf I}}_N}}\\ { - \sin \varphi {{{\bf I}}_N}}&\quad {\cos \varphi {{{\bf I}}_N}} \end{array}} \right]\ \left[ {\begin{array}{@{}*{1}{c}@{}} {{{{\bf d}}_R}}\\ {{{{\bf d}}_T}} \end{array}} \right],$$ (5)where 0 and I are zero and identity matrices, respectively, and their subscripts indicate the dimension. Then, the resultant harmonic terms due to the misalignment, $${{\bf m}}_R^\prime$$ and $${{\bf m}}_T^\prime$$, can be expressed as   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf m}}_R^\prime}\\ {{{\bf m}}_T^\prime} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}}&\quad {{{\boldsymbol 0}_{{\bf 5} \times N}}}\\ {{{\boldsymbol 0}_{{\bf 5} \times N}}}&\quad {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}} \end{array}} \right]\ \left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf d}}_R^\prime}\\ {{{\bf d}}_T^\prime} \end{array}} \right].$$ (6)Combining the eqs (5) and (6), we have   \begin{eqnarray}\left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf m}}_R^\prime}\\ {{{\bf m}}_T^\prime} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}}&\quad {{{\boldsymbol 0}_{{\bf 5}\times N}}}\\ {{{\boldsymbol 0}_{{\bf 5}\times N}}}&\quad {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}} \end{array}} \right]\nonumber\\ \times \, \left[ {\begin{array}{@{}*{2}{c}@{}} {\cos \varphi {{{\bf I}}_N}}&\quad {\sin \varphi {{{\bf I}}_N}}\\ { - \sin \varphi {{{\bf I}}_N}}&\quad {\cos \varphi {{{\bf I}}_N}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{{{\bf d}}_R}}\\ {{{{\bf d}}_T}} \end{array}} \right]. \end{eqnarray} (7)By placing the rotation matrix in front, we can further decrease the dimension from 2N × 2N to 10 × 10, and eq. (7) can be written as   \begin{eqnarray*}\ \left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf m}}_R^\prime}\\ {{{\bf m}}_T^\prime} \end{array}} \right] &=& \left[ {\begin{array}{@{}*{2}{c}@{}} {\cos \varphi {{{\bf I}}_5}}&\quad {\sin \varphi {{{\bf I}}_5}}\\ { - \sin \varphi {{{\bf I}}_5}}&\quad {\cos \varphi {{{\bf I}}_5}} \end{array}} \right]\nonumber\\ &&\times \,\left[ {\begin{array}{@{}*{2}{c}@{}} {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ -1}}{{{\bf G}}^T}}&\quad {{{\boldsymbol 0}_{{\bf 5} \times N}}}\\ {{{\boldsymbol 0}_{{\bf 5} \times N}}}&\quad {{{\left( {{{{\bf G}}^T}{{\bf G}}} \right)}^{ - 1}}{{{\bf G}}^T}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{{{\bf d}}_R}}\\ {{{{\bf d}}_T}} \end{array}} \right].\end{eqnarray*}That is,   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{{\bf m}}_R^\prime}\\ {{{\bf m}}_T^\prime} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {\cos \varphi {{{\bf I}}_5}}&\quad {\sin \varphi {{{\bf I}}_5}}\\ { - \sin \varphi {{{\bf I}}_5}}&\quad {\cos \varphi {{{\bf I}}_5}} \end{array}} \right]\ \left[ {\begin{array}{@{}*{1}{c}@{}} {{{{\bf m}}_R}}\\ {{{{\bf m}}_T}} \end{array}} \right]$$ (8)where mR and mT are the harmonic terms of RFs from the misoriented sensor, and $${{\bf m}}_R^\prime$$ and $${{\bf m}}_T^\prime$$ are the terms from the orientation correction. The eqs (7) and (8) mean that the regression and rotation are commutative. We note that solving eq. (8) is computationally more cost-effective than solving eq. (7). Using eq. (3), eq. (8) can be written as   \begin{eqnarray*} \ H_{Ri}^\prime &=& {{\bf m}}_R^\prime \left( {i,:} \right) = \cos \varphi \ {{{\bf m}}_R}\left( {i,:} \right) + \sin \varphi {{{\bf m}}_T} \left( {i,:} \right) \nonumber\\ &=& \cos \varphi \ {H_{Ri}} + \sin \varphi {H_{Ti}}\end{eqnarray*}and   \begin{eqnarray} H_{Ti}^\prime &=& {{\bf m}}_T^\prime \left( {i,:} \right) = - \sin \varphi \ {{{\bf m}}_R}\left( {i,:} \right) + \cos \varphi {{{\bf m}}_T}\left( {i,:} \right) \nonumber\\ &=& - \sin \varphi \ {H_{Ri}} + \cos \varphi {H_{Ti}}\end{eqnarray} (9)where the integer i is from 1 to 5. We then define a misfit function f (φ), which is based on the root-mean-square in a window bounded by integers M1 and M2 (1 ≤ M1 < M2 ≤ M), as   \begin{eqnarray} f\left( {\rm{\varphi }} \right)\ &=& \sqrt {\frac{1}{{{M_2} - {M_1}}}\mathop \sum \limits_{j = {M_1}}^{{M_2}} {{\left\{ {H_{Ti = 1}^\prime\left( j \right)} \right\}}^2}} {\rm{\ }} \nonumber\\ &=& \sqrt {\frac{1}{{{M_2} - {M_1}}}\mathop \sum \limits_{j = {M_1}}^{{M_2}} {{\left\{ { - \sin \varphi {{{\bf m}}_R}\left( {1,j} \right) + \cos \varphi {{{\bf m}}_T}\left( {1,j} \right)} \right\}}^2}}\nonumber\\ \end{eqnarray} (10)where φ ranges from 0° to 180°. One can determine the orientation φmin by minimizing the f (φ) through a grid-search scheme (with an increment of 0.01° in this study). The error of φmin can be estimated by bootstrapping 90 per cent random selection from a row of the data matrices dR and dT in eq. (2) without repetition. Since we stack the R- and T-RFs over a 5° interval in backazimuth, the size of randomly selected data matrices in bootstrapping is no more than 64 × M. Ambiguity between the φmin and φmin + 180° can be resolved by choosing an orientation that yields a positive polarity of a constant harmonic term in the R-RF. 3 SYNTHETIC TEST In order to examine the robustness of our proposed method, we first compute synthetic RFs (Levin & Park 1997; Frederiksen & Bostock 2000) and show results from a series of synthetic tests based on several velocity models and ranges of backazimuth data coverage. Here we select the data window of −1.0 s and 1.0 s in the T-RFs in the misfit calculation. We compare these results with those from the method of minimizing the amplitudes of a mean of the T-RFs within the data window. Also, we consider more realistic event distribution (from Korean seismic network and ocean bottom seismic network) to test the effect of slowness variation on the harmonic terms of the RFs. The event distribution from the Korean network ensures nearly complete backazimuth coverage, whereas that from the ocean bottom seismic network does not. Lastly, we examine how the level of noise and range of backazimuthal data coverage interplay in the estimation of the orientation angle. In Sections 3.1–3.3, we discuss synthetic tests performed free of noise, whereas Section 3.4 illustrates synthetic test performed with the addition of noise. 3.1 Synthetic test for six representative models with full backazimuth data coverage We perform synthetic tests from a few velocity models, which produce some notable peaks in the T-RFs near 0 s as a function of backazimuth. In this test, we assume full backazimuthal coverage of earthquakes and slowness range from 0.04 to 0.08 s km−1. Fig. 3 shows six representative layered velocity models and corresponding synthetic T-RFs in backazimuth within the data window. Anisotropy strength, its fast axis direction, and a thickness of the anisotropic layer are chosen arbitrarily in this test. We first examine the effect of a different thickness (0.5 or 3 km) of the topmost anisotropic layer with (1) a horizontal fast axis over an isotropic half space (Figs 3a and d) or with (2) a tilted (45°) fast axis over an isotropic half space (Figs 3g and j). Fig. 3(m) displays a model with two isotropic layers divided by an interface with a 10° dip, whereas Fig. 3(p) includes an additional anisotropic dipping layer on top of the isotropic medium. Also, we test the effect of a sensor misorientation of 1° using the six models. Figure 3. View largeDownload slide Six velocity structure models and corresponding synthetic T-RFs with backazimuth. The velocity model is shown in (a), (d), (g), (j), (m) and (p), and corresponding synthetic T-RFs in (b), (e), (h), (k), (n) and (q), assuming the sensor is properly oriented. The T-RFs from the sensor misoriented by 1° (clockwise rotation) are shown in (c), (f), (i), (l), (o) and (r). The thickness of a green layer (in panels m and p; above the location of a dashed line) and isotropic velocities of P- and S-waves of green and blue layers are taken from those of the upper and mid crusts in Kim et al. (2011). The north is towards the right. The RFs from the models (a, d, g and j) are calculated by the code anirec (Levin & Park 1997) and the RFs from the models (m and p) by the code raysum (Frederiksen & Bostock 2000). The angle of incidence is 25°. We note that the amplitudes of the synthetics are not scaled. The maximum amplitudes within −1 to 1 s are (b) 0.4 per cent, (c) 1.7 per cent, (e) 2.9 per cent, (f) 3.6 per cent, (h) 0.6 per cent, (i) 2.0 per cent, (k) 5.9 per cent, (l) 6.3 per cent, (n) 3.1 per cent, (o) 4.3 per cent, (q) 3.1 per cent, and (r) 4.3 per cent. Figure 3. View largeDownload slide Six velocity structure models and corresponding synthetic T-RFs with backazimuth. The velocity model is shown in (a), (d), (g), (j), (m) and (p), and corresponding synthetic T-RFs in (b), (e), (h), (k), (n) and (q), assuming the sensor is properly oriented. The T-RFs from the sensor misoriented by 1° (clockwise rotation) are shown in (c), (f), (i), (l), (o) and (r). The thickness of a green layer (in panels m and p; above the location of a dashed line) and isotropic velocities of P- and S-waves of green and blue layers are taken from those of the upper and mid crusts in Kim et al. (2011). The north is towards the right. The RFs from the models (a, d, g and j) are calculated by the code anirec (Levin & Park 1997) and the RFs from the models (m and p) by the code raysum (Frederiksen & Bostock 2000). The angle of incidence is 25°. We note that the amplitudes of the synthetics are not scaled. The maximum amplitudes within −1 to 1 s are (b) 0.4 per cent, (c) 1.7 per cent, (e) 2.9 per cent, (f) 3.6 per cent, (h) 0.6 per cent, (i) 2.0 per cent, (k) 5.9 per cent, (l) 6.3 per cent, (n) 3.1 per cent, (o) 4.3 per cent, (q) 3.1 per cent, and (r) 4.3 per cent. Four-lobed patterns with backazimuth are distinctively shown in Figs 3(b) and (e), and two-lobed patterns in Figs 3(h), (k), (n) and (q). In particular, the T-RFs from the models with the thin topmost anisotropic layer (Figs 3a, d, g and j) show distinctive polarity reversal in the vicinity of 0 s, caused by both the coupling between the P, SV and SH at the interface between the thin anisotropic and isotropic layers (Levin & Park 1998). Synthetic T-RFs from the model with an anisotropic topmost layer near 0 s, as expected, display waveform shape similar to the derivative of a Gaussian pulse (Figs 3b, e, h and k; Levin & Park 1998). Even if the sensor is misoriented by only 1°, the mean of the T-RFs from the six models is no longer zero (Figs 3c, f, i, l, o and r). Maximum amplitudes of the T-RFs within the data window in Figs 3(b), (e), (h), (k), (n) and (q) are 0.4, 2.9, 0.6, 5.9, 3.1 and 3.1 per cent (with respect to vertical P), respectively. If the sensor is misoriented by 1°, maximum amplitudes of the T-RFs in Figs 3(c), (f), (i), (l), (o) and (r) become 1.7, 3.6, 2.0, 6.3, 4.3 and 4.3 per cent (with respect to vertical P), respectively. For comparison, robust signals from the transition zone seismic discontinuities are typically observed at 2–5 per cent with respect to vertical P. This indicates that the sensitivity of T-RFs to the sensor orientation is probably on the order of 1° in this instance. 3.2 Synthetic test with non-uniform backazimuth data coverage We next test a case where the backazimuthal distribution of earthquakes is not uniform and incomplete. In this case, we set up a model with a 3 km thick topmost anisotropic layer (Fig. 3j), which produces the two-lobed pattern of T-RFs against backazimuth (Fig. 3k). Fig. 4 shows R- and T-RFs, their means, and their five harmonic terms in the case of complete and poor backazimuthal coverage. In the case of complete backazimuthal data coverage, we observe that the constant harmonic term in the T-RFs (HT1) is identical to a mean value of the T-RFs and they are both zero (Fig. 4a). In the case of incomplete backazimuthal distribution, the mean of the T-RFs is no longer zero near 0 s within the data window and can be biased in the direction where earthquakes are concentrated (Fig. 4b). We note that T-RFs with a single large data gap in backazimuth yields the largest peak in the mean of T-RFs than those with several small data gaps (e.g. Fig. 4b). However, HT1 remains zero regardless of incomplete backazimuthal distribution. See the Supporting Information (Fig. S1) for the case from the model (Fig. 3d), which shows the four-lobed pattern of T-RFs in backazimuth. Figure 4. View largeDownload slide Synthetic test results showing means and harmonic terms of both R- and T-RFs in the case of full (a) and poor (b) backazimuthal coverage of earthquakes. We use the model with the 3 km thick topmost anisotropic layer (Fig. 3j) to compute the R- and T-RFs (left panel), and their means and harmonics (right panel). The RF is sampled with an interval of 5° in backazimuth. Figure 4. View largeDownload slide Synthetic test results showing means and harmonic terms of both R- and T-RFs in the case of full (a) and poor (b) backazimuthal coverage of earthquakes. We use the model with the 3 km thick topmost anisotropic layer (Fig. 3j) to compute the R- and T-RFs (left panel), and their means and harmonics (right panel). The RF is sampled with an interval of 5° in backazimuth. 3.3 Synthetic test with the backazimuth coverage based on real earthquake distribution Using more realistic distribution of the earthquakes, we examine the effect of the variation in the slowness on the harmonic terms of the R- and T-RFs under the noise-free condition (Fig. 5). First, we consider the event distribution from the Korean seismic network, which well exceeds 90 per cent coverage of backazimuth during the period of ∼10 yr (Fig. 2). We compute synthetic R- and T-RFs, which are based on the backazimuth and slowness of the available 2430 earthquakes recorded from the station JJB (with the longest operation period) in South Korea from both P and PP waves. Then, these RFs are stacked with a bin of 5° in backazimuth. Fig. 5(a) shows that the amplitudes of the HT1 are nearly zero. Alternatively, we consider the earthquake distribution for both P and PP arrivals from the OBS J61C from the Cascadia Initiative (Toomey et al.2014), and calculate the harmonics of the RFs. We note that the backazimuth coverage for the OBS during 9 months is far less complete than that for Korean seismic network (Supporting Information Fig. S2). We still observe that the amplitudes of the HT1 are nearly zero from the OBS data (Fig. 5b). Fig. 6 shows the misfit functions f (φ) of the synthetic RFs using the model (Fig. 3j) and the earthquake distributions (e.g. Fig. 2 for the station JJB and Supporting Information Fig. S2 for the OBS J61C). Two minima of the f (φ) based on the synthetics occur at 0° even with the realistic variations in both backazimuth and slowness (Fig. 6). The ambiguity between the minima can be resolved by selecting the φ that makes the HR1 positive. Figure 5. View largeDownload slide Synthetic test results showing means and harmonic terms of both R- and T-RFs in the cases of realistic distributions of the earthquakes from (a) station JJB in South Korea (Fig. 2) and (b) OBS J61C (Cascadia Initiative). We use the model (Fig. 3j) to compute the R- and T-RFs (left panel), and their means and harmonics (right panel). The slowness is calculated with a 1-D velocity model (Kennett et al.1995). The RFs are stacked by a bin of 5° in backazimuth. Figure 5. View largeDownload slide Synthetic test results showing means and harmonic terms of both R- and T-RFs in the cases of realistic distributions of the earthquakes from (a) station JJB in South Korea (Fig. 2) and (b) OBS J61C (Cascadia Initiative). We use the model (Fig. 3j) to compute the R- and T-RFs (left panel), and their means and harmonics (right panel). The slowness is calculated with a 1-D velocity model (Kennett et al.1995). The RFs are stacked by a bin of 5° in backazimuth. Figure 6. View largeDownload slide The function f (φ) based on synthetic R- and T-RFs. The φmin with a positive sign of the mean of the HR1 is 0.0°. The variations in the backazimuth and slowness are from (a) South Korea (station JJB) and (b) Cascadia Initiative (OBS J61C). Figure 6. View largeDownload slide The function f (φ) based on synthetic R- and T-RFs. The φmin with a positive sign of the mean of the HR1 is 0.0°. The variations in the backazimuth and slowness are from (a) South Korea (station JJB) and (b) Cascadia Initiative (OBS J61C). 3.4 Case for various noise levels on synthetic data In this section, we further explore how various levels of noise and backazimuthal data coverage affect the orientation estimates. We generate synthetic seismograms using the model shown in Fig. 3(j) (Fig. 7a), considering various levels of random noise (0–200 per cent) and data gap in backazimuth (0°–300°). The 50 per cent level of noise means that the root-mean-square of the noise is close to a half of the amplitude of the vertical component. The same level of uncorrelated noise is added to all three components. A line in Fig. 7(b) roughly indicates the angle estimate deviated from the true orientation as the noise level and backazimuth coverage range vary. Two green boxes in Fig. 7(b) roughly indicate the ranges of both noise level and backazimuth of earthquakes from the Korean seismic network and the Cascadia Initiative. The upper limit of the noise level is a reciprocal of an average of the SNR in the vertical component. Figure 7. View largeDownload slide Synthetic test results showing how various noise levels and a degree of backazimuthal coverage affect the orientation estimate. (a) An example of three-component synthetic seismograms with random noise of 2.3 per cent at backazimuth of 110° where the tangential component is maximized, based on the model with the 3 km thick topmost anisotropic layer (Fig. 3j). Amplitudes of the seismograms are scaled to the maximum amplitude of the vertical-component seismogram. (b) The angle estimate (deviated from the true orientation) by varying levels of noise and backazimuth coverage range. Red dashed lines roughly indicate the uncertainty range of 0.2°–30°. Empty range in backazimuth is centered on 270° following the set up in Fig. 4(b). Green boxes are approximate ranges in cases of South Korea (SK) and the Cascadia Initiative (CI). Figure 7. View largeDownload slide Synthetic test results showing how various noise levels and a degree of backazimuthal coverage affect the orientation estimate. (a) An example of three-component synthetic seismograms with random noise of 2.3 per cent at backazimuth of 110° where the tangential component is maximized, based on the model with the 3 km thick topmost anisotropic layer (Fig. 3j). Amplitudes of the seismograms are scaled to the maximum amplitude of the vertical-component seismogram. (b) The angle estimate (deviated from the true orientation) by varying levels of noise and backazimuth coverage range. Red dashed lines roughly indicate the uncertainty range of 0.2°–30°. Empty range in backazimuth is centered on 270° following the set up in Fig. 4(b). Green boxes are approximate ranges in cases of South Korea (SK) and the Cascadia Initiative (CI). This synthetic test for the orientation estimate as a function of the noise level and the backazimuthal coverage provides useful information on determining a sensor operational period (Fig. 7b). As the operational period extends, the backazimuth data gap (horizontal axis of the plot) decreases and so is the uncertainty of the estimate, and vice versa (Fig. 7b). While the noise level is site-specific and can be decreased with filtering and applying SNR criterion, a high SNR condition can cause the poor backazimuthal coverage. 3.5 Correlation between radial and tangential harmonic terms Park & Levin (2016) analytically showed that HR2 and HR4 correlate negatively with HT3 and HT5, respectively, and HR3 and HT5 correlate positively with HT2 and HR4, respectively. The polarities of the correlations from our data are opposite to what Park & Levin (2016) presented because we define our tangential component opposite to what they used in the left-handed coordinate system in defining radial, tangential, and vertical direction. The correlations between the radial and tangential harmonic terms are observed except for the HR3, HT2, HR5 and HT4, which have nearly zero amplitudes (Fig. 4). The correlation between HR2 and HT3 is observed even in the case with the variations in slowness and backazimuth (Fig. 5). To quantify a degree of the observed correlation, we estimate a confidence range for non-randomness in the correlation between the harmonics in a time window of −2.0 to 2.0 s (Bendat & Piersol 2000). A degree of freedom is set as the length of the time window (4.0 s) times a corner frequency (2.5 and 1.5 Hz for South Korea and the Cascadia Initiative, respectively). 4 DATA ANALYSIS 4.1 Data acquisition and processing of Korean network data For the RF calculation, we collect teleseismic earthquakes of a magnitude greater than 5.5 recorded at 52 stations (which includes 30 borehole sensors) in South Korea (Fig. 1) from 2005 to 2016 (Fig. 2a). We analyse teleseismic P and PP arrivals to help mitigate uneven distribution of earthquake sources and provide a more filled backazimuthal coverage (Fig. 2a). For the P-wave RFs, the epicentral distance range of 30°–100° is chosen to avoid complex triplicated mantle P waves (less than 30°) and complication from the core-mantle boundary (distances greater than 100°). Similarly, for the PP-wave RFs, the epicentral distance range is chosen at 100°–180°. By including the PP-wave RFs, the backazimuthal gaps of P-RFs in the East Pacific Rise and the Mexico-Peru-Chile subduction zone (Fig. 2a) can be filled. The total number of earthquakes for P and PP phases is 3051 and 1651, respectively. All waveforms (neglecting the SNR) are cut to 30 s before and 180 s after P and PP arrival times before the RF calculation. The R- and T-RFs are calculated in the frequency domain with the water level of 10−2 (Langston 1979). In order to remove the high frequency noise, a Gaussian pulse with a half-width (1σ) of 2.5 Hz is convoluted with the RFs. We observe that the half-width from 1.0 to 4.0 Hz yields stable orientation estimates (Supporting Information Fig. S3). After the RF calculation, we stack both R- and T-RFs over a 5° interval in backazimuth. 4.2 Data acquisition and processing of OBS data We select earthquakes recorded from two broad-band OBSs (J61C and J39C) in the Cascadia Initiative. The OBSs J61C and J39C are deployed in deep water at 2673 m and 2656 m water depth and involve relatively lower noise level compared to those in shallow water (Lodewyk & Sumy 2015). From these two OBSs, the earthquakes with the magnitude greater than 5.5 are collected from September 2013 to June 2014. The ranges of the epicentral distance for the P and PP phases for the RF calculation are the same as those for the Korean seismic network. The numbers of earthquakes for the P and PP phases are 276 and 95, respectively (Supporting Information Fig. S2). All waveforms are cut to 50 s before and 250 s after P and PP arrival times before deconvolution. Prior to the deconvolution procedure, we apply a bandpass filter between 0.1 and 1.5 Hz, the same frequency band used by Audet (2016) for the Cascadia Initiative data and Akuhara & Mochizuki (2015) for OBSs offshore Japan. The lower cut-off frequency of 0.1 Hz is effective to remove an influence of infragravity wave. The water level in the deconvolution is set at 10−2. We stack both R- and T-RFs over a 5° interval in backazimuth. 5 RESULTS 5.1 Orientation estimates for Korean seismic network The orientations of all 52 stations in South Korea are presented in Table 1. The orientation for each station is estimated for a specific period (denoted as validation period), which is defined by a start and an end of the period when the polarity and/or the amplitude of the T-RFs in the data window shifts significantly. Alternatively, the validation period is manually set when the sensor replacement/maintenance record is available (Supporting Information Table S1). Table 1 includes a total of 65 measurements of orientation for different validation periods from two methods, which are (1) the minimization of HT1 in the data window and (2) the minimization of the amplitudes of a mean of the T-RFs. In Table 1 and Supporting Information Table S2, we include results from the method (2) for comparison. In addition, Supporting Information Table S2 shows orientation estimates based on non-stacking approach. The difference between the orientation estimates from the two methods is illustrated in Supporting Information Fig. S4, and it increases as the backazimuthal coverage decreases. However, all values are confined within and less than 1°, except for three stations that show large differences (2.3° of SEHB, 7.6° of INCN00 and 4.6° of KOHB; not shown in Supporting Information Fig. S4). We note that differences in orientations between two cases (with and without stacking) are overall small for Korean network data except for the stations INCN00, KOHB and NPR (Supporting Information Table S2). Table 1. Orientation correction estimates for the stations in South Korea. STA  Orientation  Error  %Backazimuth  Orientation  Validation period  Previously reported data (°)    (°; from harmonic decomposition)  (°; 1σ)  coverage  (°; from mean of T-RF)  From  To    BAR  5.5  0.4  94  5.5  2007–12  2015–01        BGD  359.5  0.6  100  359.6  2005–01  2015–01  359.8a  357.2b    BOSB  31.5  0.6  87  31.3  2012–08  2015–01  33.6c      BUS2  9.3  0.3  99  9.4  2005–01  2015–01        CHC2  3.8  0.6  90  3.9  2012–08  2015–01        CHJ2  10.5  0.4  99  10.6  2005–01  2015–01        CHNB  19.4  0.6  100  19.3  2005–01  2015–01  13.5a  15.3b    DACB  268.4  0.3  94  268.5  2010–08  2015–01  267.6c      DAG2  3.5  0.3  93  3.5  2010–08  2015–01        DGY2  9.6  0.4  94  9.5  2010–08  2015–01        EMSB  338.2  0.6  93  338.2  2012–12  2015–01  337.0c      EURB  104.4  0.9  86  105.1  2012–06  2014–08  106.0c      GAHB  128.5  0.3  94  128.6  2008–10  2015–01  130.7c      GKP1  5.4  0.5  97  5.4  2005–01  2015–01  4.2a  3.4b    GOCB  234.5  0.5  89  234.3  2012–08  2015–01  235.0c      GSU  146.2  0.6  99  146.5  2005–01  2015–01  142.6a  147.7c  140.8b  GWYB  120.5  0.7  90  120  2012–08  2015–01  123.3c      HALB  355.4  1.1  92  355.1  2012–08  2015–01  221.3c      HAMB  121.0  0.5  90  120.8  2012–08  2015–01  125.0c      HAWB  246.7  0.6  89  247.1  2012–08  2015–01  246.2c      HKU  8.2  0.5  96  8.3  2005–01  2015–01  9.1a  5.6b    HSB  144.7  0.8  83  144.5  2005–01  2007–01        HSB  228.5  0.4  96  228.6  2008–01  2015–01  224.3a  226.2c  221.6b  HWCB  3.3  0.3  93  3.4  2010–08  2015–01  2.4c      HWSB  60.9  0.9  96  60.7  2008–01  2012–10  59.0a  66.9c  59.0b  HWSB  –  –  –  –  2014–01  2014–06c  155.6c      IMWB  196.7  0.4  92  197.1  2012–08  2015–01  196.1c      INCN00  335.3  5.0  92  327.7  2009–10  2012–10  356d  359d    INCN00  4.6  1.1  90  4  2013–10  2016–01  358d      INCN10  359.2  0.6  94  359.2  2009–10  2012–10  355d  357d  359d  INCN10  0.1  0.4  90  359.7  2013–10  2016–01  0d      JEO2  2.7  0.4  89  3  2012–08  2015–01        JJB  237.1  0.7  94  236.5  2007–10  2012–10  235.8a  237.6b    JJB  233.8  1.2  82  233.3  2014–01  2016–01        JJU  1.9  0.6  100  1.8  2005–01  2015–01        JRB  194.7  0.3  99  194.7  2006–10  2015–01  193.7a  191.6c  192.5b  JSB  245.0  0.3  99  245.2  2006–10  2015–01  245.5a  245.3c  243.0b  KOHB  256.3  4.7  93  251.7  2009–08  2013–12  190.4c      KOHB  –  –  –  –  2014–04  2014–06c  3.5c      KSA  4.4  0.3  99  4.2  2005–01  2015–01  0.7a  1.2b    MGB  134.7  0.5  96  134.8  2005–03  2009–08  134.2a  133.2b    MGB  10.8  0.9  87  10.8  2010–11  2012–12        MGB  359.2  0.7  73  359.2  2013–06c  2014–06c  356.0c      NAWB  148.6  0.7  90  148.8  2012–08  2015–01  151.2c      NPR  2.9  1.0  100  3.2  2005–01  2015–01  4.3a  4.2b    OKCB  251.1  0.5  92  250.9  2012–08  2015–01  253.5c      OKEB  35.2  1.1  96  34.8  2011–10  2015–01  18.2c      SEHB  264.5  0.5  92  264.2  2008–11  2011–01        SEHB  3.5  1.7  59  1.2  2011–02  2012–07        SEHB  77.2  0.4  93  77.2  2012–08  2015–01  74.9c      SEO  358.1  0.3  100  358.1  2005–01  2015–01        SEO2  0.9  0.4  93  1.2  2010–08  2015–01        SES2  7.3  0.6  89  7.2  2012–08  2015–01        SHHB  326.3  0.4  92  326.2  2011–10  2015–01  325.6c      SMKB  355.0  3.8  63  355.9  2013–09c  2014–03c  2.4c      SMKB  172.6  0.9  44  172.3  2014–04c  2014–06c  165.4c      SND  0.9  0.2  99  0.9  2005–01  2015–01  359.0a  358.4b    SNU  3.8  0.6  92  3.8  2005–01  2007–09  5.5a,†  5.2b,†    SNU  –  –  –  –  2007–11  2008–01        SNU  3.5  0.3  96  3.9  2008–03  2015–01        TJN  1.6  0.4  99  1.5  2005–01  2015–01  0.1a  355.0b    ULJ2  358.0  0.3  94  358.1  2010–08  2015–01        ULL  9.7  0.9  96  9.4  2005–01  2009–12        ULLB  10.4  0.5  99  10.2  2006–10  2015–01  9.4a  20.0c  2.5b  YKB  194.9  0.3  97  194.9  2006–10  2015–01  191.0a  187.0c  195.7b  YNCB  309.0  0.7  90  309.1  2008–10  2010–08        YNCB  21.0  0.5  94  20.8  2011–01  2015–01  28.4c      YSB  13.6  0.4  94  13.8  2009–02  2015–01  13.8c      STA  Orientation  Error  %Backazimuth  Orientation  Validation period  Previously reported data (°)    (°; from harmonic decomposition)  (°; 1σ)  coverage  (°; from mean of T-RF)  From  To    BAR  5.5  0.4  94  5.5  2007–12  2015–01        BGD  359.5  0.6  100  359.6  2005–01  2015–01  359.8a  357.2b    BOSB  31.5  0.6  87  31.3  2012–08  2015–01  33.6c      BUS2  9.3  0.3  99  9.4  2005–01  2015–01        CHC2  3.8  0.6  90  3.9  2012–08  2015–01        CHJ2  10.5  0.4  99  10.6  2005–01  2015–01        CHNB  19.4  0.6  100  19.3  2005–01  2015–01  13.5a  15.3b    DACB  268.4  0.3  94  268.5  2010–08  2015–01  267.6c      DAG2  3.5  0.3  93  3.5  2010–08  2015–01        DGY2  9.6  0.4  94  9.5  2010–08  2015–01        EMSB  338.2  0.6  93  338.2  2012–12  2015–01  337.0c      EURB  104.4  0.9  86  105.1  2012–06  2014–08  106.0c      GAHB  128.5  0.3  94  128.6  2008–10  2015–01  130.7c      GKP1  5.4  0.5  97  5.4  2005–01  2015–01  4.2a  3.4b    GOCB  234.5  0.5  89  234.3  2012–08  2015–01  235.0c      GSU  146.2  0.6  99  146.5  2005–01  2015–01  142.6a  147.7c  140.8b  GWYB  120.5  0.7  90  120  2012–08  2015–01  123.3c      HALB  355.4  1.1  92  355.1  2012–08  2015–01  221.3c      HAMB  121.0  0.5  90  120.8  2012–08  2015–01  125.0c      HAWB  246.7  0.6  89  247.1  2012–08  2015–01  246.2c      HKU  8.2  0.5  96  8.3  2005–01  2015–01  9.1a  5.6b    HSB  144.7  0.8  83  144.5  2005–01  2007–01        HSB  228.5  0.4  96  228.6  2008–01  2015–01  224.3a  226.2c  221.6b  HWCB  3.3  0.3  93  3.4  2010–08  2015–01  2.4c      HWSB  60.9  0.9  96  60.7  2008–01  2012–10  59.0a  66.9c  59.0b  HWSB  –  –  –  –  2014–01  2014–06c  155.6c      IMWB  196.7  0.4  92  197.1  2012–08  2015–01  196.1c      INCN00  335.3  5.0  92  327.7  2009–10  2012–10  356d  359d    INCN00  4.6  1.1  90  4  2013–10  2016–01  358d      INCN10  359.2  0.6  94  359.2  2009–10  2012–10  355d  357d  359d  INCN10  0.1  0.4  90  359.7  2013–10  2016–01  0d      JEO2  2.7  0.4  89  3  2012–08  2015–01        JJB  237.1  0.7  94  236.5  2007–10  2012–10  235.8a  237.6b    JJB  233.8  1.2  82  233.3  2014–01  2016–01        JJU  1.9  0.6  100  1.8  2005–01  2015–01        JRB  194.7  0.3  99  194.7  2006–10  2015–01  193.7a  191.6c  192.5b  JSB  245.0  0.3  99  245.2  2006–10  2015–01  245.5a  245.3c  243.0b  KOHB  256.3  4.7  93  251.7  2009–08  2013–12  190.4c      KOHB  –  –  –  –  2014–04  2014–06c  3.5c      KSA  4.4  0.3  99  4.2  2005–01  2015–01  0.7a  1.2b    MGB  134.7  0.5  96  134.8  2005–03  2009–08  134.2a  133.2b    MGB  10.8  0.9  87  10.8  2010–11  2012–12        MGB  359.2  0.7  73  359.2  2013–06c  2014–06c  356.0c      NAWB  148.6  0.7  90  148.8  2012–08  2015–01  151.2c      NPR  2.9  1.0  100  3.2  2005–01  2015–01  4.3a  4.2b    OKCB  251.1  0.5  92  250.9  2012–08  2015–01  253.5c      OKEB  35.2  1.1  96  34.8  2011–10  2015–01  18.2c      SEHB  264.5  0.5  92  264.2  2008–11  2011–01        SEHB  3.5  1.7  59  1.2  2011–02  2012–07        SEHB  77.2  0.4  93  77.2  2012–08  2015–01  74.9c      SEO  358.1  0.3  100  358.1  2005–01  2015–01        SEO2  0.9  0.4  93  1.2  2010–08  2015–01        SES2  7.3  0.6  89  7.2  2012–08  2015–01        SHHB  326.3  0.4  92  326.2  2011–10  2015–01  325.6c      SMKB  355.0  3.8  63  355.9  2013–09c  2014–03c  2.4c      SMKB  172.6  0.9  44  172.3  2014–04c  2014–06c  165.4c      SND  0.9  0.2  99  0.9  2005–01  2015–01  359.0a  358.4b    SNU  3.8  0.6  92  3.8  2005–01  2007–09  5.5a,†  5.2b,†    SNU  –  –  –  –  2007–11  2008–01        SNU  3.5  0.3  96  3.9  2008–03  2015–01        TJN  1.6  0.4  99  1.5  2005–01  2015–01  0.1a  355.0b    ULJ2  358.0  0.3  94  358.1  2010–08  2015–01        ULL  9.7  0.9  96  9.4  2005–01  2009–12        ULLB  10.4  0.5  99  10.2  2006–10  2015–01  9.4a  20.0c  2.5b  YKB  194.9  0.3  97  194.9  2006–10  2015–01  191.0a  187.0c  195.7b  YNCB  309.0  0.7  90  309.1  2008–10  2010–08        YNCB  21.0  0.5  94  20.8  2011–01  2015–01  28.4c      YSB  13.6  0.4  94  13.8  2009–02  2015–01  13.8c      Notes: The orientation angle ranges from 0° to 360° with clockwise direction, and our values indicate the angles that we need to correct from the misaligned north. We note that we apply this correction scheme to previous estimates (Ekström & Busby 2008; Shin et al.2009; Lee & Rhie 2015; Lee & Sheen 2015). ‘00’ and ‘10’: location code given for the station INCN (http://ds.iris.edu/mda/IU/INCN) by the Data Management Center (DMC) of the Incorporated Research Institutions for Seismology (IRIS). aShin et al. (2009). bLee & Rhie (2015). cLee & Sheen (2015). dEkström & Busby (2008); http://www.ldeo.columbia.edu/∼ekstrom/Projects/WQC/COMB_QC/POL_IU_S_ALL_EPO.html. †Their periods of used data (January 2007 to September 2008) contain the replacement dates of sensors: 23 October 2007 and 18 February 2008 (Supporting Information Table S1). View Large Uncertainty estimates (1σ) of the orientation from the method (1) are shown in Table 1. The 57 out of 65 measurements have an uncertainty of less than 1° (Table 1). Table 1 also includes previous estimates for some stations (Ekström & Busby 2008; Shin et al.2009; Lee & Rhie 2015; Lee & Sheen 2015) for comparisons. We note that estimates from station SMKB (for a validation period of April 2014–June 2014) involve relatively large uncertainty value (Table 1) due to poor backazimuthal coverage (44 per cent). Lee & Sheen (2015) previously reported changes in the orientation for that station, which have not been reported by the seismic network operator. We confirm that significant temporal variations of the polarity and amplitude of the T-RFs likely had stemmed during sensor replacement/maintenance (e.g. Supporting Information Table S1). Of the 52 stations examined, we present results from three stations (HSB, MGB and SEHB) that show substantial changes in both polarity and amplitude of the T-RFs. First, we illustrate such changes by showing T-RFs stacked over an interval of two months without overlap (Supporting Information Figs S5a, c and e for stations HSB, MGB and SEHB, respectively). Second, we show the maximum and minimum amplitudes of the stacked T-RFs in the data window (Supporting Information Figs S5b, d and f for stations HSB, MGB and SEHB, respectively; Supporting Information Fig. S6 for the other stations). Our results show that six borehole stations experience sudden orientation change ranging from 12° (at station MGB) to 178° (at station SMKB) over the period of 10 yr (Table 1). Fig. 8 shows the misfit function f (φ) of stations HSB, MGB and SEHB at different periods. The data from station HSB shows that f (φ) for the period of January 2005–January 2007 is minimized at 144.7° with a plus sign for the mean of HR1 (Fig. 8a). The orientation for the period of January 2008– January 2015 is estimated as 228.5° (48.5° + 180°) because f (φ) is minimized at 48.5° with a minus sign for the mean of HR1 (Fig. 8b). The orientation estimates for station MGB are 134.7°, 10.8° and 359.2° (179.2° + 180°) for three different periods (Figs 8c–e). The estimates for station SEHB are 264.5° (84.5° + 180°) and 77.2° in two different periods (Figs 8f and g). Figure 8. View largeDownload slide The function f (φ) for stations HSB, MGB and SEHB plotted against φ for different periods. The φmin for the different periods are 144.7° and 48.5° for HSB (a and b), 134.7°, 10.8° and 179.2° for MGB (c, d, e), and 84.5° and 77.2° for SEHB (f and g). Figure 8. View largeDownload slide The function f (φ) for stations HSB, MGB and SEHB plotted against φ for different periods. The φmin for the different periods are 144.7° and 48.5° for HSB (a and b), 134.7°, 10.8° and 179.2° for MGB (c, d, e), and 84.5° and 77.2° for SEHB (f and g). We show temporal changes of f (φ) and φmin for stations HSB, MGB and SEHB with a moving window of one year (Fig. 9, background colour and solid black lines, respectively). The orientation angles for the stations and their validation periods (Table 1) are also marked as references (Fig. 9; white ticked lines). We find periods that involve a sudden change in φmin, which are 2007–2008 for HSB (Fig. 9a) and October 2009–April 2010, June 2010–December 2010 and October 2012–February 2013 for MGB (Fig. 9b). These correspond to the time when the orientation had changed, which is consistent with a record of the sensor replacement/maintenance (Supporting Information Table S1). Besides the record, we also identify other periods when either maximum or minimum amplitude of the T-RFs is statically changed (Supporting Information Figs S5b, d and f for stations HSB, MGB and SEHB, respectively). We report that the period of 2011–2012 for station SEHB involves a significant change in φmin (Fig. 9c), which is considered to be the orientation change but not noted by the seismic network operator. In addition, we observe that both f (φ) and φmin slightly deviate from the reference during the periods of August 2005–December 2006, October 2008–August 2009 and June 2013–October 2014 for HSB (Fig. 9a) and during the periods of February 2007–June 2007 and October 2011– October 2012 for MGB (Fig. 9b). Such periods moderately coincide with the temporal fluctuation of backazimuthal coverage of earthquakes (Supporting Information Figs S5a and b). Figure 9. View large Download slide The f (φ) and φmin as a function of time for stations HSB, MGB and SEHB. The f (φ) is calculated with an interval of two months and a moving window of 1 yr. The f (φ) and φmin are represented by colour and a solid black line, respectively. The validation periods for our orientation estimates (Table 1) are indicated by white ticked lines. Blue arrows indicate reported replacement date of sensors (Supporting Information Table S1). Grey background colour means no data. Figure 9. View large Download slide The f (φ) and φmin as a function of time for stations HSB, MGB and SEHB. The f (φ) is calculated with an interval of two months and a moving window of 1 yr. The f (φ) and φmin are represented by colour and a solid black line, respectively. The validation periods for our orientation estimates (Table 1) are indicated by white ticked lines. Blue arrows indicate reported replacement date of sensors (Supporting Information Table S1). Grey background colour means no data. We next show how the harmonic terms of R- and T-RFs, and also means both R- and T-RFs, for the three stations change with the orientation correction based on the method (1) (Fig. 10). The amplitudes of HT1 for the different periods are minimized and essentially nearly zero after the orientation correction (Fig. 10). We now observe strong coherence in waveform between the harmonic terms for different validation periods, and between the means after the correction. This correlation is expected unless there is any temporal change in structure and significant change in the distribution of the earthquakes. We observe that HR2 and HR4 correlate negatively with HT3 and HT5, respectively, and HR3 and HT5 correlate positively with HT2 and HR4, respectively (Fig. 10; Supporting Information Figs S7–S9). We note that the amplitudes of HR4, HR5, HT4, and HT5 are small. Thus, we do not see a clear correlation in waveform between the different harmonics (see Supporting Information Figs S7–S9). Figure 10. View large Download slide Harmonic terms of R- and T-RFs, and their means for the stations HSB, MGB and SEHB prior to (left panel) and after the correction (right panel). Note that the amplitudes of the HT1 are minimized to nearly zero by the orientation correction. See the Supporting Information (Figs S7–S9) for HR4, HR5, HT4 and HT5 for the stations. Figure 10. View large Download slide Harmonic terms of R- and T-RFs, and their means for the stations HSB, MGB and SEHB prior to (left panel) and after the correction (right panel). Note that the amplitudes of the HT1 are minimized to nearly zero by the orientation correction. See the Supporting Information (Figs S7–S9) for HR4, HR5, HT4 and HT5 for the stations. Lastly, we show T-RFs and HT1 for the selected station HSB with and without the orientation correction (Fig. 11). The similar two-lobed patterns of the T-RFs with backazimuth are revealed by the correction (Figs 11c and d). See Supporting Information Fig. S10 for the stations MGB and SEHB, respectively. Figure 11. View large Download slide Stacked T-RFs plotted according to backazimuths for station HSB for each different validation period, prior to (a and b) and after the orientation correction (c and d). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). See Supporting Information Fig. S10 for stations MGB and SEHB. Figure 11. View large Download slide Stacked T-RFs plotted according to backazimuths for station HSB for each different validation period, prior to (a and b) and after the orientation correction (c and d). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). See Supporting Information Fig. S10 for stations MGB and SEHB. 5.2 Orientation estimates for OBS network Our estimates based on the harmonic decomposition are compared against previous estimates based on the minimization of stacked T-RFs at 0.0–5.0 s (Janiszewski & Abers 2015) and surface wave arrival angles (Lodewyk & Sumy 2015; Table 2). Our estimates are similar to the estimates by Lodewyk & Sumy (2015) and the estimates by Janiszewski & Abers (2015). We note that slight deviation between the estimates might stem from the different data selection criteria and processing. Table 2. Orientation correction estimates for the stations from Cascadia Initiative (Toomey et al.2014). STA  Orientation (°; from harmonic decomposition)  Error (°; 1σ)  %Backazimuth coverage  Orientation (°; from mean of T-RF)  Previously reported data (°; 1σ)  J61C  190.2  6.6  79  185.9  200 ± 2a  187 ± 9b  J39C  89.0  4.1  79  97.1  105 ± 1a  92 ± 9b  STA  Orientation (°; from harmonic decomposition)  Error (°; 1σ)  %Backazimuth coverage  Orientation (°; from mean of T-RF)  Previously reported data (°; 1σ)  J61C  190.2  6.6  79  185.9  200 ± 2a  187 ± 9b  J39C  89.0  4.1  79  97.1  105 ± 1a  92 ± 9b  Notes: The orientation angle ranges from 0° to 360° with clockwise direction, and our values indicate the angles that we need to correct from the misaligned north. We note that we apply this correction scheme to previous estimates (Janiszewski & Abers 2015; Lodewyk & Sumy 2015). For the uncertainty estimates, Janiszewski & Abers (2015) used the 95 per cent confidence bounds from the F-test, with degrees of freedom determined from the net filter response of the signal. Our uncertainty estimates for J61C and J39 from the F-test are 3.3° and 2.5°, respectively. The difference in the two estimates for each station may arise from slightly different dataset and filtering. aJaniszewski & Abers (2015). bLodewyk & Sumy (2015). View Large We observe very small positive amplitudes of HR1 near 0 s even after the orientation corrections (Fig. 12). Such small amplitude of HR1 makes the range of f (φ) very narrow (Supporting Information Fig. S11). From the data from J61C, the orientation is estimated as 190.2° (10.2° + 180°) (Supporting Information Fig. S11a). The orientation estimate for J39C is 89.0° (Supporting Information Fig. S11b). Figure 12. View largeDownload slide Harmonic terms of R- and T-RFs, and their means for the OBSs (a) J61C and (b) J39C prior to (first two columns) and after the correction (third and fourth columns). Figure 12. View largeDownload slide Harmonic terms of R- and T-RFs, and their means for the OBSs (a) J61C and (b) J39C prior to (first two columns) and after the correction (third and fourth columns). The correlations between HR2 and HT3 and between HR3 and HT2 of the J61C are observed strong, whereas those between HR4 and HT5 and between HR5 and HT4 are weak (Fig. 12a). The correlations between the harmonic terms for J39C are not all visibly clear because of small amplitudes at and near 0 s (Fig. 12b). Fig. 13 shows the T-RFs before and after the orientation correction for the two OBSs. Figure 13. View largeDownload slide Stacked T-RFs plotted according to backazimuths for the OBSs J61C and J39C, prior to (a and b) and after the orientation correction (c and d). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). Figure 13. View largeDownload slide Stacked T-RFs plotted according to backazimuths for the OBSs J61C and J39C, prior to (a and b) and after the orientation correction (c and d). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). 6 DISCUSSION Most indirect seismic methods in determining the orientation inevitably involve inaccuracy due to uncertain medium seismic properties beneath the receiver. The medium is assumed to be isotropic and homogeneous in the horizontal direction (e.g. Ekström & Busby 2008; Shin et al.2009; Zha et al.2013; Lee & Rhie 2015; Lee & Sheen 2015; Wang et al.2016). If the anisotropy of the medium were strong, the particle motion of the P wave is no longer parallel to its propagation direction. The same applies to the Rayleigh wave. Furthermore, the minimum-time path between the source and receiver might not coincide with the shortest-distance path if the medium is laterally heterogeneous. Although our estimates based on the teleseismic RFs can also be influenced by the uncertainty in the nature of the medium, the degree of its dependency can be much less than that of the studies based on the particle motions. Our method works well in the presence of anisotropy and/or dipping structures beneath the receiver (e.g. Fig. 3). Nevertheless, to make a single measurement of orientation with an uncertainty less than 1°, we require a relatively large dataset with nearly full backazimuthal coverage of teleseismic earthquakes (e.g. Fig. 2). 6.1 Orientation estimates of the Korean seismic network Most of the orientations for stations in South Korea are measured with over 85 per cent of backazimuthal coverage (Table 1). Using both P and PP phases, the 70, 80 and 90 per cent of the backazimuthal coverage are guaranteed by over 8, 10 and 21 months of station operational periods (Figs 2b and 14) or 250, 350 and 750 earthquakes, respectively (Supporting Information Fig. S12a). Even in the case of relatively poor backazimuthal coverage at station MGB (73 per cent), our value (359.2° ± 0.7°; Table 1) still agrees reasonably well with the previous estimate (Lee & Sheen 2015; a difference of 3.2°). It is difficult to statistically compare our estimate with the previous estimate because different types of data (e.g. teleseismic body-wave phase versus regional body-wave phase versus) and different data processing (e.g. passband) are involved in determining the orientation. The measurement for station SEHB includes the smallest backazimuthal coverage (59 per cent), and our estimate for this station is 3.5° ± 1.7°. Figure 14. View largeDownload slide Global distribution of the earthquakes and the period for the 70 per cent backazimuthal coverage. (a) Distribution of the earthquakes in 14 months from November 2007 with magnitude larger than 5.5. (b) The variation of the period for the 70 per cent of the backazimuthal coverage using P and PP (distance ranging from 30° to 180°). Figure 14. View largeDownload slide Global distribution of the earthquakes and the period for the 70 per cent backazimuthal coverage. (a) Distribution of the earthquakes in 14 months from November 2007 with magnitude larger than 5.5. (b) The variation of the period for the 70 per cent of the backazimuthal coverage using P and PP (distance ranging from 30° to 180°). For Korean network data (the case with the nearly full backazimuthal coverage), the backazimuthal coverage at 20 months reaches to about 90 per cent (Fig. 15a; 750 earthquakes, Supporting Information Fig. S12a), and the orientation (φmin) converges to a steady value after 20 months (Fig. 15c). The orientation angle at 20 months slightly deviates from the best estimate by 0.3°, which is smaller than the uncertainty. We note that just over three months the measured orientation deviates from the best angle only by 6° with about 60 per cent backazimuthal coverage. Figure 15. View largeDownload slide Convergence of the orientation (φmin) and corresponding backazimuthal coverage, and correlation (in −2.0 to 2.0 s) between the harmonic terms with cumulative data (a, b, c) for the station CHJ2 in South Korea and (d, e, f) for the OBS J61C of the Cascadia Initiative. The confidence level for non-randomness (Bendat & Piersol 2000) is marked as pink colour. Figure 15. View largeDownload slide Convergence of the orientation (φmin) and corresponding backazimuthal coverage, and correlation (in −2.0 to 2.0 s) between the harmonic terms with cumulative data (a, b, c) for the station CHJ2 in South Korea and (d, e, f) for the OBS J61C of the Cascadia Initiative. The confidence level for non-randomness (Bendat & Piersol 2000) is marked as pink colour. The harmonics HR2 and HT3 and the HR3 and HT2 correlate negatively and positively, respectively (Fig. 15b). The observed correlation between the harmonic terms is quantified in terms of a confidence range for randomness in the correlation (Bendat & Piersol 2000) (Fig. 15b, pink colour). While the correlation value between the HR2 and HT3 reaches to −1 just after one month, the correlation between the HR3 and HT2 slowly converges to 0.6 (Fig. 15b). The observed slow convergence and relatively low correlation (still within the confidence range of 95 per cent) for the HR3 and HT2 are primarily due to small amplitudes of the HR3 and HT2. Similarly, the correlation between the HR4 and HT5 and HR5 and HT4 are scattered within the confidence range during the whole period of 10 yr. Assuming backazimuthal coverage of earthquakes is good, strong correlation (beyond the confidence range of 95 per cent) between the orthogonal harmonic terms may further supports the reliability of our orientation estimates (Figs 10 and 15). The orientation angle is deviated by about 5° from the final estimate during the period of 2005−mid-2006, although a degree of correlations between HR3 and HT2 and between HR5 and HT4 in 2005−mid-2006 is observed high (Fig. 15b). We observe that the harmonic terms from most stations (except for the stations in oceanic islands) in South Korea share similar patterns in waveforms within the data window near 0 s (Fig. 10), and we suspect that this is in part due to the structure. For instance, a negative pulse just before 0 s and a positive pulse just after 0 s are shown in the HR2 (conversely in the HT3) of the stations (Fig. 10). The common features of the harmonic terms might indicate an approximately southward dipping of the fast axis in the anisotropic top layer beneath South Korea, as illustrated in Figs 3(g) and (j). A separate analysis of P coda phases is required for a precise investigation on the crustal seismic structure of Korea. 6.2 Applicability of the method on OBS data One clear difference in the harmonic terms of the RFs between the land-based stations and OBSs is the amplitude of the HR1 within the data window near and at 0 s. Observed small amplitudes of the HR1 (∼0.1 per cent) in the OBS data can be a source of a problem in the orientation correction for the OBS (Fig. 12). Small amplitudes of the R-RFs typically arise from a nearly vertical incidence of the teleseismic P phase due to the presence of the sedimentary layer with low P-wave velocity (Kawakatsu & Abe 2016). The pseudo vertical incidence of P can yield a very shallow range of the f (φ), which makes it difficult to determine φmin in the presence of noise since the orientation estimate depends on the amplitude of the HR1 near 0 s. Despite the high background noise embedded in the OBS data (Webb 1998; Lin et al.2010), the RFs at OBS station J61C show clear correlation between HR2 and HT3 and between HR3 and HT2 (Fig. 12a), as expected in an anisotropic or dipping media (Park & Levin 2016). Our estimates generally agree with previous estimates based on slightly different data processing scheme (Janiszewski & Abers 2015) and different seismic phase (Lodewyk & Sumy 2015). Because of shorter operation period for the OBSs (∼9 months), the backazimuthal coverage reaches to about 80 per cent in 9 months (Fig. 15d; 380 earthquakes in Supporting Information Fig. S12b) and the orientation estimate (φmin) just begins to converge (Fig. 15f). If the noise level is assumed at 80 per cent in the OBS data, about 70 per cent of backazimuthal coverage is required to obtain the orientation estimate with 5° uncertainty (Fig. 7b). For example, this can be achieved approximately within 4−7 months of the deployment in Pacific, Atlantic, Indian Oceans, and North America and 9−14 months in southeastern Asia, South America and Africa (Fig. 14). 7 CONCLUSIONS We calculate teleseismic receiver functions from the P and PP phases using the harmonic decomposition method and determine the orientation of a seismometer by minimizing a constant term in a harmonic expansion of T-RFs in backazimuth near and at 0 s. This approach is applied to data from 52 stations in South Korea and two OBSs from the Cascadia Initiative project. Our method works well in the presence of anisotropy and/or dipping structures beneath the receiver. Also, the method yields a much more robust estimate than those obtained by minimizing the average amplitude of the T-RFs near 0 s. In particular, the method effectively determines the orientation in the case of deficient backazimuthal coverage of earthquakes (i.e. from the OBSs). The station operational period of about one year guarantees ∼80 per cent of the backazimuthal coverage with the usage of both P and PP phases from earthquakes with a magnitude over 5.5 recorded in the seismic network in South Korea. Our analysis shows that the operational period of the station of 20 months (recording about 750 earthquakes) is required for a single measurement of orientation with an uncertainty less than 1°. We keep track of the history of the orientation of a total 52 broad-band velocity seismometers for the period of 2005–2016 by detecting its period when the polarity and the amplitude of the T-RFs are largely changed. Of the 30 borehole stations in South Korea, six stations are confirmed to experience a significant orientation change (12° to ∼178°) during the 10 yr period. Our method can be applied to any broad-band data, regardless of sensor types and orientation, for seismic studies (i.e. receiver functions). ACKNOWLEDGEMENTS Two sets of seismic data from Korea Institute of Geoscience and Mineral Resources (KIGAM) and Korea Meteorological Administration (KMA) are available with permissions at http://quake.kigam.re.kr (last accessed September 2014) and http://necis.kma.go.kr (last accessed April 2015), respectively. The data of INCN (Fig. 1) can be obtained from the IRIS Data Management Center at www.iris.edu (last accessed May 2015). Earthquake catalogue was obtained at US Geological Survey (http://earthquake.usgs.gov/earthquakes/search; last accessed March 2016). Receiver function calculation was performed using processRFmatlab (https://github.com/iwbailey/processRFmatlab; see Bailey et al. (2012) for detail). Most of the figures are plotted by Generic Mapping Tools version 4.5.9 (http://www.soest.hawaii.edu/gmt; Wessel & Smith 1998). YK and HL acknowledge supports from the Korea Meteorological Administration Research and Development Program under Grant KMIPA2015-7020, and from National Research Foundation of Korea Grant funded by the Korean Government (NRF-2014S1A2A2027609). TAS acknowledges the support from the Natural Environment Research Council (NE/P001378/1). Finally, authors thank the editor and two reviewers for helpful comments which improved the manuscript. REFERENCES Akuhara T., Mochizuki K., 2015. Hydrous state of the subducting Philippine Sea plate inferred from receiver function image using onshore and offshore data, J. geophys. Res. , 120( 12), 8461– 8477. Google Scholar CrossRef Search ADS   Audet P., 2015. Layered crustal anisotropy around the San Andreas fault near Parkfield, California, J. geophys. Res. , 120( 5), 3527– 3543. Google Scholar CrossRef Search ADS   Audet P., 2016. Receiver functions using OBS data: promises and limitations from numerical modelling and examples from the Cascadia Initiative, Geophys. J. Int. , 205( 3), 1740– 1755. https://doi.org/10.1093/gji/ggw111 Google Scholar CrossRef Search ADS   Bailey I.W., Miller M.S., Liu K., Levander A., 2012. Vs and density structure beneath the Colorado Plateau constrained by gravity anomalies and joint inversions of receiver function and phase velocity data, J. geophys. Res. , 117, B02313, doi:10.1029/2011JB008522. Bendat J.S., Piersol A.G., 2000. Random Data Analysis and Measurement Procedures , 3rd edn, pp. 108– 111, Wiley. Bianchi I., Park J., Agostinetti N.P., Levin V., 2010. Mapping seismic anisotropy using harmonic decomposition of receiver functions: an application to Northern Apennines, Italy, J. geophys. Res. , 115( B12), doi:10.1029/2009JB007061. https://doi.org/10.1029/2009JB007061 Cassidy J., 1992. Numerical experiments in broadband receiver function analysis, Bull. seism. Soc. Am. , 82, 1453– 1474. Ekström G., Busby R.W., 2008. Measurements of seismometer orientation at USArray transportable array and backbone stations, Seismol. Res. Lett. , 79( 4), 554– 561. https://doi.org/10.1785/gssrl.79.4.554 Google Scholar CrossRef Search ADS   Farra V., Vinnik L., 2000. Upper mantle stratification by P and S receiver functions, Geophys. J. Int. , 141( 3), 699– 712. https://doi.org/10.1046/j.1365-246x.2000.00118.x Google Scholar CrossRef Search ADS   Frederiksen A.W., Bostock M.G., 2000. Modelling teleseismic waves in dipping anisotropic structures, Geophys. J. Int. , 141( 2), 401– 412. https://doi.org/10.1046/j.1365-246x.2000.00090.x Google Scholar CrossRef Search ADS   Girardin N., Farra V., 1998. Azimuthal anisotropy in the upper mantle from observations of P-to-S converted phases: application to southeast Australia, Geophys. J. Int. , 133( 3), 615– 629. https://doi.org/10.1046/j.1365-246X.1998.00525.x Google Scholar CrossRef Search ADS   Janiszewski H.A., Abers G.A., 2015. Imaging the plate interface in the Cascadia seismogenic zone: new constraints from offshore receiver functions, Seismol. Res. Lett. , 86, doi:10.1785/0220150104. Kawakatsu H., Abe Y., 2016. Significance of sediment reverberations on receiver functions of broadband OBS data: Comments on Olugboji et al. (2016) nature of the seismic lithosphere-asthenosphere boundary within normal oceanic mantle from high-resolution receiver functions, Geochem. Geophys. Geosyst. , 17( 8), 3488– 3492. https://doi.org/10.1002/2016GC006418 Google Scholar CrossRef Search ADS   Kennett B.L.N., Engdahl E.R., Buland R., 1995. Constraints on seismic velocities in the Earth from traveltimes, Geophys. J. Int. , 122( 1), 108– 124. https://doi.org/10.1111/j.1365-246X.1995.tb03540.x Google Scholar CrossRef Search ADS   Kim S., Rhie J., Kim G., 2011. Forward waveform modelling procedure for 1-D crustal velocity structure and its application to the southern Korean Peninsula, Geophys. J. Int. , 185( 1), 453– 468. https://doi.org/10.1111/j.1365-246X.2011.04949.x Google Scholar CrossRef Search ADS   Langston C.A., 1979. Structure under Mount Rainier, Washington, Inferred from Teleseismic Body Waves, J. geophys. Res. , 84( B9), 4749–4762. https://doi.org/10.1029/JB084iB09p04749 Larson E.W., 2000. Measuring Refraction and Modeling Velocities of Surface Waves, PhD dissertation , Harvard University, Cambridge, Massachusetts. Larson E.W., Ekström G., 2002. Determining surface wave arrival angle anomalies, J. geophys. Res. , 107( B6), doi:10.1029/2000JB000048. https://doi.org/10.1029/2000JB000048 Laske G., 1995. Global observation of off-great-circle propagation of long-period surface waves, Geophys. J. Int. , 123( 1), 245– 259. https://doi.org/10.1111/j.1365-246X.1995.tb06673.x Google Scholar CrossRef Search ADS   Laske G., Masters G., 1996. Constraints on global phase velocity maps from long-period polarization data, J. geophys. Res. , 101( B7), 16 059–16 075. https://doi.org/10.1029/96JB00526 Google Scholar CrossRef Search ADS   Lee H., Sheen D., 2015. A study on determination of orientation of borehole seismometer, J. Geol. Soc. Korea , 51( 1), 93– 103 (in Korean). https://doi.org/10.14770/jgsk.2015.51.1.93 Google Scholar CrossRef Search ADS   Lee S.-J., Rhie J., 2015. Determining the orientations of broadband stations in South Korea using ambient noise cross-correlation, Geophys. Geophys. Explor. , 18( 2), 85– 90 (in Korean) https://doi.org/10.7582/GGE.2015.18.2.085 Google Scholar CrossRef Search ADS   Levin V., Park J., 1997. P-SH conversions in a flat-layered medium with anisotropy of arbitrary orientation, Geophys. J. Int. , 131( 2), 253– 266. https://doi.org/10.1111/j.1365-246X.1997.tb01220.x Google Scholar CrossRef Search ADS   Levin V., Park J., 1998. P-SH conversions in layered media with hexagonally symmetric anisotropy: a cookbook, in Geodynamics of Lithosphere & Earth's Mantle , pp. 669– 697, Springer, doi:10.1007/s000240050136. Lin C.R., Kuo B.Y., Liang W.T., Chi W.C., Huang Y.C., Collins J., Wang C.Y., 2010. Ambient noise and teleseismic signals recorded by ocean-bottom seismometers offshore eastern Taiwan, Terr. Atmos. Ocean. Sci. , 21( 5), 743– 755. https://doi.org/10.3319/TAO.2009.09.14.01(T) Google Scholar CrossRef Search ADS   Lodewyk J., Sumy D., 2015. Cascadia Amphibious Array Ocean Bottom Seismograph Horizontal Component Orientations, 2013–2014 OBS Deployments, OBSIP Management Office Incorporated Research Institutions for Seismology. Long M.D., Becker T.W., 2010. Mantle dynamics and seismic anisotropy, Earth planet. Sci. Lett. , 297( 3–4), 341– 354. https://doi.org/10.1016/j.epsl.2010.06.036 Google Scholar CrossRef Search ADS   Long M.D., Silver P.G., 2009. Shear wave splitting and mantle anisotropy: measurements, interpretations, and new directions, Surv. Geophys. , 30( 4–5), 407– 461. https://doi.org/10.1007/s10712-009-9075-1 Google Scholar CrossRef Search ADS   Park J., Levin V., 2016. Anisotropic shear zones revealed by backazimuthal harmonics of teleseismic receiver functions, Geophys. J. Int. , 207( 2), 1216– 1243. https://doi.org/10.1093/gji/ggw323 Google Scholar CrossRef Search ADS   Schulte-Pelkum V., Mahan K.H., 2014. A method for mapping crustal deformation and anisotropy with receiver functions and first results from USArray, Earth planet. Sci. Lett. , 402, 221– 233. https://doi.org/10.1016/j.epsl.2014.01.050 Google Scholar CrossRef Search ADS   Schulte-Pelkum V., Masters G., Shearer P.M., 2001. Upper mantle anisotropy from long-period P-polarization, J. geophys. Res. , 106( B10), 21 917–21 934. https://doi.org/10.1029/2001JB000346 Google Scholar CrossRef Search ADS   Shin J.S., Sheen D.-H., Shin I.C., 2009. Orientation correction for borehole seismic stations in South Korea, J. Geol. Soc. Korea , 45, 47– 54 (in Korean). Stachnik J.C., Sheehan A.F., Zietlow D.W., Yang Z., Collins J., Ferris A., 2012. Determination of New Zealand ocean bottom seismometer orientation via Rayleigh-wave polarization, Seismol. Res. Lett. , 83( 4), 704– 713. https://doi.org/10.1785/0220110128 Google Scholar CrossRef Search ADS   Toomey D.R. et al.  , 2014. The Cascadia Initiative: a sea change in seismological studies of subduction zones, Oceanography , 27( 2), 138– 150. https://doi.org/10.5670/oceanog.2014.49 Google Scholar CrossRef Search ADS   Vinnik L., 1977. Detection of waves converted from P to SV in the mantle, Phys. Earth planet. Inter. , 15( 1), 39– 45. https://doi.org/10.1016/0031-9201(77)90008-5 Google Scholar CrossRef Search ADS   Vinnik L., Kiselev S., Weber M., Oreshin S., Makeyeva L., 2012. Frozen and active seismic anisotropy beneath southern Africa, Geophys. Res. Lett. , 39( 8), doi:10.1029/2012GL051326. https://doi.org/10.1029/2012GL051326 Wang X., Chen Q.F., Li J., Wei S., 2016. Seismic sensor misorientation measurement using P-wave particle motion: an application to the NECsaids Array, Seismol. Res. Lett. , 87( 4), 901– 911. https://doi.org/10.1785/0220160005 Google Scholar CrossRef Search ADS   Webb S.C., 1998. Broadband seismology and noise under the ocean, Rev. Geophys. , 36( 1), 105– 142. https://doi.org/10.1029/97RG02287 Google Scholar CrossRef Search ADS   Wessel P., Smith W.H.F., 1998. New, improved version of the generic mapping tools released, EOS, Trans. Am. geophys. Un. , 79, 579. Google Scholar CrossRef Search ADS   Yoshizawa K., Yomogida K., Tsuboi S., 1999. Resolving power of surface wave polarization data for higher-order heterogeneities, Geophys. J. Int. , 138, 205– 220. https://doi.org/10.1046/j.1365-246x.1999.00861.x Google Scholar CrossRef Search ADS   Zha Y., Webb S.C., Menke W., 2013. Determining the orientations of ocean bottom seismometers using ambient noise correlation, Geophys. Res. Lett. , 40, 3585– 3590. https://doi.org/10.1002/grl.50698 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Table S1. Replacement date of velocity sensor of KG network in South Korea. Table S2. Comparison of the orientations based on different methods. Figure S1. Synthetic test results showing means and five harmonic terms of both R- and T-RFs in the case of full (a) and poor (b) backazimuthal coverage of earthquakes. We use the model with the 3 km thick topmost anisotropic layer (Fig. 3d) to compute the R- and T-RFs (left panel), and their means and harmonic terms (right panel). The RF is sampled with an interval of 5° in backazimuth. Figure S2. Distribution of teleseismic earthquakes in spatial and temporal scales, recorded from the OBS J61C (Cascadia Initiative; Toomey et al. 2014). (a) The numbers of earthquakes for P and PP phases are 276 and 95, respectively, in 2013–2014. The location of the station is indicated by a red rectangle. (b) Cumulative coverage in backazimuth ray coverage for nine months, discretized in 72 bins, after September 2013. (c) Distribution of the backazimuth and slowness of the earthquake in (b). Figure S3. Orientation estimations using various half-widths in Gaussian low-pass filter for the station HSB. (a) Period of January 2005−January 2007 and (b) period of January 2008−January 2015. The infinity symbol means no filter. Figure S4. Difference between the orientation estimates from two methods, which are (1) the minimization of HT1 and (2) the minimization of the mean of T-RF near 0 s (Table S2). The stations SEHB, INCN00 and KOHB are excluded in this figure because of large differences (2.3°, 7.6° and 4.6°, respectively) in this figure. Figure S5. Stacked T-RFs and temporal change of their maximum and minimum amplitudes for stations HSB, MGB and SEHB. (a, c, e) Stacked T-RFs (lower panel) and a histogram showing the number of earthquakes in each bin with a two-month window (upper panel) for stations HSB, MGB and SEHB. Blue arrows indicate reported replacement date of sensors (Table S1). (b, d, f) Temporal change of the maximum and minimum amplitudes (black and blue lines, respectively) of the T-RFs within −1.2 to 1.2 s for stations HSB, MGB and SEHB. Coloured background indicates the number of earthquakes in time and backazimuth. The fluctuation of either maximum or minimum amplitude can be partially explained by heterogeneous space–time occurrence of earthquakes. Grey background colour means no data. Figure S6. Temporal change of the maximum and minimum amplitudes (black and blue lines, respectively) of the T-RF within −1.2 to 1.2 s for 51 stations (including the station INCN ‘00’ and ‘10’). See Figs S5b, d and f for stations HSB, MGB and SEHB. Coloured background indicates the number of earthquakes in time and backazimuth. A blue arrow indicates reported replacement date of a sensor (Table S1). The fluctuation of either maximum or minimum amplitude can be partially explained by heterogeneous space–time occurrence of earthquakes. Grey background colour means no data. Figure S7. Harmonic terms of R- and T-RFs, and their means for the station HSB during the two different validation periods. In each panel, first two columns show the results prior to the correction, and the third and fourth columns after the correction. Figure S8. Harmonic terms of R- and T-RFs, and their means for the station MGB during the three different validation periods. In each panel, first two columns show the results prior to the correction, and the third and fourth columns after the correction. Figure S9. Harmonic terms of R- and T-RFs, and their means for the station SEHB during the two different validation periods. In each panel, first two columns show the results prior to the correction, and the third and fourth columns after the correction. Figure S10. Stacked T-RFs plotted according to backazimuths for stations HSB, MGB and SEHB for each different validation period, prior to (a, c, e, g, i, k, m) and after the orientation correction (b, d, f, h, j, l, n). A solid green horizontal line marks a time of 0 s. The histogram in upper panel shows the number of earthquakes in each bin with a 5° window. The right panel shows a constant harmonic term of the T-RF (blue traces indicating not corrected waveforms). Figure S11. The function f (φ) for OBSs (a) J61C and (b) J39C plotted against φ. Note that the range of the function f (φ) is quite limited because of small amplitudes of the T-RFs at and near 0 s for the OBS data. Regardless of such small range in f (φ), the φmin can be determined as shown in the figure. Figure S12. Backazimuthal coverage plotted against the cumulative numbers of earthquakes for (a) the station JJB and (b) OBS J61C. Figure S13. Histogram showing our uncertainties and differences in the orientation estimates shown in Table 1. (a) The error based on the bootstrapping method (the third column in Table 1). (b) Difference between the orientations determined by the minimization of HT1 (the second column) and the mean of T-RFs (the fifth column) (see also Fig. S4). Difference between ours (the second column) and previous estimates, (c) Shin et al. (2009) (the eighth column), (d) Lee & Sheen (2015) (the ninth column), (e) Lee & Rhie (2015) (the tenth column), respectively. The differences shown in (b), (c), (d) and (e) are root-mean-square values, and we exclude a few values which exceed 10°. Dif.—difference; col.—column; stddev—standard deviation. Figure S14. Synthetic test results showing means and harmonic terms of both R- and T-RFs in the cases of (a and b) for Korean seismic network and (c and d) the Cascadia Initiative when the RFs are stacked by 5° of bin in backazimuth (a and c) and not stacked (b and d). Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial