The popular methods used in modeling the deep crustal and lithospheric structure of the marginal sea and oceanic ridges are gravity and seismic methods,the latter of which is much more reliable and precise. The main techniques are passive source tomography and active source reflection and refraction experiments for the seismic method. The passive source tomography is mainly based on surface waves observation and has low accuracy and resolution due to the lack of stations in the oceanic basin. Deploying broadband ocean bottom seismometers(OBSs)on the seabed for long-term observation of earthquakes is an effective method to study the oceanic crustal and lithospheric structures by accumulating data for 3D tomography of one area or by using a single-station observed teleseismic data to model 1D structures beneath the station. For example,the teleseismic S wave splitting technique and receiver functions are used to invert the lithospheric anisotropy and the S wave velocity structure beneath the station,respectively. In recent years,broadband OBS observations in the sea basin in the South China Sea(SCS)have been conducted,resulting in the accumulation of certain volume data and achievement of some concrete results(Ruan et al.,2012; Liu et al.,2014). However,a teleseismic receiver functions inversion method has not been popularly used in lithospheric structures. It is documented that a teleseismic receiver functions inversion method based on OBS passive observation data has scarcely been used or mentioned by researchers(Menke et al.,2005). The technique has problems that are possibly related to various aspects,such as abyssal ooze,the height of the sphere shelf feet from the seabed,and the ocean bottom currents,which result in the observed P-wave amplitude and S-wave receiving effectiveness,and signal-to-noise ratios are different from those observed on the continental bedrock. Therefore,it is necessary to study the OBS receiver functions inversion method intensively. In this paper,we present the lithospheric S wave velocity structures beneath three stations at the southwestern subbasin in the SCS using the OBS receiver functions method,based on the observed teleseismic data from Dec. 2010 to Mar. 2011 obtained using a Chinese domestically I-4C type broadband OBS successfully(3 stations are shown in Fig. 1).

Fig. 1

Fig. 1 Layout of the 3D seismic array of OBSs at the southwestern subbasin in the SCS and the distribution of seismic events

The broadband OBS of I-4C type is developed by the Institute of Geology and Geophysics,Chinese Academy of Sciences(Yao and You,2011). The design of the OBS meets the international development trend and adopts some new techniques,such as a rechargeable lithium battery,blue tooth for data input or output,USB for data reading and GPS communication method,to follow the technology development in the world.The main part of the OBS includes a status controllable geophone of three components with a band of 60 s~50Hz(equipped in a container),a hydrophone and a compass for the azimuthal measurement of the horizontal component. On the input interface of the preamplifier circuit of a recorder,a low-pass anti-alias LC filter of the first order without a source is added,and the amplifier circuit of the equipment consists of very low noise,high accuracy and double calculation amplifier for high ability of anti-jamming. The A/D converter adopts ∑−∆ fourth-order delta modulator with a dynamic range > 120 dB,a data capacity 16 G,a total power consumption < 0.3 W and an inner clock accuracy better than 5×10⁻⁸ s. The data sampling interval is 8 ms. The posture controlling system uses a stepper motor to adjust the positions of the geophones to be horizontal and uses a decelerating motor to press geophones to contact with capsule closely to ensure a good coupling between them. The capsule ball is fixed on an iron heaving-coupling holder by four cables to maintain good coupling between them,even when the capsule ball is inclined. The iron heaving-coupling holder consists of five feet for inserting into the seabed to ensure good coupling between the equipment and the seabed,especially for horizontal coupling.

