Numerical Modeling on Ocean-Bottom Seismograph P-Wave Receiver Function to Analyze Influences of Seawater and Sedimentary Layers

Abstract

It is challenging to apply the receiver function method to teleseisms recorded by ocean-bottom seismographs (OBSs) due to a specific working environment that differs from land stations. Teleseismic incident waveforms reaching the area beneath stations are affected by multiple reflections generated by seawater and sediments and noise resulting from currents. Furthermore, inadequate coupling between OBSs and the seabed basement and the poor fidelity of OBSs reduce the signal-to-noise ratio (SNR) of seismograms, leading to the poor quality of extracted receiver functions or even the wrong deconvolution results. For instance, the poor results cause strong ambiguities regarding the Moho depth. This study uses numerical modeling to analyze the influences of multiple reflections generated by seawater and sediments on H-kappa stacking and the neighborhood algorithm. Numerical modeling shows that seawater multiple reflections are mixed with the coda waves of the direct P-wave and slightly impact the extracted receiver functions and can thus be ignored in subsequent inversion processing. However, synthetic seismograms have strong responses to the sediments. Compared to the waveforms of horizontal and vertical components, the sedimentary responses are too strong to identify the converted waves clearly. The extracted receiver functions correspond to the above influences, resulting in divergent results of H-kappa stacking (i.e., the Moho depth and crustal average V_P/V_S ratio are unstable and have great uncertainties). Fortunately, waveform inversion approaches (e.g., the neighborhood algorithm) are available and valid for obtaining the S-wave velocity structure of the crust–upper mantle beneath the station, with sediments varying in thickness and velocity.

1. Introduction

When conducting a seismic observation using ocean-bottom seismographs (OBSs) to detect Earth’s deep structures, the signals are mainly from onboard air guns (active sources) and earthquakes (passive sources). Passive sources usually have massive power, so they are able to reveal much deeper and farther structures reaching the inner core with lower cost [1]. The receiver function (RF) is widely used in passive seismic observations and is an effective way to obtain the crust–upper mantle structures employing data recorded at a single seismic station [2,3,4,5]. Due to the approximately vertical incidence of the ray path of a teleseismic event, the vertical component waveform is roughly assumed to be the source function. Seismic responses from near station structures can be extracted by performing deconvolution in frequency [6] or time [7] domains to eliminate the source function from the radial or tangential component. This result is the so-called RF. By applying post-processing techniques to the RFs, deep structures below the station(s) can be obtained. The H-kappa stacking (H-k stacking; [8]) and Common Converted Point Stacking (CCP Stacking; [9]) are used to calculate the Moho depth, average crustal V_P/V_S ratio and main discontinuities in Earth’s interior. Inversion techniques (e.g., the linear inversion approach, genetic algorithm and simulated annealing; [10,11,12]) are utilized to derive a one-dimensional (1-D) S-wave velocity (V_S) structure.

Even though the RF method is widely applied to data from land stations, it is still hard to use in marine seismic observations. Firstly, an OBS is deployed on the seafloor, so oceanic currents and other noise from complex sources (such as tides, electric current signals of an OBS and signals produced by tilted instruments) reduce the waveform quality recorded by the OBS. Secondly, when teleseismic waves arrive at the seafloor, the wave propagation behaviors generate transmissions and reflections. The up and down multiple reflections in the seawater layer are recorded by an OBS mixed with the converted waves of the direct P-wave. Thirdly, when deploying an OBS, it free-falls through seawater and then lands on the seafloor, which leads to insufficient contact and poor coupling between the OBS and the seabed basement. Lastly, conversion from shear waves to Scholte waves near the OBS and sometimes by the OBS itself results in poor vector fidelity, even when the OBS is well coupled [13]. According to the above-mentioned issues, deriving effective RFs from OBS passive seismic data is still a great challenge. Even for data recorded by the same seismic station, SNRs change over time periods. Hu et al. [2] proposed that suppressing non-stationary noise is critical for extracting RFs from OBS recordings by analyzing huge amounts of data. A previous study on data recorded by seven OBSs deployed at the ultra-spreading Southwest Indian Ridge shows that only one OBS successfully extracted three RFs from a month of recordings [5]. When comparing the results of the high-quality RFs from [2,5] with the results from active source inversions [14,15], the Moho depth determined by the RF method is 2–3 km deeper than that from the active source inversion. The differences are interpreted as indicating that the Moho discontinuity is a transition zone with a certain thickness [16,17,18,19,20]. When Moho-reflected waves of an active source reach the upper interface of this transition zone, they will return immediately. However, the plane waves from teleseismic incidents piercing the upper mantle reach the lower Moho interface, generating the converted waves at once. The above interpretation is a reliable case for the different Moho depths of RFs and active source inversions. Alternatively, different seismic inversion algorithms have different parameterization schemes and constraints, which may also lead to slight variabilities in the final results. For example, when using OBS RFs to obtain the lithosphere structures in the South China Sea Basin, the Moho depths estimated by the H-k stacking method and the neighborhood algorithm (NA) [21,22,23] were quite similar beneath some stations, but there are a few differences below other stations [2].

Previous studies have been carried out to try to solve parts of the above-mentioned issues. Audet [24] used a series of geological models to compute their seismic responses and RFs and then analyzed which factors affect the OBS RFs, such as the noise of seawater and the reflections from the seawater or sedimentary layer. The results found that the sedimentary layer leads to unstable OBS RFs. Akuhara et al. [25] used an advanced non-linear waveform analysis technique accompanied by the simulated annealing algorithm to analyze OBS array recordings from the water layer response. Yang et al. [26] determined the influence of the sedimentary layer on OBS RFs through numerical modeling to constrain gravity inversion. However, these studies have not completely resolved the challenges associated with OBS RFs. Here, we perform some numerical modeling focusing on the seawater layer with different thicknesses and the sedimentary layer with both different thicknesses and a V_S that affect the inversion results of the widely used H-k stacking method and the NA. Because the models are composed of a horizontally stratified medium, we compute synthetic seismograms by using the wavenumber integration method [27] based on modules from Computer Programs in Seismology (CPS330; [28]). To enhance the SNRs of the RFs, we use the time-domain iterative deconvolution method [7] to extract the RFs. Finally, we use the H-k stacking method to estimate the Moho depth and average crustal V_P/V_S ratio and use the non-linear NA to invert the V_S structure beneath the station.