The passive seismic records used in this study were observed by the OBS stations deployed at the abandoned spreading ridge of the southwest subbasin with a water depth of ~3580 m. The southwestern subbasin in the SCS is a small sea basin formed by seafloor spreading and appears as a trumpet northeast forward opening. The ridge axis area is covered by the Zhongnan-Changlong seamount chain and has a lift valley with thick sediments, showing slow spreading characteristics. The magnetic anomaly indicates that the seafloor spreading time was 23~16 Ma(Briais et al.,1993; Sun et al.,2009),and the seamount distribution indicates some magmatism occurred after the cessation of seafloor spreading. The chronology analysis of basalt rock dredged from the seamount indicates that the age of rock was 14~3.5 Ma(Yan et al.,2006; Wang et al.,2009)and magmatism was not simultaneously stopped after the cessation of seafloor spreading. Signal analysis shows that,according to the frequency characteristics,submarine waves observed by OBS can be divided into two types: one type is a long wave with wavelength ~50 s,and the other is a short wave with higher frequency,which is possibly due to the height of iron heaving-coupling holder from the seabed. This short wave is prominently observed by the hydrophone component(Liu et al.,2012).

The receiver function method,which is based on the equivalent source hypothesis(Langston,1977),introduces a water level and Gauss low-pass filter function in the frequency domain to extract the so-called receiver functions from the observed data(Langston,1979). This method can calculate a 1D deep crustal structure effectively and has become a widely used technique in studies of Moho depth,V_p/V_s ratio,asthenosphere structure,lithospheric anisotropy and the interior of the Earth(Zhu and Kanamori,2000; Liu and Niu,2012; Tian,2008; Xu et al.,2011). This method has also been modified to invert the fine crustal structure based on islands of broadband seismometer teleseismic data(Leahy and Park,2005). Although some proposals have been reported(Menke et al.,2005),the application of the receiver function method in OBS data is not reported in the literature. We use the same scheme adopted by Zhu and Kanamori(2000)to extract the receiver functions from the OBS observed data and use the inversion approach adopted by Sambridge 1999a,Sambridge 2001b,Sambridge 2001),Sambridge and Mosegaard 2002to simulate the S wave velocity structure beneath the seafloor.

Seismic data often contain various types of noise: some types of noise present a stable long period feature with strong amplitude or a short period feature of high frequency,whereas other types consist of random nonstationary signals; this variability introduces difficulties in discriminating the real seismic phases. Typically,the seismic signals are assumed to be stationary signals and are processed by removing the mean value,removing the trending variation,and performing signal tapering and high/low pass Fourier transform filtering. Nevertheless,we discovered that the OBS data signal-to-noise ratio was not improved after processing using the typical approach. Therefore,we use the Fourier transform first to perform high/low pass filtering and to perform bandpass filtering to suppress the stationary disturbances,and then use the time-frequency features of the wavelet transform to conduct threshold filtering in the wavelet domain to suppress the non-stationary noise. For the signal(or function)f(t)2L2(R),its continuous wavelet transform(CWT)is defined as

(1)

where a is a scaling factor that reflects the frequency change of the wavelet; b is a time factor that reflects the position of wavelet in the time window. The CWT of signal f(t)with ψ(t)is the transvection between f(t)and wavelet ψ(t). The CWT is a quantitative representation of similarity degree between the signal and wavelet series. The inverse transform of CWT is

(2)

where which is the original function admissible conditions of WT_f(a,b). The discrete wavelet transform(DWT)of function f(t)∈L²(R)is

(3)

The inverse DWT is

(4)

The wavelet transform of all f(t)must satisfy the following condition:

(5)

The discrete function series { m,n; m,n 2 Z} is the so-called “frame" in mathematics.

To explain the difference between one-dimensional Fourier transform filtering and wavelet transform filtering,we use the soft threshold method(Donoho,1995)on the collected data as an example. Fig. 2a shows three-components of the ZRT waveforms,which were cut from an earthquake event without filtering and lasted 150 s. From Fig. 2a,the radial component(red waveform)and the tangential component(green waveform)have strong low frequency features. Particularly in the tangential component,the effective wave is masked by noise.In addition,all the components are masked by many types of high frequency noise. Fig. 2b shows the signal after one-dimensional Fourier transform filtering with waveform analysis. To obtain a valid signal,we attempted to combine the following filters: high or low pass filter,bandpass filter,and notch filter. However,the result was not satisfactory. Fig. 2b shows that the noise is suppressed but the effective wave is lost at the same time.Further,the wave polarity is destroyed at some regions. Furthermore,we first used the Fourier transform to suppress the stationary noise,and then used the one-dimensional wavelet adaptive threshold filtering method to suppress the non-stationary noise. The result shown in Fig. 2c indicates that the latter approach is better at suppressing noise and maintaining the waveform. As a result,the signal-to-noise ratio is improved drastically.

Fig. 2

Fig. 2 The comparison between Fourier transform filtering and wavelet transform filtering

The steps of the extraction of the receiver functions from the actual records are as follows. First,three components of the seismic records are rotated,based on the earthquake azimuth relative to the station(as determined by the OBS compass),to obtain the radial component,the tangential component,and the vertical component. Next,the time domain signals are transformed into the frequency domain,and then the vertical component is used to divide two horizontal components for two receiver functions obtaining(Langston,1979):

(6)

where X^r(ω),X^t(ω)and X^v(ω)are the radial,tangential and vertical component spectra,respectively,and R^r(ω)and R^t(ω)are the radial and tangential receiver functions spectra,respectively.

To ensure the stability of the receiver function calculation,a water level is introduced(Zhu and Kanamori,2000):

(7)

where R(ω)is the horizontal component spectrum,S(ω)is the equivalent source spectrum represented by vertical component spectrum,S^*(ω)is the complex conjugate of S(ω), is a Gaussian low-pass filter,and Cσ_O² is the water level. The factor(1+c)is used to compensate for the amplitude loss caused by the water level.

Finally,we obtained the receiver functions from the OBS observed data,as shown in Fig. 3. From these receiver functions,the P-s converted phases are easily recognized as the waves appearing several seconds after the P wave arrivals,indicating the relatively thin crust of the oceanic basin. Although some individual receiver functions extracted from various events at the same station have pronounced diversities that are not used in structure inversion,most receiver functions are quite consistent. These diversities possibly result from various aspects of the differences in the waves,such as focal depth,wave traveling path and event magnitude,and the seismic waves arriving at one station carry complicated random noise and path homogeneities. These problems remain to be further studied in the future. Because the individual receive functions have substantial differences,a single receiver function and stacked receiver functions are used in the inversion processing. The best fitted receiver function and the corresponding inverted model are selected as the final results.

Fig. 3

Fig. 3 The teleseismic receiver functions extracted from station observations

To overcome some of the issues of linear inversion methods,various nonlinear inversion techniques have been widely used,such as genetic algorithm(GA)(Shibutani et al.,1996),simulated annealing(SA)(Ingber et al.,1989),wavelet transformation(Wu et al.,2003),and the complex receiver function spectrum ratio(Liu and Li,1996). These methods are based on solving the global optimization problem by searching for the most optimized model from various models. We adopt the nonlinear inversion method known as the neighborhood algorithm(NA),which was proposed by Sambridge et al.(Sambridge,1999a; Sambridge,1999b; Sambridge,2001; Sambridge and Mosegaard,2002). This algorithm first preferentially samples a collection of models from a better fitting data domain and then determines the most optimized model in this collection. This procedure can improve the inversion precision compared with ordinary random searching to obtain a single optimal model. Correspondingly,the inversion code used in this study is NA-Sampler,in which the forward problem of the differential seismogram is based on the generalized reflection transmission matrix theory proposed by Kennett(1979). Furthermore,we use the H-k stacking method(Zhu and Kanamori,2000)to calculate the Moho depth and the V_p/V_s ratio as a test of the inverted model.

An active and passive synchronous seismic experiment was conducted using a 3D array of OBSs from Dec.2010 to Mar. 2011 at the abandoned spreading ridge of the southwestern subbasin in the SCS(Fig. 1)(Ruan et al.,2012). The passive data includes six earthquakes of magnitude ≥ M_s7.0 recorded by three Chinese made I-4C type broadband OBSs. The selected data have high signal-to-noise ratios and epicenter distances of 30° to 63°. The OBS station coordinates and the earthquake catalog are presented in Table 1 and Table 2,respectively.