2. RF-Related Methods

2.1. Receiver Function

For a distant earthquake, the incident angle near the station is approximately vertical. Due to the close ray path at the source side, the influence of the source side on different components of one seismograph could be ignored. Therefore, the recorded three-component waveforms can be treated as the convolution of the instrument response, the source function and the medium response beneath the station, which can be expressed in the time domain:

(1)$\left\{\begin{array}{c}{D}_V\left(t\right)=I\left(t\right) \ast S\left(t\right) \ast E_V\left(t\right)\\ D_R\left(t\right)=I\left(t\right) \ast S\left(t\right) \ast E_R\left(t\right)\\ D_T\left(t\right)=I\left(t\right) \ast S\left(t\right) \ast E_T\left(t\right)\end{array}\right.$

where $S\left(t\right)$ represents the effective source time function of the incident plane wave; $I\left(t\right)$ represents the instrument response; and $E_V\left(t\right)$, $E_R\left(t\right)$ and $E_T\left(t\right)$ represent the medium structure responses in the vertical, radial and tangential components, respectively. Under the equivalent source assumption of the vertical component [29], $E_V\left(t\right)$ is suggested to be an impulse function. Therefore, the two horizontal components perform deconvolutions with the vertical component in the frequency domain:

(2)$\left\{\begin{array}{c}{E}_R\left(\omega \right)={\displaystyle \frac{D_R\left(\omega \right)}{I\left(\omega \right)S\left(\omega \right)}}\approx {\displaystyle \frac{D_R\left(\omega \right)}{D_V\left(\omega \right)}}\\ E_T\left(\omega \right)={\displaystyle \frac{D_T\left(\omega \right)}{I\left(\omega \right)S\left(\omega \right)}}\approx {\displaystyle \frac{D_T\left(\omega \right)}{D_V\left(\omega \right)}}\end{array}\right.$

$E_R\left(t\right)$ and $E_T\left(t\right)$ can be obtained by the inverse Fourier transform of $E_R\left(\omega \right)$ and $E_T\left(\omega \right)$, respectively. These are so-called radial and tangential RFs. Different from the deconvolution algorithm in the frequency domain, the time-domain iterative deconvolution method based on the assumption of the equivalent source minimizes the difference between the observed horizontal component and a predicted signal generated by the convolution of the iteratively predicted pulse time series with the vertical component to obtain the final RF. As depicted in Figure 1, two kinds of RFs (P- and S-wave RFs) are based on their incident wave type. Using their post-processing, the Moho depth, average crustal V_P/V_S ratio and V_S structure of the medium beneath the station can be estimated.

2.2. H-k Stacking

The H-k stacking method is widely used to determine the Moho depth (H) and the average crustal V_P/V_S ratio (k) [8]. When a certain average crustal P-wave velocity (V_P) beneath the station is given, the travel-time differences between multiple phases (i.e., Ps, PpPs and PpSs + PsPs) and the direct P-wave can be calculated. For example, the relationship of the Ps phase between the travel time, velocity and thickness is expressed by the following Equation (3).

(3)$\mathrm{H}={\displaystyle \frac{{\mathrm{t}}_{\mathrm{P}\mathrm{s}}}{\sqrt{\frac{1}{{\mathrm{v}}_{\mathrm{S}}^2}-{\mathrm{p}}^2}-\sqrt{\frac{1}{{\mathrm{v}}_{\mathrm{P}}^2}-{\mathrm{p}}^2}}}$

where H represents the Moho depth; ${\mathrm{t}}_{\mathrm{P}\mathrm{s}}$ is the travel-time difference between the Ps phase and the direct P-wave; V_S is the average crustal S-wave velocity (can be represented by ${\mathrm{V}}_{\mathrm{P}}/\mathrm{k}$); and $\mathrm{p}$ is the ray parameter. According to the estimation of [8], H is much more sensitive to the variation in k than to that in V_P. By grid searching for H and k in regular grids; calculating the travel times of Ps, PpPs and PsPs + PpSs; and stacking the corresponding amplitudes on RFs, the optimal value of H and k can be determined according to the maximum stacked amplitude:

(4)$\mathrm{s}\left(\mathrm{H},\mathrm{k}\right)={\mathrm{W}}_1\mathrm{r}\left({\mathrm{t}}_{\mathrm{P}\mathrm{s}}\right)+{\mathrm{W}}_2\mathrm{r}\left({\mathrm{t}}_{\mathrm{P}\mathrm{p}\mathrm{P}\mathrm{s}}\right)-{\mathrm{W}}_3\mathrm{r}\left({\mathrm{t}}_{\mathrm{P}\mathrm{p}\mathrm{S}\mathrm{s}+\mathrm{P}\mathrm{s}\mathrm{P}\mathrm{s}}\right)$

where $\mathrm{s}\left(\mathrm{H},\mathrm{k}\right)$ represents the stacked amplitude with the specific H and k; ${\mathrm{W}}_1$, ${\mathrm{W}}_2$ and ${\mathrm{W}}_3$ are the stack weights for the different multiple phases (the weights used in this study are 0.6, 0.3 and 0.1, respectively); and $\mathrm{r}$ is the amplitude of the RF waveform corresponding to the travel times.

2.3. Neighborhood Algorithm

With the development of computer techniques, non-linear methods have been rapidly used in geophysical inversions, such as the genetic algorithm [12], simulated annealing [11] and Neural Network [30]. These non-linear inversion algorithms can search model parameters globally to overcome the inherent limitations of linear approximation, avoid solutions falling into local extremum and significantly reduce the reliance on the initial model. However, these algorithms contain a lot of empirical and random parameters, which lead to complicated inversion processes and huge computation costs. Here, we use the NA proposed by Sambridge [21,22,23] for the inversion of RFs. This algorithm parameterizes the model randomly. The best velocity model can be obtained by minimizing the waveform difference between synthetic and observed RFs, which is able to infer the geological information beneath the station. Similar to the genetic algorithm, the NA can search the whole model space globally. However, the NA incorporates the geometric concept of the “Voronoi diagram”, which offers robust and adaptive searching capabilities. Initially, the model space is partitioned into numbers of Voronoi cells by randomly sampled points. Subsequently, mismatch functions are constructed within each Voronoi cell to find a series of cells with acceptable mismatch values. The reliable model is inverted by sparsely and densely resampled cells with relatively higher and lower mismatches.

3. Seismic Responses of Seawater Multiple Reflections

Differing from land stations, OBSs are deployed on the seafloor. Due to the strong reflection boundaries of the free surface and seafloor, incident P-waves from distant events generate multiple reflections in the seawater layer. Because the incident waves are not exactly perpendicular to the horizontal components, when the vertical component records these multiples, their polarized waves are recorded by the horizontal components simultaneously. Therefore, although these added multiples destroy the assumption of the equivalent source of the vertical component, the highly correlated polarized waves on the horizontal components can mitigate these influences on the RFs. The following numerical modeling is conducted to identify seawater multiple reflections and analyze their influences to verify the above inference.

3.1. Half-Space Model with Overlying Seawater Layer

The average thickness of the oceanic crust is approximately 6 km (hereafter, it is referred to as the normal oceanic crust) [31,32,33], but an extremely thick oceanic crust is found at the Southwest Indian Ridge, where the crust thickness reaches 10 km (we call it the thick oceanic crust hereafter) [34]. Therefore, two models with thicknesses of 6 km and 10 km, respectively, are considered in the synthetic modeling. The thickness of the seawater layer varies from 0 to 5 km with a step of 1 km. The parameters of each layer of the 1-D stratified model are listed in Table S1, and the velocity structures are illustrated in Figure 2a and Figure 3a.

3.2. Synthetic Results of Normal Oceanic Crust Models

Figure 2b–d illustrate the synthetic vertical and radial component seismograms and their corresponding RFs for different seawater layer thicknesses. Notably, the seismic phases of P, Ps, PpPs and PsPs + PpSs can be clearly identified in the RF waveforms (Figure 2d), which are almost not affected by the thickness of the seawater layer. In terms of the vertical (Figure 2b) and radial (Figure 2c) seismograms and RFs (Figure 2d), waveforms within 0–5 s have very high consistency. In the absence of a seawater layer, the RF waveforms after 5 s are nearly zero, without other phases. However, when the models have a seawater layer, some seawater multiple reflections are contained in the RF waveforms after 5 s and result in alternating positive and negative waveforms, respectively. Fortunately, their amplitudes are significantly lower than those of the Moho reflections. Furthermore, three obvious anomalies are found in Figure 2d: (1) when the thickness of the seawater layer is 1 km, the RF has a negative peak around 1 s, and the amplitude of the PpPs phase is relatively stronger than that of the seawater layer thickness equal to 0 km; (2) when the seawater layer thickness is 2 km, the amplitude of the PpPs phase is comparatively the smallest; and (3) when the seawater layer thickness is 3 km, there is a negative peak following the PsPs + PpSs phases, which is considerably stronger than no seawater layer. These three anomalies are obviously a result of seawater multiples.

3.3. Analysis of Abnormal Characteristics of Synthetic Results

To interpret the above-mentioned anomalies of RFs, we analyze the travel-time characteristics of seawater multiple reflections on horizontal and vertical components. The vertical component shows that the odd-numbered waveforms have positive peaks (see the blue solid arrows in Figure 2b), and the even-numbered waveforms have negative peaks (see the blue dashed arrows in Figure 2b), which are consistent with the theoretical predictions of seawater multiple reflections. At the same time, corresponding positive and negative peaks can be found on the radial component at many places. However, most places are hard to identify visually because polarized waves are too weak and mixed with other seismic waves (Figure 2c). Compared to a synthetic RF calculated from the model without a seawater layer, the places differing from those of models with seawater layers are consistent with the predicted reflections. Therefore, we conclude that the seawater layer has certain influences on OBS RFs, but the influences are quite slight. Generally, for a thin seawater layer (e.g., 1 km), the travel times of multiple reflections are relatively short; thus, they appear within 2.5–10 s. However, when the seawater layer is thick (e.g., 5 km), multiple reflections are not visible before 5 s, and the RF is the same as the RF calculated from the model without a seawater layer in this period. After 5 s, there are alternatively positive and negative phases. Whatever the thickness variation in the seawater layer, the influences are reduced with time. Their impacts on the RFs can be ignored after 20 s. In terms of the seawater thickness, a thinner seawater layer has a relatively stronger impact on RFs than a thicker seawater layer.

3.4. Influence of Seawater Multiple Reflections on H-k Stacking Inversion

Theoretically, the travel times of the reflected seismic waves from the Moho discontinuity are less than 5 s on the synthetic RFs, and seawater multiple reflections would affect the application of the H-k stacking method to the inversion of the Moho depth and average crustal V_P/V_S ratio when the seawater layer is thin. However, the phases of Ps, PpPs and PsPs + PpSs dominate the main waveforms of the RFs; thus, the disturbances induced by seawater multiple reflections are quite slight. According to Figure 2, it is reasonable to deduce that the influence of seawater multiple reflections on H-k stacking inversion could be ignored, and the corresponding results are shown in Figure 4. Based on the RFs calculated from the normal oceanic crust model (the oceanic crust thickness and the average crustal V_P/V_S ratio are 6 km and 1.714, respectively), the H-k stacking method can be applied to these RFs to achieve the results. Comparing Figure 4a with Figure 4b, the calculated crust thickness equals the presupposed thickness (6 km), and the computed V_P/V_S ratio (1.71) is close to the presupposed value (1.714). The calculated crust thicknesses are the same whether the seawater thickness is 0 km or 5 km. In addition, because the grid-searching step of the V_P/V_S ratio is 0.01, this tiny difference has no significance in synthetic tests and real data practices. It is worth noting that the grid search of the H-k stacking method should exclude the initial arrival time to avoid the strong amplitude of the direct P-wave, which is not used in stacking. We have to adjust the search range, especially when the optimal point is close to the grid edge.

3.5. Synthetic Results of Thick Oceanic Crust Models

To analyze the impact of seawater multiple reflections on the RF in the case of a thick oceanic crust model (10 km), we also calculated the corresponding RFs (Figure 3d) and analyzed the influence of seawater multiple reflections on the application of the H-k stacking method (Figure 5). They are extremely similar to those of normal oceanic crust models (Figure 2 and Figure 3). The waveforms of the vertical, radial and RF components are close to those of the normal oceanic crust models. Because of the thicker oceanic crust, the longer travel times of the Moho-reflected waves extend the time window of seawater multiple reflections up to ~10 s (Figure 3d), in which the Ps, PpPs and PsPs + PpSs phases are affected by seawater multiple reflections. However, the P, Ps, PpPs and PsPs + PpSs phases dominate the main waveforms (Figure 3d), which are consistent with the results of the normal oceanic crust models. Although seawater multiple reflections still have impacts on the RFs of the thick oceanic crust models, they do not affect the application of the H-k stacking method (see Figure 5). Compared with Figure 2d, due to the thicker oceanic crust, the converted phases of the Moho are distributed sparsely on the waveforms, which simplifies the selection of the H-k stacking parameters.

4. Influences of Sedimentary Layer Thickness on OBS RFs

Except for the influence of the seawater layer, oceanic basins are covered by marine sediments, especially in marginal seas. The seismic wave propagation and dissipation in sedimentary layers may directly affect the quality of the RFs extracted from the observed seismic data. Therefore, it is necessary to evaluate the impact of sedimentary layers for most oceanic basins, except for regions of newly formed oceanic crust or seamounts.

4.1. Half-Space Model with Seawater and Sedimentary Layers

To estimate the detailed impact of a sedimentary layer on OBS RFs, the model parameters of the synthetic tests are configured in Table S2. The velocity structures of the normal and thick oceanic crust are shown in Figure 6a and Figure 7a, respectively. According to Table S2, these models set the seawater thickness to 5 km. The normal and thick oceanic crust thicknesses are 6 km and 10 km, respectively. The thickness of the sedimentary layer varies from 0 to 2.5 km with a step of 0.25 km. The V_P and V_S of the sediments are 2.1 km/s and 0.8 km/s, respectively. The density of the sedimentary layer is 2.0 g/cm³.

4.2. Influences of Sedimentary Layers on OBS RFs

The synthetic seismograms and their RFs are shown in Figure 6b–d and Figure 7b–d, respectively. The figures show that the thickness variation in the sedimentary layer would drastically affect the waveform of the RFs. When the sedimentary layer is thin, the converted waves (Pbs; Figure 6d and Figure 7d blue dashed line) of the layer have a larger amplitude than that of the direct P-wave (Figure 6d and Figure 7d black dashed line), and their arrival times are close. The thin sedimentary layer with a relatively lower speed causes the delayed arrival of the direct P-wave and results in the overlap of their waveforms on the RFs. When the thickness of the sedimentary layer increases, the amplitude of the direct P-wave becomes weaker, and the phases of Ps, PpPs and PsPs + PpSs are difficult to discriminate distinctly. Multiple reflected phases overlap with each other, leading to the Ps, PpPs and PsPs + PpSs phases not occupying a significant dominant position on the RFs. These phenomena indicate that submarine sediments significantly influence the effective information of RFs, which correspond to the geological structures below the seismic station array.

4.3. Influences of Sedimentary Layers on H-k Stacking Inversion

The synthetic results of the models with sediments indicate that the seismic responses are quite complicated compared to the models without sedimentary layers. The waveform overlap of multiple converted phases reduces or even destroys the convergence of the H-k stacking results, which leads to difficulty in determining the oceanic crust thickness and average crustal V_P/V_S ratio accurately (see Figure 8 and Figure 9). When the oceanic crust is 6 km thick, thickness variations in the sedimentary layer lead to the estimated oceanic crust being too thick or too thin. Most of them are thicker than the synthetic models. Similarly, thickness variations in the sedimentary layer also result in lower or higher V_P/V_S ratios (see Table S3 for more details). When the oceanic crust is 10 km thick, the thickness variations in the sedimentary layer cause the same results (see Figure 10 and Table S3). These tests suggest that the existence of a sedimentary layer causes a significant difference in structures between the H-k stacking results and the presupposed model, and the variation in the sedimentary layer leads to a much more complicated situation.

4.4. Influences of Sedimentary Layers on RFs Waveform Inversion

The above numerical tests show that the sedimentary layer (regardless of its thickness) will lead to significant uncertainties in the estimated structures by the H-k stacking method due to the overlap and disordered waveform of multiple reflected phases. Therefore, the effectiveness of the inversion using RF waveforms to obtain the V_S structure beneath the station through non-linear inversion algorithms becomes very important, because the V_S structure can serve as an alternative to the H-k stacking results.

This study adopts the NA in RF waveform inversions for different models, including thickness variations in the oceanic crust, seawater and sedimentary layer (see Table S3 for details). Figure 11 illustrates two results obtained from the NA inversion for the 6 and 10 km thick crust with overlying 0.5 km thick sediments, respectively. The left panel (Figure 11a) shows the results of the normal oceanic crust (6 km) with a sedimentary layer of 0.5 km (i.e., the Moho depth is 6.5 km). The inverted optimal V_S and V_P/V_S ratio are very close, with very slight differences. In contrast, when the oceanic crust thickness is 10 km with a sediment layer of 0.5 km thick (i.e., a Moho depth of 10.5 km), the results of the NA inversion are approximately the same as the presupposed model. In addition, the other inversion results based on the parameters listed in Table S3 are also perfect. In practical work, the seawater depth can be obtained using acoustic techniques (e.g., multi-beam systems), and the V_P of the seawater can be taken as a constant of 1500 m/s. Therefore, the synthetic tests mentioned above suggest that the influence of the seawater layer on the inversion results is almost negligible. Furthermore, the RF waveform inversion (such as the NA used in this study) can effectively overcome the negative effects caused by the sedimentary layer. Thus, the reliable structures below the stations are available by using the post-processing of the RF mentioned above.

5. Influences of V_S Variation in the Sedimentary Layer on OBS RFs

5.1. Half-Space Model with Varied V_S of the Sedimentary Layer

Different from the almost constant seismic velocity values of seawater, both the seismic velocity and thickness of the sedimentary layer should be considered in the RF waveform inversion. Following, we simulate the influence of V_S variations in the sedimentary layer on OBS RFs. According to Table S4, these models set the seawater and sediment thickness to 5 km and 0.5 km, respectively. The sediment V_S varies from 0.1 to 1.0 km/s with a step of 0.1 km/s, the V_P is always set to 2.1 km/s and the density is 2.0 g/cm³. Table S4 presents the parameters for each layer of the presupposed models, and the velocity structures are illustrated in Figure 12a and Figure 13a.

5.2. Influences of the Sedimentary Layer with Different Vs on RFs and H-k Stacking Inversion

Using the same synthetic test routines as in Section 3.2 and Section 3.3, we employ models with the parameters shown in Table S4 to calculate the synthetic seismograms for extracting theoretical RFs. Then, these RFs will be used to test the H-k stacking method and RF waveform inversion. Figure 12 and Figure 13 show the theoretical seismogram and RFs computed from the models with different Vs values for the sedimentary layer for both normal and thick oceanic crust thicknesses (6 km and 10 km thick, respectively). Similar to Figure 6 and Figure 7, the sedimentary layer leads to waveform overlap among the converted multiple phases. When the Vs decreases, the arrival times of the main converted phases of the RFs are delayed, causing a significant deviation between the H-k stacking results and the presupposed models (see Figure 14 and the detailed information in Table S5). In general, no matter how the Vs changes, the H-k stacking results of the Moho depth and crustal V_P/V_S ratio are not ideal and are quite different from the presupposed values.

5.3. Influence of the Sedimentary Layer Different Vs on the RF Waveform Inversion

To further analyze the impact of Vs variations in the sedimentary layer on the RF waveform inversion and to verify the effectiveness of the waveform inversion, we adopted the NA to estimate the model parameters by using the RFs obtained under different Vs values for the sedimentary layer. In the presupposed models, the thicknesses of the oceanic crust are 6 km and 10 km, the seawater layer is 5 km thick and the sedimentary layer is 0.5 km thick, with its Vs varying from 0.1 km/s to 1.0 km/s. Figure 15 shows the inversion results of the normal and thick oceanic crust models (6 km and 10 km, respectively), in which the Vs of the sedimentary layer is 0.4 km/s. It is shown that the inversion results are aligned with that of the presupposed models, and the RF waveform fittings are also acceptable. These simulations indicate that changes in the Vs of the sedimentary layer have no significant impact on the RF inversion when the NA algorithm is used.

6. Conclusions

This paper focuses on some issues in applying the RF method to passive seismic observation data of OBSs. Numerical tests are performed to estimate the influences of seawater multiple reflections and the sedimentary layer on the RFs. Two post-processing methods (H-k stacking and the NA) are used to invert the Moho depth, average crustal V_P/V_S ratio and V_S structure beneath the station. The main conclusions are as follows:

(1) Seawater multiple reflections would break the equivalent source assumption of the vertical component, but the deconvolution of the radial and vertical components helps to reduce the influences caused by seawater multiples. This is because the polarization waves on the radial component are closely correlated with multiple reflections on the vertical component. As a result, the negative effect of seawater multiple reflections on RFs can be significantly eliminated or almost negligible through deconvolution. Consequently, the H-k stacking method and the non-linear algorithm (e.g., the NA) can be adopted to estimate an accurate Moho depth, an average crustal V_P/V_S ratio and the V_S of the crust–upper mantle structure.

(2) The sedimentary layer significantly impacts synthetic seismograms, and the dispersion effects of seismic wave propagation in sediments would cause challenges in extracting RFs from OBS recordings. The first arrival and multiple reflections of the Moho are weakened, and it is difficult to identify them with vast variations in their waveforms. These changes lead to significant uncertainties in estimating the Moho depth and the crustal average V_P/V_S ratio when using the H-k stacking method.

(3) The waveform-related inversion methods (e.g., the NA) can be effectively applied to OBS RFs to obtain the V_S structures from the crust to the upper mantle beneath the station whenever the sedimentary layer varies in thickness and velocity.

Author Contributions

Conceptualization, W.G., H.H. and A.R.; methodology, W.G. and H.H.; validation, W.G.; formal analysis, W.G. and H.H., investigation, W.G.; writing—original draft preparation, W.G. and H.H.; writing—review and editing, W.G., A.R., H.H., X.N., W.W. and Y.T.; visualization, W.G.; supervision, A.R., H.H. and Y.T. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data will be made available on request.

Acknowledgments

The forward modeling and calculating the receiver functions use the Computer Programs in Seismology (CPS330; [28]). Performing the H-k stacking uses the codes provided by Lupei Zhu [8]. The NA sampler program [21,22,23] is used to invert the 1-D Vs structure. Some figures in this study were plotted by Generic Mapping Tools (GMT; [35]).

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Illustration of P- and S-wave receiver functions. Different colors correspond to different ray paths of multiple phases. (a) P-wave receiver function and the related ray paths; (b) the same as (a) but for S-wave receiver function.

Figure 2. Normal oceanic crust models and their synthetic seismograms and RFs. (a) A layered velocity model for calculating synthetic seismograms and corresponding RFs. (b,c) Vertical and radial seismograms are calculated from different seawater thicknesses. (d) Synthetic RFs are computed from the waveforms shown in (b,c). The blue solid and dashed arrows represent positive and negative peaks caused by seawater multiple reflections.

Figure 3. The same as Figure 2 but for models with a thick oceanic crust. (a) A layered velocity model for calculating synthetic seismograms and corresponding RFs. (b,c) Vertical and radial seismograms are calculated from different seawater thicknesses. (d) Synthetic RFs are computed from the waveforms shown in (b,c). The blue solid and dashed arrows represent positive and negative peaks caused by seawater multiple reflections.

Figure 4. H-k stacking results of the normal oceanic crust models, where k is the average crustal VP/VS ratio; HM is the Moho depth; and HW is the thickness of the seawater layer; the blue dot represents the optimal value of Moho depth and the average crustal VP/VS ratio. (a) H-k stacking results of the normal oceanic crust model without a seawater layer. (b) The same as (a) but for a normal oceanic crust model with a 5 km thick seawater layer.

Figure 5. The same as Figure 4 but for models with a thick oceanic crust, where k is the average crustal VP/VS ratio; HM is the Moho depth; and HW is the thickness of the seawater layer; the blue dot represents the optimal value of Moho depth and the average crustal VP/VS ratio. (a) H-k stacking results of the thick oceanic crust model without a seawater layer. (b) The same as (a) but for a thick oceanic crust model with a 5 km thick seawater layer.

Figure 6. The same as Figure 2 but for models with a water layer (5 km), varying sedimentary layer thicknesses and a normal oceanic crust. (a) A layered velocity model for calculating synthetic seismograms and corresponding RFs. (b,c) Vertical and radial seismograms are calculated from different sedimentary layer thicknesses. (d) Synthetic RFs are computed from the waveforms shown in (b,c).

Figure 7. The same as Figure 2 but for models with a seawater layer (5 km), varying sedimentary layers thicknesses and a thick oceanic crust. (a) A layered velocity model for calculating synthetic seismograms and corresponding RFs. (b,c) Vertical and radial seismograms are calculated from different sedimentary layer thicknesses. (d) Synthetic RFs are computed from the waveforms shown in (b,c).

Figure 8. H-k stacking results of the normal oceanic crust models with a seawater layer (5 km) and sedimentary layers, where k represents the average crustal VP/VS ratio; HM represents the Moho depth; HC represents the presupposed oceanic crust thickness; and HS represents the submarine sedimentary layer thickness; the blue dot represents the optimal value of Moho depth and the average crustal VP/VS ratio. (a) The thickness of the sedimentary layer is 0.25 km. (b) The thickness of the sedimentary layer is 1.75 km.

Figure 9. H-k stacking results of the thick oceanic crust models with a seawater layer (5 km) and sedimentary layers, where k represents the average crustal VP/VS ratio; HM represents the Moho depth; HC represents the presupposed oceanic crust thickness; and HS represents the submarine sedimentary layer thickness; the blue dot represents the optimal value of Moho depth and the average crustal VP/VS ratio. (a) The thickness of the sedimentary layer is 0.25 km. (b) The thickness of the sedimentary layer is 1.75 km.

Figure 10. Comparison of the H-k stacking results between the presupposed and estimated values of the normal and thick oceanic crust model with a seawater thickness of 5 km and different sedimentary layer thicknesses. The presupposed crustal thickness excludes the sediment thickness. (a,c) The crustal thicknesses and their differences under models with the normal oceanic crust (6 km thick) and the thick oceanic crust (10 km thick), respectively; (b,d) the crustal VP/VS ratios and their differences under models with the normal oceanic crust (6 km thick) and the thick oceanic crust (10 km thick), respectively.

Figure 11. Results of the NA inversion structures. In the synthetic tests, the thicknesses of the seawater and sediments are 5 km and 0.5 km, respectively. The VS of the sediments is 0.8 km/s. (a) The inversion results of the normal oceanic crust (6 km thick). The red and blue solid lines represent the best VS structure and VP/VS ratio; the orange dashed line indicates the Moho depth. (b) The same as (a) but for a model with a thick oceanic crust (10 km thick). (c) A comparison of the waveforms between the synthetic RFs (blue lines) and the RFs calculated from the inversion model (red lines). The solid lines correspond to the model with the normal oceanic crust, and the dashed lines correspond to the model with the thick oceanic crust.

Figure 12. The same as Figure 2 but for models with a different Vs of the sedimentary layer. The oceanic crust, seawater and sedimentary layer thicknesses are 6 km, 5 km and 0.5 km, respectively. (a) A layered velocity model for calculating synthetic seismograms and corresponding RFs. (b,c) Vertical and radial seismograms are calculated from different sedimentary layer VS. (d) Synthetic RFs are computed from the waveforms shown in (b,c).

Figure 13. The same as Figure 3 but for models with a different Vs for the sedimentary layer. The oceanic crust, seawater and sedimentary layer thicknesses are 10 km, 5 km and 0.5 km, respectively. (a) A layered velocity model for calculating synthetic seismograms and corresponding RFs. (b,c) Vertical and radial seismograms are calculated from different sedimentary layer VS. (d) Synthetic RFs are computed from the waveforms shown in (b,c).

Figure 14. The same as Figure 10 but for models with a different Vs for the sedimentary layer. The seawater and sedimentary layer thicknesses are 5 km and 0.5 km, respectively. (a,c) The crustal thicknesses and their differences under models with the normal oceanic crust (6 km thick) and the thick oceanic crust (10 km thick), respectively; (b,d) the crustal VP/VS ratios and their differences under models with the normal oceanic crust (6 km thick) and the thick oceanic crust (10 km thick), respectively.

Figure 15. The same as Figure 11 but for models with different properties. The seawater thickness and sedimentary layer thickness are 5 km and 0.5 km, respectively. The Vs of the sedimentary layer is 0.4 km/s. (a) The inversion results of the normal oceanic crust (6 km thick). The red and blue solid lines represent the best VS structure and VP/VS ratio; the orange dashed line indicates the Moho depth. (b) The same as (a) but for a model with a thick oceanic crust (10 km thick). (c) A comparison of the waveforms between the synthetic RFs (blue lines) and the RFs calculated from the inversion model (red lines). The solid lines correspond to the model with the normal oceanic crust, and the dashed lines correspond to the model with the thick oceanic crust.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/jmse12112053/s1, Table S1: Parameters of seawater–crust–mantle models for synthetic OBS RFs; Table S2: Parameters of seawater–sediments–crust–mantle models for synthetic OBS RFs; Table S3: Comparison of H-k stacking results between presupposed and estimated values of normal and thick oceanic crust model with a seawater thickness of 5 km; Table S4: Parameters of seawater–sediments–crust–mantle models with different V_S for synthetic OBS RFs; Table S5: The same as Table S3 but for models with a different V_S for the sedimentary layer.

References

1. Ruan, A. Ocean Bottom Seismic Theory and Application; Science Press: Beijing, China, 2020.

2. Hu, H.; Ruan, A.G.; You, Q.Y.; Li, J.B. Using OBS teleseismic receiver funtions to invert lithospheric structure—A case study of the southwestern subbasin in the South China Sea. Chin. J. Geophys.; 2016; 59, pp. 1426-1434.

3. Huang, H.; Qiu, X.; Zhang, J.; Hao, T. Low-velocity layers in the northwestern margin of the South China Sea: Evidence from receiver functions of ocean-bottom seismometer data. J. Asian Earth Sci.; 2019; 186, 104090. [DOI: https://dx.doi.org/10.1016/j.jseaes.2019.104090]

4. Hung, T.D.; Yang, T.; Le, B.M.; Yu, Y.; Xue, M.; Liu, B.; Liu, C.; Wang, J.; Pan, M.; Huong, P.T. et al. Crustal Structure Across the Extinct Mid-Ocean Ridge in South China Sea from OBS Receiver Functions: Insights Into the Spreading Rate and Magma Supply Prior to the Ridge Cessation. Geophys. Res. Lett.; 2021; 48, e2020GL089755. [DOI: https://dx.doi.org/10.1029/2020GL089755]

5. Ruan, A.; Hu, H.; Li, J.; Niu, X.; Wei, X.; Zhang, J.; Wang, A. Crustal structure and mantle transition zone thickness beneath a hydrothermal vent at the ultra-slow spreading Southwest Indian Ridge (49°39′ E): A supplementary study based on passive seismic receiver functions. Mar. Geophys. Res.; 2017; 38, pp. 39-46. [DOI: https://dx.doi.org/10.1007/s11001-016-9298-8]

6. Clayton, R.W.; Wiggins, R.A. Source shape estimation and deconvolution of teleseismic bodywaves. Geophys. J. R. Astron. Soc.; 1976; 47, pp. 151-177. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1976.tb01267.x]

7. Ligorría, J.P.; Ammon, C.J. Iterative deconvolution and receiver-function estimation. Bull. Seismol. Soc. Amer.; 1999; 89, pp. 1395-1400. [DOI: https://dx.doi.org/10.1785/BSSA0890051395]

8. Zhu, L.P.; Kanamori, H. Moho depth variation in southern California. J. Geophys. Res.; 2000; 105, pp. 2969-2980. [DOI: https://dx.doi.org/10.1029/1999JB900322]

9. Zhu, L.P. Crustal structure across the San Andreas Fault, southern California from teleseismic converted waves. Earth Planet. Sci. Lett.; 2000; 179, pp. 183-190. [DOI: https://dx.doi.org/10.1016/S0012-821X(00)00101-1]

10. Ammon, C.J.; Randall, G.E.; Zandt, G. On the nonuniqueness of receiver function inversions. J. Geophys. Res. Solid Earth; 1990; 95, pp. 15303-15318. [DOI: https://dx.doi.org/10.1029/JB095iB10p15303]

11. Bertsimas, D.; Tsitsiklis, J. Simulated Annealing. Stat. Sci.; 1993; 8, pp. 10-15. [DOI: https://dx.doi.org/10.1214/ss/1177011077]

12. Holland, J.H. Genetic algorithms. Scholarpedia; 2012; 7, 1482. [DOI: https://dx.doi.org/10.4249/scholarpedia.1482]

13. Paffenholz, J.; Docherty, P.; Shurleff, R.; Hays, D. Shear Wave Noise on OBS Vz Data-Part II Elastic Modeling of Scatterers in the Seabed. Proceedings of the 68th Conference and Exhibition; Vienna, Austria, 12–15 June 2006.

14. Zhang, J.; Li, J.; Ruan, A.; Ding, W.; Niu, X.; Wang, W.; Tan, P.; Wu, Z.; Yu, Z.; Wei, X. et al. Seismic Structure of a Postspreading Seamount Emplaced on the Fossil Spreading Center in the Southwest Subbasin of the South China Sea. JGR Solid Earth; 2020; 125, e2020JB019827. [DOI: https://dx.doi.org/10.1029/2020JB019827]

15. Zhao, M.; Qiu, X.; Li, J.; Sauter, D.; Ruan, A.; Chen, J.; Cannat, M.; Singh, S.; Zhang, J.; Wu, Z. et al. Three-dimensional seismic structure of the Dragon Flag oceanic core complex at the ultraslow spreading Southwest Indian Ridge (49°39′ E). Geochem. Geophys. Geosystems; 2013; 14, pp. 4544-4563. [DOI: https://dx.doi.org/10.1002/ggge.20264]

16. Ciazela, J.; Koepke, J.; Dick, H.J.B.; Botcharnikov, R.; Muszynski, A.; Lazarov, M.; Schuth, S.; Pieterek, B.; Kuhn, T. Sulfide enrichment at an oceanic crust-mantle transition zone: Kane Megamullion (23° N, MAR). Geochim. Cosmochim. Acta; 2018; 230, pp. 155-189. [DOI: https://dx.doi.org/10.1016/j.gca.2018.03.027]

17. Herman, J.; Muentener, O.; Guenther, D. Differentiation of mafic magma in a continental crust-to-mantle transition zone. J. Petrol.; 2001; 42, pp. 189-206. [DOI: https://dx.doi.org/10.1093/petrology/42.1.189]

18. Zanon, V.; Silva, R.; Goulart, C. The crust-mantle transition beneath the Azores region (central-north Atlantic Ocean). Contrib. Mineral. Petrol.; 2023; 178, 50. [DOI: https://dx.doi.org/10.1007/s00410-023-02036-6]

19. Song, T.; Shen, X.; Mei, X. Constraining Moho characteristics with frequency-dependence of receiver function and its application. Acta Seismol. Sin.; 2020; 42, pp. 135-150. [DOI: https://dx.doi.org/10.11939/jass.20190149]

20. Kaban, M.K.; Flovenz, O.G.; Palmason, G. Nature of the crust-mantle transition zone and the thermal state of the upper mantle beneath Iceland from gravity modelling. Geophys. J. Int.; 2002; 149, pp. 281-299. [DOI: https://dx.doi.org/10.1046/j.1365-246X.2002.01622.x]

21. Sambridge, M. Geophysical inversion with a neighborhood algorithm—I, Searching a parameter space. Geophys. J. Int.; 1999; 138, pp. 479-494. [DOI: https://dx.doi.org/10.1046/j.1365-246X.1999.00876.x]

22. Sambridge, M. Geophysical inversion with a neighbourhood algorithmö—II. Appraising the ensemble. Geophys. J. Int.; 1999; 138, pp. 727-746. [DOI: https://dx.doi.org/10.1046/j.1365-246x.1999.00900.x]

23. Sambridge, M. Finding acceptable models in nonlinear inverse problems using a neighbourhood algorithm. Inverse Probl.; 2001; 17, pp. 387-403. [DOI: https://dx.doi.org/10.1088/0266-5611/17/3/302]

24. Audet, P. Receiver functions using OBS data: Promises and limitations from numerical modelling and examples from the Cascadia Initiative. Geophys. J. Int.; 2016; 205, pp. 1740-1755. [DOI: https://dx.doi.org/10.1093/gji/ggw111]

25. Akuhara, T.; Mochizuki, K.; Kawakatsu, H.; Takeuchi, N. Non-linear waveform analysis for water-layer response and its application to high-frequency receiver function analysis using OBS array. Geophys. J. Int.; 2016; 206, pp. 1914-1920. [DOI: https://dx.doi.org/10.1093/gji/ggw253]

26. Yang, T.; Xu, Y.; Du, N.; Xu, T.; Cao, D.; Nan, F.; Chu, W.; Liang, C.; Hao, T. Gravity inversion constrained by OBS receiver function reveals crustal structure in Ryukyu Trench. Front. Earth Sci.; 2023; 11, 1187683. [DOI: https://dx.doi.org/10.3389/feart.2023.1187683]

27. Kennett, B.L.N.; Kerry, N.J. Seismic waves in a stratified half space. Geophys. J. Int.; 1979; 57, pp. 557-583. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1979.tb06779.x]

28. Herrmann, R.B. Computer Programs in Seismology: An Evolving Tool for Instruction and Research. Seismol. Res. Lett.; 2013; 84, pp. 1081-1088. [DOI: https://dx.doi.org/10.1785/0220110096]

29. Langston, C.A. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. J. Geophys. Res. Solid Earth; 1979; 84, pp. 4749-4762. [DOI: https://dx.doi.org/10.1029/JB084iB09p04749]

30. Bishop, C.M. Neural networks and their applications. Rev. Sci. Instrum.; 1994; 65, pp. 1803-1832. [DOI: https://dx.doi.org/10.1063/1.1144830]

31. Chen, Y.J. Oceanic crustal thickness versus spreading rate. Geophys. Res. Lett.; 1992; 19, pp. 753-756. [DOI: https://dx.doi.org/10.1029/92GL00161]

32. Christeson, G.L.; Goff, J.A.; Reece, R.S. Synthesis of Oceanic Crustal Structure from Two-Dimensional Seismic Profiles. Rev. Geophys.; 2019; 57, pp. 504-529. [DOI: https://dx.doi.org/10.1029/2019RG000641]

33. Hu, H.; Zhao, D.; Lin, J.; Pilia, S. A Slab Window Beneath North Sumatra Revealed by P-Wave Mantle Tomography. JGR Solid Earth; 2023; 128, e2022JB025976. [DOI: https://dx.doi.org/10.1029/2022JB025976]

34. Niu, X.; Ruan, A.; Li, J.; Minshull, T.A.; Sauter, D.; Wu, Z.; Qiu, X.; Zhao, M.; Chen, Y.J.; Singh, S. Along-axis variation in crustal thickness at the ultraslow spreading Southwest Indian Ridge (50° E) from a wide-angle seismic experiment. Geochem. Geophys. Geosyst.; 2015; 16, pp. 468-485. [DOI: https://dx.doi.org/10.1002/2014GC005645]

35. Wessel, P.; Luis, J.F.; Uieda, L.; Scharroo, R.; Wobbe, F.; Smith, W.H.F.; Tian, D. The Generic Mapping Tools Version 6. Geochem. Geophys. Geosyst.; 2019; 20, pp. 5556-5564. [DOI: https://dx.doi.org/10.1029/2019GC008515]