Table 1

Table 1 Coordinates and recording times of the OBS stations

Table 1 Coordinates and recording times of the OBS stations

Table 2

Table 2 Catalog of earthquakes (Ms ≥ 7.0) selected

Table 2 Catalog of earthquakes (Ms ≥ 7.0) selected

Using the above-described data processing,receiver function calculation and inversion methods,we obtained the lithospheric 1D S-wave velocity structure(beneath the seafloor)of three OBS stations at the southwestern subbasin in the SCS,as shown in Fig. 4. The results show that the uppermost crust is a 1~2 km thick low velocity layer,being especially thin at the OBS06 station,which is close to the ridge center axis. We propose that the shallow low velocity layer comprises sediments mixed with clastic rock from magma eruption that occurred after the ceasing of sea floor spreading. However,this layer is difficult to subdivide. Previous studies also indicated that,along both sides of the Changlong seamounts,the sediment thickness varies from 0.5 to 1 km at the central southwest subbasin and gradually increases northwestward and southeastward to 2~2.5km(Yao,1998). Seismic reflection profiles revealed that the spreading ridge has thick sediments(Gao et al.,2009). This prior observation is consistent with our results. The three 1D S-wave velocity models all indicate the presence of an abrupt velocity increasing discontinuity at a depth of 11~12 km,which is suggested to be Moho discontinuity. However,the active seismic 2D P-wave velocity profile indicated that the Moho depth of the same area is at approximately 10~12 km(below sea level)(Zhang,2013). Another 2D P-wave velocity profile also indicated that the Moho depth near our research area is at approximately 10 km(Pichot et al.,2014). Obviously,the data show that S-wave velocity model is different from the P-wave velocity model,the latter being shallower by 3~4 km. This is an important scientific discovery. We interpreted this result in two aspects. First,the Moho prefers to be a thick layer rather than a simple layer(Coleman,1977). Second,the S-wave velocity is more sensitive to melting and high temperature conditions than the P-wave velocity.

Note that the 1D S-wave velocity model beneath station OBS06 has a low velocity layer at depths of 6~12 km(above the Moho),different from the usual velocity downward increase,as observed at the other two stations. This difference may be explained by station OBS06 being closest to the spreading ridge axis(rift valley)where magma upwelled. Therefore,we suggest that this low velocity layer near the axis consists of melting magma or chamber. The numerical simulation of the thermodynamic structure showed that at depths of 8~12 km,the strain rate has an obvious change; in addition,at that depth,the oceanic crust is relatively weak and undergoes large tectonic deformation (Zhang and Shi,2008). Based on a 2D P-wave velocity model, Zhang(2013) suggested that there exist the combined effects of a detachment fault at the shallow crust and mantle melting at the lower crust.

Another remarkable velocity anomaly in the 1D models of OBS06 and OBS34 occurs at the depths from 17 km to 30 km,where the mantle has a small negative S-wave velocity gradient,whereas at the OBS35 station,at the depth of 26 km,the S-wave velocity decreases sharply. We speculate that a hot material supply source is present at the deep mantle beneath the spreading center axis. Correspondingly,previous studies have shown the following: 1)the SCS magmatism mainly occurred after the ceasing of seafloor spreading,and 2)the seamount age is young(14~3.5 Ma)(Yan et al.,2006; Wang et al.,2009). The heat flow measurements also showed that the southwest subbasin has a high heat flow value of 100~150 mW·m⁻²(Chen and Lin,1997;Zhang et al.,2001). As a result,we propose that,at the current tectonic and stress environment,although seafloor spreading has stopped,upper mantle melting still exists. This hypothesis merits future study because it has a great significance for our understanding of the deep hydrothermal state of the southwest subbasin in the SCS.

Fig. 4

Fig. 4 Density diagram of the best 1000 S velocity models generated by the NA code

We used the evaluation method of NA-Sampler to evaluate the above 1D S-wave velocity models. The Gauss factor is equal to 5,and the water level is equal to 0.001 in the forward calculation. The result is shown in Fig. 5. From Fig. 5,the receiver functions of all three models(blue solid line)are found to be very close to detected receiver functions(red solid line),especially within a few seconds after P wave. Pms and other Moho multiple wave are clear in the figure. The latter multiple waves are also clear. Thus,our result is reliable.

Fig. 5

Fig. 5 Comparison between the observed receiver function and the predicted receiver function

On the Moho discontinuity,many multiple P-s converted waves exist,for which the arrivals are very sensitive to the ratio variation of the Moho depth versus V_p/V_s. Thus,we can use the relationship between the travel times of multiple waves and the Moho depth along with the H-k stack procedure(Zhu and Kanamori,2000)to test the Moho depth inverted(Fig. 6). Within the depth range from 7 to 25 km,the V_p/V_s ratio ranges from 1.5 to 2 and the weighting factors are 0.7,0.2 and 0.1(Huang et al.,2011),the testing results of the grid searching stack shows that the Moho depths beneath OBS06,OBS34 and OBS35 are 12.5 km,10 km and 12.5 km,respectively,and the V_p/V_s ratios are 1.85,1.69 and 1.77,respectively. These results indicate that the obtained 1D velocity models have a good stability. Consequently,we obtained the crust mean Poisson ratios of 0.2936,0.2306 and 0.2656 for the three models of OBS06,OBS34,and OBS35,respectively,by using the following equation:

Fig. 6

Fig. 6 The test results of the H-k stack procedure

(8)

where σ is the crustal Poisson ratio,and k is the V_p/V_s ratio.

We combined the use of the Fourier transform and the use of the wavelet transform to suppress the non-stationary disturbances in the OBS data successfully to acquire clear seismic phase records with high signal-to-noise ratio,and then conducted an analysis of the lithospheric structure of the southwest subbasin ridge in the SCS by using the receiver functions inversion technique for teleseismic data. The main conclusions of this study are as follows:

(1) The receiver functions inversion technique is feasible for use in analyzing OBS passive data,and suppressing the non-stationary disturbances is a key in data processing.

(2) The Moho depth of the southwest subbasin determined by S wave is 10~12 km(from the seabed)with a crust thickness of 6~8 km. This depth is 3 km deeper than the result determined by active seismic P-wave analysis. The shallow crust has a low velocity zone with a sediment thickness of 1~2 km; this crust possibly consists of sediments mixed with rock debris and fracture that resulted from magma eruption after sea floor spreading ceased.

(3) The models show that there is a low S-wave velocity zone at the depth of 6~12 km,just above the Moho discontinuity. We propose that the spreading center has a melting lower crust or magma chamber. At 17~30 km,the S-wave velocity is found to have a negative gradient. We speculate that a hot material supply source is present at the deep mantle beneath the spreading center axis. Although the seafloor spreading has stopped in the current tectonic and stress environment,upper mantle melting can still exist.

We thank all colleagues who joined this OBS experimental investigation of southwestern subbasin in the SCS in 2010–2011: Wu Z L,Dong C Z,Niu X W and Wei X D from the Second Institute of Oceanography,SOA,Zhu J J and Huang H B from the South China Sea Institute of Oceanology,Chinese Academy of Sciences,and Liang J W,Wang L J,Qiu M C,and Xu W F from Taiwan Ocean University. The investigation was conducted by separately using at different times in the following research vehicles: R/V “Fendou No.7",R/V “Kan 407",and R/V “hiyan No. 2". We used GMT software(Wessel and Smith,1995) for plotting some of the figures,the code provided by Lupei Zhu(Zhu and Kanamori,2000) for extraction of the receiver functions,and the code of NA-sampler(Sambridge et al.,2001,2002)for receiver functions inversion.