# Receiver function HV ratio: a new measurement for reducing non-uniqueness of receiver function waveform inversion

Receiver function HV ratio: a new measurement for reducing non-uniqueness of receiver function... Summary It is known that a receiver function has relatively weak constraint on absolute seismic wave velocity, and that joint inversion of the receiver function with surface wave dispersion has been widely applied to reduce the trade-off of velocity with interface depth. However, some studies indicate that the receiver function itself is capable for determining the absolute shear-wave velocity. In this study, we propose to measure the receiver function HV ratio which takes advantage of the amplitude information of the receiver function to constrain the shear-wave velocity. Numerical analysis indicates that the receiver function HV ratio is sensitive to the average shear-wave velocity in the depth range it samples, and can help to reduce the non-uniqueness of receiver function waveform inversion. A joint inversion scheme has been developed, and both synthetic tests and real data application proved the feasibility of the joint inversion. Joint inversion, Body waves, Crustal structure 1 INTRODUCTION Receiver function analysis is a widely used technique for resolving crustal and upper-mantle structures. It is sensitive to the local structure beneath a seismic station, and is almost insensitive to the earthquake source and propagation effects through the mantle (Burdick & Langston 1977; Lang 1979). It has been applied to obtain the crustal structure beneath a seismic station by waveform inversion (Owens et al. 1984; Randall 1989; Kind et al. 1995; Zhao et al. 1996; Sambridge 1999; Chang et al. 2004; Lawrence & Wiens 2004; Vinik et al. 2004) or image the seismic interfaces by converted phases or multiples found within the receiver functions (Zhu & Kanamori 2000; Chen et al. 2005; Rondenay 2009; Cheng et al. 2016). It has been generally accepted that a receiver function has a non-uniqueness issue in constraining the absolute velocity because of the strong trader-off between the velocity and interface depth (Ammon et al. 1990). A popular way to release the non-uniqueness is through joint inversion of the receiver function waveform (RFWF) and surface wave dispersion (e.g. Last et al. 1997; Özalaybey et al. 1997; Du & Foulger 1999; Juliá et al. 2000; Chang et al. 2004; Lawrence & Wiens 2004; Tkalčić et al. 2006; Moorkamp et al. 2010; Bodin et al. 2012; Shen et al. 2013; Chen & Niu 2016; Kim et al. 2016b). Recently, Chong et al. (2016) tried to invert receiver function with Rayleigh wave ellipticity, and found that the joint inversion had better constraints on the velocity gradient than the absolute velocity. Some studies, however, indicate that the receiver function itself is capable of resolving the absolute shear-wave velocity (Ammon 1991; Svenningsen & Jacobsen 2007) and the non-uniqueness of the waveform inversion can be reduced by using other aspects of receiver function (Wu et al. 2007; Jacobsen & Svenningsen 2008; Li et al. 2017). Wu et al. (2007) and Li et al. (2017) used multifrequency receiver functions to reduce the non-uniqueness of RFWF inversion. Jacobsen & Svenningsen (2008) stratified the velocity model in quasi-continuous delay time to reduce the non-uniqueness. Svenningsen & Jacobsen (2007) used the apparent shear-wave velocity (Vs,app) measured from RFWF to constrain shear-wave velocity structure. Inspired by this Vs,app approach, Kieling et al. (2011) and Peng et al. (2012) developed empirical relations between the Vs,app and the shear-wave velocity structure to construct starting models for the kernel based RFWF inversion, and they claimed that the non-uniqueness was reduced. Schiffer et al. (2015, 2016) applied the joint inversion of the RFWF and Vs,app to study crustal and upper-mantle structure. Hannemann et al. (2016) studied oceanic lithospheric shear-wave velocity structure with Vs,app measured from the P wave recorded by an OBS (ocean bottom seismometer). All these studies show that P wave particle motion is an important parameter for determining crustal structure, and it can help to improve waveform inversion from receiver functions. Although the Vs,app method provides constraints on absolute Vs, it seems to have weak sensitivity to deep discontinuities. Instead, the time-domain RFWF may resolve discontinuities better. A joint inversion with Vs,app and receiver function might improve inversion of both absolute Vs and interface positions. In this study, we propose to measure the time and frequency dependent amplitude information of receiver function to determine the crustal shear-wave velocity. Since the proposed measurement contains the amplitude ratio between the radial and vertical P wave train, and we name it receiver function HV ratio. The receiver function HV ratio contains not only the information of P wave amplitude ratio (particle motion) but also the amplitude information of other phases such as P-to-S conversion and its multiples. We investigated the sensitivity of the receiver function HV ratio to the velocity structure, and the feasibility of joint inversion of the HV ratio and waveform of the receiver function. Both synthetic tests and real data applications indicate that the receiver function HV ratio can be used to reduce the non-uniqueness of RFWF inversion. 2 RECEIVER FUNCTION HV RATIO (RFHV) Studies indicate that P wave particle motion carries velocity structure information beneath a seismic station (Kruger 1994; Svenningsen & Jacobsen 2007; Ni et al. 2014). For a half space model, P wave particle motion on the surface is determined by the ray parameter and the subsurface shear-wave velocity (Svenningsen & Jacobsen 2007), which can be derived from the relation between the apparent and true incident angle of P wave (Nuttli & Whitmore 1961). The relation can be written as following:   $$\mathop V\nolimits_S = \frac{{\sin [{\textstyle{1 \over 2}}\overline {{i_p}} ]}}{p},$$ (1)where $${\bar{i}_p}$$is the apparent incident angle of P wave, and p is the ray parameter. For the teleseismic P wave, the ray parameter can be calculated based on a reference earth model, thus the apparent incident angle or amplitude ratio can be used to determine the subsurface shear-wave velocity. Direct measuring of the apparent incident angle could be difficult due to the complicated seismic P waveform (Kruger 1994). Svenningsen & Jacobsen (2007) proposed a convenient way to measure the apparent incident angle from receiver function. They measure the frequency dependent apparent incident angle of P wave, and derive the apparent shear-wave velocity Vs, app(T) to determine the crustal velocity structure. Here, we propose to measure a transformation of the RFWF, which contains the amplitude ratio between the radial and vertical P wave train. The amplitude ratio is measured for different time window, with a squared-cosine smooth filter (weighting function) centred at t = 0 s following Svenningsen & Jacobsen (2007). It is named the ‘receiver function HV ratio’ (shortened as RFHV hereafter), which can be measured as following:   $${\rm{RFHV(}}T{\rm{)}} = {{{\int_{{{\rm{ - }}T}}^{T}{{{R_{RF}}(\tau ){{\cos }^2}\left({{{\pi \tau } \over {2T}}}\right){\rm{d}}\tau }}} \over {\int_{{ - T}}^{T}{{{Z_{RF}}(\tau ){{\cos }^2}\left({{{\pi \tau } \over {2T}}}\right){\rm{d}}\tau }}}}},$$ (2)where T is half length of the time window, RRF is the radial RFWF. And ZRF is ‘Z’ receiver function (Svenningsen & Jacobsen 2007) obtained from deconvolution of vertical P wave train by itself, which is ideally the Gaussian filtering function used in calculating receiver function. In this equation, the radial receiver function is integrated in a variable time window, with a time window depended squared-cosine weighting function which smooth the waveform for integration. It is then scaled with a scaling factor which is calculated from the ZRF with a integration in the same way. Since in the calculation of receiver function, a Gaussian filtering function is applied, thus the scaling factor helps to remove the effect of the filtering function. The RFHV contains information on P wave particle motion and also amplitude information for P-to-S converted phases or multiples from Moho and other interfaces. Thus we can invert for the velocity structure beneath the seismic station by fitting the RFHV ratio curve, without constructing a starting velocity model for inversion using empirical relations (Kieling et al. 2011; Peng et al. 2012). Before applying RFHV ratio to invert for the crustal structure, we take a one layered crustal model (Fig. 1a) as an example to analyse the characteristics and sensitivity of this new measurement to the velocity structure. Synthetic receiver functions are calculated from the ray parameter range from 0.04 to 0.08 s km–1 with an increment of 0.005 s km–1. The Gaussian filtering factor is chosen to be 2.5, which is a commonly used value for crustal structure study. Here only radial receiver functions are displayed (Fig. 1b). The RFHV ratio curves are then calculated following eq. (2), and are shown in Fig. 1(c). We found that for this synthetic example, the effects of the Gaussian factor on the RFHV ratio is fairly small in the commonly applied range (1.0–4.0) due to the scaling factor in the eq. (2). As we can see in Fig. 1(c), for each ray parameter, the RFHV ratio curve varies mostly in the time window of T < 20 s and varies slightly after this time window. This is because most of the multiples generated in the crust end before 20 s. For the first few seconds, the measurement is related to P wave HV ratio, and each of amplitude jumps on the RFHV ratio curves correspond to the contribution of a new phase coming into the integral time window. Another important feature is that the RFHV ratio increases almost linearly with the ray parameter for a given T (integral parameter), as shown in Fig. 2(a). Here, the RFHV ratios with integral parameter T of 1.0, 6.0, 10.0 and 20.0 s are taken as an example to demonstrate their linear relationship with the ray parameter. Because of this relationship, the real data measurements can be fitted with a line in the RFHV and ray parameter domain for each filter-parameter to depress the anomalies, and RFHV ratio for given ray parameters can be obtained for further analysis. Since the RFHV maybe unevenly sampled in the ray parameter domain, linear fitting of the measurements in RFHV and ray parameter domain helps to avoid the bias, thus reduce the uncertainty of the observations. For sensitivity analysis, we calculate the sensitivity kernel numerically with finite difference method. We find that the RFHV ratio is more sensitive to Vs, and less sensitive to density and Vp/Vs, while being almost insensitive to Vp, as shown in Fig. 2(b). The longer time window measurement is more sensitive to the deeper structure (Fig. 2c). Also, the sensitivity kernel for Vs is mostly positive, which means the RFHV ratio is sensitive to the average shear-wave velocity in the depth range it samples. Figure 1. View largeDownload slide One layered crustal model, synthetic radial RFWFs and the corresponding RFHV ratio curves. (a) The one layered crustal model, with the Moho located at a depth of 30 km, the Vp/Vs ratio is 1.732. (b) The synthetic P wave RFWFs (radial component) with the Gaussian factor of 2.5, displayed according to the ray parameters. (c) The synthetic RFHV ratio curves with ray parameters marked. Figure 1. View largeDownload slide One layered crustal model, synthetic radial RFWFs and the corresponding RFHV ratio curves. (a) The one layered crustal model, with the Moho located at a depth of 30 km, the Vp/Vs ratio is 1.732. (b) The synthetic P wave RFWFs (radial component) with the Gaussian factor of 2.5, displayed according to the ray parameters. (c) The synthetic RFHV ratio curves with ray parameters marked. Figure 2. View largeDownload slide Features of RFHV ratio. (a) Almost linear relation between the RFHV ratio and the ray parameter for each integral parameter T. Examples of four different integral parameters (T = 1.0, 6.0, 10.0 and 20.0 s) are plotted with different symbols. (b) The sensitivity kernels respected to Vp, Vs, density and Vp/Vs ratio for T = 6.0 s. (c) Variation of the sensitivity kernel (Ks) with integral parameter T. Figure 2. View largeDownload slide Features of RFHV ratio. (a) Almost linear relation between the RFHV ratio and the ray parameter for each integral parameter T. Examples of four different integral parameters (T = 1.0, 6.0, 10.0 and 20.0 s) are plotted with different symbols. (b) The sensitivity kernels respected to Vp, Vs, density and Vp/Vs ratio for T = 6.0 s. (c) Variation of the sensitivity kernel (Ks) with integral parameter T. From the analysis, we see that the RFHV ratio is sensitive to the average shear-wave velocity, thus it may be combined with the RFWF to reduce the non-uniqueness of the inversion. 3 JOINT INVERSION OF RFHV AND RFWF As demonstrated in former section, the RFHV ratio is sensitive to the average shear-wave velocity, which could help to reduce the non-uniqueness of RFWF inversion. Here we developed a tool to invert for crustal velocity structure with these two measurements and investigate the feasibility of the joint inversion method. We will describe the inversion method, and present numerical tests and real data application to demonstrate the feasibility of the joint inversion. 3.1 Inversion method The crustal model is parameterized into sublayers with shear-wave velocity (Vs), ratio of P-wave velocity and shear-wave velocity (Vp/Vs), and layer thickness (hi). The density of each layer is scaled to Vp using the relation by Brocher (2005). A subroutine by Shibutani et al. (1996) is modified to calculate the radial and ‘Z’ receiver functions, and the RFHV ratio is calculated following eq. (2) in Section 2. The objective function for fitting both the RFWF and RFHV ratio is defined as following:   \begin{eqnarray} ObjFun &=& {{{{W_{{\rm{RFWF}}}}} \over {{N_{{\rm{RFWF}}}}}}}\sum\limits_{i = 1}^{{N_{{\rm{RFWF}}}}}\nonumber\\ &&{\sqrt {{{1 \over {NPT{S_i}}}}\sum\limits_{j = 1}^{NPT{S_i}} {{{{{{\left( {RFWF_{\rm obs}^i(j) - RFWF_{\rm syn}^i(j)} \right)}^2}} \over {\sigma _{{\rm{RFWF}}}^2(j)}}}} } }\nonumber\\ && + \,{{{{W_{{\rm{RFHV}}}}} \over {{N_{{\rm{RFHV}}}}}}}\sum\limits_{k = 1}^{{N_{{\rm{RFHV}}}}}\nonumber\\ &&{\sqrt {{{1 \over {N{T_k}}}}\sum\limits_{l = 1}^{N{T_k}} {{{{{{\left( {RFHV_{\rm obs}^k(l) - RFHV_{\rm syn}^k(l)} \right)}^2}} \over {\sigma _{{\rm{RFHV}}}^2(l)}}}} } } ,\nonumber\\ && +\, {W_S}*Smoothness \end{eqnarray} (3)where WRFWF and WRFHV are the weights for RFWF and RFHV, respectively, and WS is the weight for model smoothness factor. RFWFobs and RFWFsyn represent observed and synthetic radial receiver functions, while RFHVobs and RFHVsyn denote observed and synthetic receiver function HV ratios. NPTS and NT are the number of points in the RFWF and RFHV ratio. Uncertainties of the RFWF and the RFHV ratio (σRFWF and σRFHV) are used to weight each measurement. To avoid abnormal velocity jumps in the inversion, a smoothness factor is applied to the objective function, which is defined as following,   $$Smoothness = \sqrt {{{1 \over {{N_{\rm layer}}{\rm{ - 1}}}}}\sum\limits_{n = 1}^{{N_{\rm layer}}{\rm{ - 1}}} {{{{{{\left( {{V_j} - {V_{j + 1}}} \right)}^{\rm{2}}}} \over {V_{\rm mean}^2}}}} }$$ (4)Or   $$Smoothness {=} \sqrt {{{1 \over {{N_{\rm layer}}{\rm{ - 2}}}}}\sum\limits_{n = 1}^{{N_{\rm layer}}{\rm{ - 2}}} {{{{{{\left( {{{1 \over 2}}\left( {{V_j} {+} {V_{j {+} {\rm{2}}}} - {V_{j + {\rm{1}}}}} \right)} \right)}^2}} \over {V_{\rm mean}^2}}}} } ,$$ (5)where Nlayer is the number of layers in the model, and Vj is the velocity of the jth layer. Eqs (4) and (5) are for the first and second order velocity changes, respectively, with the former favouring half-space like models, while the latter favours velocity models with gradients (Chong et al. 2016). A fast simulated annealing algorithm is employed as the driver for the inversion. It is a non-linear inversion method, which searches for the optimal solution in a given model space. Just like other Monte Carlo non-linear inversion methods, it only needs forward modelling of each measurement without calculating the sensitivity kernels (Rothman 1986; Ingber 1989). More description of the fast simulated annealing algorithm can be found in section 4 of Chong et al. (2015). 3.2 Synthetic tests In order to test the feasibility of the joint inversion method, we carried out numerical inversion tests before applying it to real data. Similar to the inversions in fig. 4 of Chong et al. (2016), we use a velocity model composed of smoothly varying velocity layers under a low velocity layer for the synthetic inversion tests. Since RFWF may fail to determine the smoothly varying velocity structure (Chong et al. 2016), this kind of model will be a good choice for investigating the stability of the inversion. RFWF and RFHV are calculated with a Gaussian factor of 2.5 and ray parameter of 0.06 s km–1. The uncertainty of the radial receiver function is set to be 5 per cent of its maximum amplitude, and for the RFHV ratio it is set to be 5 per cent of its mean value. Numerical tests are carried out for inversions with only single data set and combination of the two data sets. During the inversions, the Vp/Vs ratio is fixed at 1.732, the same with the value of the true model, in which way only the shear-wave velocity and layer thickness are inverted. As one may expect, the inversion with RFWF alone is able to resolve features of the sharp velocity jumps, however, it cannot determine the absolute shear-wave velocity (Vs), as in Figs 3(a) and (b). Although RFWF is well fitted, the optimal velocity model from the inversion cannot predict the RFHV ratio curve (Fig. 3c). From the comparison, the optimal model predicted RFHV ratio is larger than the synthetic data, which indicates that the inverted optimal model is faster in shear-wave velocity than the true model. On the contrary, the inversion with RFHV ratio alone is able to determine a model which matches the average velocity of the true model, however, the sharp velocity jumps like the thin low velocity layer and the Moho interface are not captured (Figs 3d and e). Consequently, the optimal model inverted with RFHV ratio alone cannot predict RFWF (Fig. 3f). The single data set inversion tests suggest that there is complementarity information between the RFHV and RFWF, and that the combination should do a better job. Indeed, the joint inversion of the two data sets is able to resolve a better crustal model, and both the velocity jumps and the absolute velocity are well determined, as in Figs 3(g)–(i). Figure 3. View largeDownload slide Numerical inversion tests for a smooth variation model with a low velocity layer. (a–c) Inversion with receiver function only. The input model (blue) and optimal model (red) are displayed in (a), while the grey lines indicate the 100 best models from the inversion. (b) Waveform comparison between the data (thick black curve) and the optimal model prediction (red curve). Thin black curves indicate the uncertainty of the data. (c) Comparison between the RFHV ratio data (black) and the optimal model prediction (red). (d–f) Inversion of RFHV ratio only. (d) Same as (a) but for the result of RFHV ratio inversion. (e) Comparison between the data (thick black curve) and the optimal model predicted RFHV ratio (red curve), uncertainty is shown with the black bars. (f) Comparison between the RFWF data (black) and the optimal model prediction (red). (g–i) Joint inversion of the two data sets. (g) Same as (a) but for joint version. (h) Waveform fitting of the RFWF. (i) Fitting of the RFHV ratio curve. Notice that the data and the optimal model prediction overlap with each other. Figure 3. View largeDownload slide Numerical inversion tests for a smooth variation model with a low velocity layer. (a–c) Inversion with receiver function only. The input model (blue) and optimal model (red) are displayed in (a), while the grey lines indicate the 100 best models from the inversion. (b) Waveform comparison between the data (thick black curve) and the optimal model prediction (red curve). Thin black curves indicate the uncertainty of the data. (c) Comparison between the RFHV ratio data (black) and the optimal model prediction (red). (d–f) Inversion of RFHV ratio only. (d) Same as (a) but for the result of RFHV ratio inversion. (e) Comparison between the data (thick black curve) and the optimal model predicted RFHV ratio (red curve), uncertainty is shown with the black bars. (f) Comparison between the RFWF data (black) and the optimal model prediction (red). (g–i) Joint inversion of the two data sets. (g) Same as (a) but for joint version. (h) Waveform fitting of the RFWF. (i) Fitting of the RFHV ratio curve. Notice that the data and the optimal model prediction overlap with each other. 3.3 Application to the Indian craton In this section, we apply this method to station HYB in the Indian craton to obtain the crustal structure. Since the crustal structure beneath HYB is relatively simple and has been extensively studied (Zhou et al. 2000; Sarkar et al. 2003; Kiselev et al. 2008; Julia et al. 2009; Chong et al. 2016; Li et al. 2017), it is a suitable place to investigate the feasibility of the method. As done for the synthetic tests, single data set inversions and two data set joint inversions are performed. 3.3.1 Data processing High quality teleseismic waveforms from 785 earthquakes that occurred between 2000 and 2007 with epicentral distances between 30° and 90° and magnitudes between M6 and M7, were collected for the station HYB (Fig. 4). The instrument responses, trends and means of each component were removed from the raw data. All the waveforms are then selected according to the signal-to-noise ratio of P waves on the three components (ENZ). After the waveform selection, P wave trains were obtained by applying a 100-s time window which starts from 10 s before the onset of P wave. The E–N components were then rotated to radial and tangential components. Receiver functions are calculated with an iterative deconvolution method (Kikuchi & Kanamori 1982; Ligorría & Ammon 1999). Figure 4. View largeDownload slide Station location of HYB in Indian Craton and distribution of the earthquakes. Black triangle for the station and circles for the earthquakes. Figure 4. View largeDownload slide Station location of HYB in Indian Craton and distribution of the earthquakes. Black triangle for the station and circles for the earthquakes. In total, 403 high quality RFWFs were selected for further processing, Fig. 5(a). We first remove the outliers according to the mean and standard deviation of the RFWF, this process helps to remove waveforms with long period noise. Then the waveforms are further selected according to the cross-correlation coefficients between the individual waveform and the mean RFWF. For the waveform inversion, we need to stack the radial receiver functions to reduce the number of records and improve the signal-to-noise ratio. We carried out two-step stacking of the radial receiver functions. For the first step, the raw radial receiver functions are stacked using a ray parameter bin of 0.002 s km–1 with an overlap of 0.002 s km–1. This procedure generates stacked receiver functions that are nearly evenly spaced in the ray parameter domain, avoiding large contributions from some ranges of the ray parameter. In the second step, the radial receiver functions are further processed with a ray parameter based stacking method proposed by Chen & Niu (2013), which corrects for the moveouts and amplitudes for all ray parameters based on a one layer crustal model obtained with H–K stacking method (Zhu & Kanamori 2000). And the method is modified to calculate the uncertainties for the stacked waveform, which will be used for the data weighting during the inversion. The radial receiver functions are stacked to three reference ray parameters (0.05, 0.06 and 0.07 s km–1) with the ray parameter based stacking method (Fig. 5b). Figure 5. View largeDownload slide RFWF and RFHV ratio curves for station HYB. (a) Radial RFWF measured from the records of earthquakes in Fig. 4. (b) Stacked radial RFWFs with ray parameter correction. Uncertainties of the stacked waveform are shown with thin lines. (c) RFHV ratio curves. (d) RFHV ratio curves for three ray parameters (0.05, 0.06 and 0.07 s km−1) calculated from the linear fitting of the measurements in (c), the uncertainties are shown with vertical bars. Details of the linear fitting can be found in Fig. 6. Figure 5. View largeDownload slide RFWF and RFHV ratio curves for station HYB. (a) Radial RFWF measured from the records of earthquakes in Fig. 4. (b) Stacked radial RFWFs with ray parameter correction. Uncertainties of the stacked waveform are shown with thin lines. (c) RFHV ratio curves. (d) RFHV ratio curves for three ray parameters (0.05, 0.06 and 0.07 s km−1) calculated from the linear fitting of the measurements in (c), the uncertainties are shown with vertical bars. Details of the linear fitting can be found in Fig. 6. RFHV ratios are computed from the raw radial receiver functions shown in Fig. 5(a) according to eq. (2) in Section 2. After eliminating the outliers, we obtained 252 RFHV ratio curves, as shown in Fig. 5(c). The RFHV ratios are then fitted in the ray parameter domain with a line for each integral parameter (T) according to the linear relation between the RFHV ratio and the ray parameter, as shown in Figs 6(a)–(c). And we used the RFHV values on the fitted line of three ray parameters (0.05, 0.06 and 0.07 s km−1) for inversion (Fig. 5d). Figure 6. View largeDownload slide (a–c) Linear fitting of the RFHV ratios measured for station HYB. Linear fitting of RFHV ratios with integral parameters of T = 6, 10 and 20 seconds are shown in (a–c), respectively. Open triangles for data, thick black lines for fitted lines of RFHV ratio of each integral parameter. Thin grey lines indicate the one standard deviation of linear fitting, and half of the standard deviation is used for the uncertainty of measurement. The linear fitting is carried out with the function ROBUSTFIT in MATLAB. It performs for a robust multilinear regression, and in the data processing we used default algorithm with iteratively reweighted least squares with a bi-square weighting function. Figure 6. View largeDownload slide (a–c) Linear fitting of the RFHV ratios measured for station HYB. Linear fitting of RFHV ratios with integral parameters of T = 6, 10 and 20 seconds are shown in (a–c), respectively. Open triangles for data, thick black lines for fitted lines of RFHV ratio of each integral parameter. Thin grey lines indicate the one standard deviation of linear fitting, and half of the standard deviation is used for the uncertainty of measurement. The linear fitting is carried out with the function ROBUSTFIT in MATLAB. It performs for a robust multilinear regression, and in the data processing we used default algorithm with iteratively reweighted least squares with a bi-square weighting function. 3.3.2 Inversion results Inversions with single data sets and the combination of the two data sets are carried out. From the inversion of only RFWF, we get an optimal model with the Moho located at the depth of about 29 km, and the uppermost mantle shear-wave velocity of about 4.0 km s–1 (Fig. 7a). However, the Moho depth and shear-wave velocity of the uppermost mantle of the optimal model are quite different from previous studies (Zhou et al. 2000; Sarkar et al. 2003; Kiselev et al. 2008; Julia et al. 2009; Li et al. 2017). Although the optimal model is able to fit the RFWF, the RFHV ratio curves are not well-predicted (Figs 7b and c). From the comparison, the amplitude of predicted RFHV ratio is smaller than the data, which means that the shear-wave velocity of the optimal model is slower than that of the real model. The inversion with only RFHV ratio curves gives us a model in which velocity increases smoothly with the depth (Figs 7d and e), but the Moho interface has not been captured. Thus the RFWF, especially the P-to-S conversion and multiples from the Moho, are not well predicted (Fig. 7f). The joint inversion of these two data sets, gives us a crustal model with a Moho located at depth of about 32 km, and the shear-wave velocity increasing smoothly from about 3.5 km s–1 near the surface to roughly 3.8 km s–1 above the Moho (Figs 7g–i). From the comparison with other studies (Sarkar et al. 2003; Kiselev et al. 2008; Zhou et al. 2000; Julia et al. 2009; Li et al. 2017), we found the crustal model from the joint inversion of RFWF and RFHV is similar to other models in the lower crust, but slightly slower in the upper and middle crust than most previous velocity models except for the one by Kiselev et al. (2008) (Fig. 8(a)). The shear-wave velocity of the uppermost mantle is about 4.3 km s–1, which is close to the model by Li et al. (2017), but a little bit slower than others. Figure 7. View largeDownload slide Inversions for the station HYB in Indian Craton. (a–c) Inversion with RFWF alone. (d–f) Inversion with RFHV ratio alone. (g–i) Jointly inversion of the two data sets. Figure 7. View largeDownload slide Inversions for the station HYB in Indian Craton. (a–c) Inversion with RFWF alone. (d–f) Inversion with RFHV ratio alone. (g–i) Jointly inversion of the two data sets. Figure 8. View largeDownload slide Figures for comparing the results of this study with previous studies. (a) Crustal model comparison of the station HYB in Indian Craton. (b) Crustal model comparison with the one by Chong et al. (2016) from the joint inversion of RFWF and Rayleigh wave ellipiticity (SWZH). (c) Comparison between the observed Rayleigh wave ellipiticity (black) and prediction (red) by the optimal model of this study. (d) Comparison between the observed RFHV ratio curves (black) and prediction (red) by the model of Chong et al. (2016). Figure 8. View largeDownload slide Figures for comparing the results of this study with previous studies. (a) Crustal model comparison of the station HYB in Indian Craton. (b) Crustal model comparison with the one by Chong et al. (2016) from the joint inversion of RFWF and Rayleigh wave ellipiticity (SWZH). (c) Comparison between the observed Rayleigh wave ellipiticity (black) and prediction (red) by the optimal model of this study. (d) Comparison between the observed RFHV ratio curves (black) and prediction (red) by the model of Chong et al. (2016). 4 DISCUSSION The inversion of RFWF usually has difficulties in constraining the absolute shear-wave velocity. The RFHV technique, which is based on a integral transformation of the RFWFs, provides a way to constrain the absolute shear-wave velocity. Sensitivity analysis and numerical inversions indicate that RFHV provides constraints on the average shear-wave velocity (Figs 2b–c and Figs 3d–e). Since RFWF provides constraints to the sharp interfaces, the RFHV and RFWF can be combined to put tighter constraints on velocity. The numerical inversion tests and real data application indicate that joint inversion is more reliable than inversions with single data set (Figs 3 and 7). Both the RFWF and RFHV ratio are single station measurements with localized sensitivity, and the joint inversion can be applied to regions with either a sparse network or dense array. In real data applications, the joint inversion may be affected by the background noise. We carried out two numerical tests to investigate the effect of noise on the joint inversion of RFHV and RFWF. The crustal model in Fig. 3 is used for the tests. Radial and vertical receiver functions are calculated with a Gaussian factor of 2.5 and ray parameter of 0.06 s km–1. White Gaussian noise is generated with the function AWGN in MATLAB (see data and resources) and added to the radial receiver function but not to the ‘Z’ receiver function, because the noise has been cancelled on the vertical component during the deconvolution. And the RFHV ratio curve is then measured from the radial and ‘Z’ receiver functions. Numerical tests are carried out with the signal-to-noise-ratio of 5 and 10 for radial receiver function. Same inversion parameters are used as in the joint inversion in Fig. 3. From the tests, we find that the joint inversion with a SNR of 10 is quite stable and the result is close to that without noise (Fig. 9a). When the noise level is increased, absolute velocity of the middle crust and the uppermost mantle is affected, as in the inversion with SNR of 5 (Fig. 9d). From the comparison, we find that noise mainly affects the RFHV ratio of larger integral parameter T, while the part of short integral time window is less affected. As a result, the upper crustal velocity is better constrained compared to the deeper part (Fig. 9d). We have to point out that the numerical tests with artificial noise may not be representative of the real noise effect on the Earth. Real noise may be composed by both short and long periods coming from the background. One practical way to suppress the noise is to use more data for stacking of RFWF and linear fitting of RFHV ratio. In addition, we should use the RFHV ratio of short integral time window which is less affected, if the noise level is relatively high. Figure 9. View largeDownload slide Numerical inversion tests with noise. (a–c) Joint inversion of RFHV with RFWF test for the signal to noise ratio of 10. (d–f) Joint inversion of RFHV with RFWF test for the signal-to-noise ratio of 5. Figure 9. View largeDownload slide Numerical inversion tests with noise. (a–c) Joint inversion of RFHV with RFWF test for the signal to noise ratio of 10. (d–f) Joint inversion of RFHV with RFWF test for the signal-to-noise ratio of 5. The joint inversion of RFHV and RFWF may also be improved with surface wave measurements. Surface wave dispersion may help the RFHV ratio to determine the shear-wave velocity in the middle or lower crust, in case the long period RFHV is affected by background noise. And Rayleigh wave ellipiticity may also be added to the joint inversion to help to constrain the shear-wave velocity gradient. In previous study, we found that Rayleigh wave ellipiticity can help determine the smooth long period velocity change with depth, which is complementary to RFWF (Chong et al. 2016). And the joint inversion of Rayleigh wave ellipiticity and RFWF has been applied to obtain a crustal model of the station HYB. From comparing the crustal model of station HYB by Chong et al. (2016) and this study, we found there are some difference in the velocity gradient, as shown in Fig. 8(b). The optimal model of this study does not predict the observed Rayleigh wave ellipiticity very well (Fig. 7c), suggesting that the velocity gradient determination could be improved by Rayleigh wave ellipiticity. The model by Chong et al. (2016) cannot predict the observed RFHV ratio in this study, the predicted RFHV is smaller than data in the range of T < 5 s as in Fig. 8(d). This suggests the absolute velocity of the model by Chong et al. (2016) is slower than the true model, which has been demonstrated in by Chong et al. (2016) and proves that the RFHV ratio can help to determine the absolute shear-wave velocity. 5 CONCLUSIONS We proposed a new seismic measurement: the RFHV ratio. It makes use of the time and frequency dependence of the amplitude information of receiver functions. The sensitivity analysis indicates that it provides integral constraints to the average velocity in the depth range it samples, thus can be jointly inverted with the RFWF. Numerical tests indicate that the RFHV ratio alone can give us a smooth velocity model with the average velocity well determined, but the sharp interfaces are not well resolved. However, the joint inversion of both data sets better constrains both absolute velocity and velocity jumps. Both synthetic and real data applications proves that the RFHV ratio and RFWF can be jointly inverted, and non-uniqueness reduced. The joint inversion can be used for sparse networks and dense arrays, since they are single station measurements and have localized sensitivities. The RFHV ratio might also be useful for studying the sedimentary or subsurface shear-wave velocity structure. From previous studies, the amplitude ratio between the vertical direct P and the radial P wave or sedimentary converted Ps phase can be used to estimate the shallow shear-wave velocity structure (Li et al. 2014; Ni et al. 2014; Kim et al. 2016a). Thus, in a similar way, the receiver function HV ratio technique is possible for determining shallow shear-wave velocity structure. It may also be useful for determining the crustal azimuthal anisotropy by studying azimuthal variation of RFHV ratio (Cristiano et al. 2016). A possible limitation is that, the joint inversion may be affected by the background noise and structural generated noise, as a result, the absolute velocity in deeper part may not be well constrained. So, surface wave measurements, like dispersion or Rayleigh wave ellipiticity, which can sample the deeper structure, may be combined to put tighter constraints to the velocity structure. 5.1 Data and resources The seismic waveform data of the station HYB were obtained from the Incorporated Research Institutions for Seismology Data Management Center (http://ds.iris.edu/ds/nodes/dmc/, last accessed January 2017). The deconvolution code for receiver function is provided by Lupei Zhu of Saint Louis University (http://www.eas.slu.edu/People/LZhu/home.html, last accessed January 2017). The function AWGN and ROBUSTFIT in MATLAB is used to add white Gaussian noise to the receiver function and perform the linear fitting (www.mathworks.com/products/matlab, last accessed January 2017). ACKNOWLEDGEMENTS This study is supported by National Science Foundation of China (Grant Number 41504075, 41461164003, 41474049 and 41661164035). Many thanks to two anonymous reviewers for careful review and constructive suggestions. REFERENCES Ammon C.J., 1991. The isolation of receiver effects from teleseismic P waveforms, Bull. seism. Soc. Am. , 81( 6), 2504– 2510. Ammon C.J., Randall G.E., Zandt G., 1990. On the nonuniqueness of receiver function inversions, J. geophys. Res. , 95( B10), 15 303–15 318. https://doi.org/10.1029/JB095iB10p15303 Google Scholar CrossRef Search ADS   Bodin T., Sambridge M., Tkalčić H., Arroucau P., Gallagher K., Rawlinson N., 2012. Transdimensional inversion of receiver functions and surface wave dispersion, J. geophys. Res. , 117, B02301, doi:10.1029/2011JB008560. https://doi.org/10.1029/2011JB008560 Google Scholar CrossRef Search ADS   Brocher T.M., 2005. Empirical relations between elastic wavespeeds and density in the Earth's crust, Bull. seism. Soc. Am. , 95( 6), 2081– 2092. https://doi.org/10.1785/0120050077 Google Scholar CrossRef Search ADS   Burdick L.J., Langston C.A., 1977. Modeling crustal structure through the use of converted phases in teleseismic body-waveforms, Bull. seism. Soc. Am. , 67, 677– 691. Chang S.-J., 2004. Joint analysis of teleseismic receiver functions and surface wave dispersion using the genetic algorithm, Bull. seism. Soc. Am. , 94( 2), 691– 704. https://doi.org/10.1785/0120030110 Google Scholar CrossRef Search ADS   Chen Y., Niu F., 2013. Ray-parameter based stacking and enhanced pre-conditioning for stable inversion of receiver function data, Geophys. J. Int. , 194, 1682– 1700. Google Scholar CrossRef Search ADS   Chen Y., Niu F., 2016. Joint inversion of receiver functions and surface waves with enhanced preconditioning on densely distributed CNDSN stations: crustal and upper mantle structure beneath China, J. geophys. Res. , 121, 743– 766. Google Scholar CrossRef Search ADS   Chen L., Wen L., Zheng T., 2005. A wave equation migration method for receiver function imaging: 2. Application to the Japan subduction zone, J. geophys. Res. , 110, B11310, doi:10.1029/2005JB003666. Cheng C., Bodin T., Allen R.M., 2016. Three-dimensional pre-stack depth migration of receiver functions with the fast marching method: a Kirchhoff approach, Geophys. J. Int. , 205( 2), 819– 829. https://doi.org/10.1093/gji/ggw062 Google Scholar CrossRef Search ADS   Chong J., Ni S., Zhao L., 2015. Joint inversion of crustal structure with the Rayleigh wave phase velocity dispersion and the ZH ratio, Pure appl. Geophys. , 172( 10), 2585– 2600. https://doi.org/10.1007/s00024-014-0902-z Google Scholar CrossRef Search ADS   Chong J., Ni S., Chu R., Somerville P., 2016. Joint inversion of body-wave receiver function and rayleigh-wave ellipticity, Bull. seism. Soc. Am. , 106( 2), 537– 551. https://doi.org/10.1785/0120150075 Google Scholar CrossRef Search ADS   Cristiano L., Meier T., Krüger F., Keers H., Weidle C., 2016. Teleseismic P-wave polarization analysis at the Gräfenberg array, Geophys. J. Int. , 207( 3), 1456– 1471. https://doi.org/10.1093/gji/ggw339 Google Scholar CrossRef Search ADS   Du Z.J., Foulger G.R., 1999. The crustal structure beneath the northwest fjords, Iceland, from receiver functions and surface waves, Geophys. J. Int. , 139( 2), 419– 432. https://doi.org/10.1046/j.1365-246x.1999.00945.x Google Scholar CrossRef Search ADS   Hannemann K., Krüger F., Dahm T., Lange D., 2016. Oceanic lithospheric S-wave velocities from the analysis of P-wave polarization at the ocean floor, Geophys. J. Int. , 207( 3), 1796– 1817. https://doi.org/10.1093/gji/ggw342 Google Scholar CrossRef Search ADS   Ingber L., 1989. Very fast simulated re-annealing, Math. Comput. Model. , 12( 8), 967– 973. https://doi.org/10.1016/0895-7177(89)90202-1 Google Scholar CrossRef Search ADS   Jacobsen B.H., Svenningsen L., 2008. Enhanced uniqueness and linearity of receiver function inversion, Bull. seism. Soc. Am. , 98( 4), 1756– 1767. https://doi.org/10.1785/0120070180 Google Scholar CrossRef Search ADS   Juliá J., Ammon C.J., Herrmann R.B., Correig A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations, Geophys. J. Int. , 143( 1), 99– 112. https://doi.org/10.1046/j.1365-246x.2000.00217.x Google Scholar CrossRef Search ADS   Juliá J., Jagadeesh S., Rai S.S., Owens T.J., 2009. Deep crustal structure of the Indian shield from joint inversion of P wave receiver functions and Rayleigh wave group velocities: implications for Precambrian crustal evolution, J. geophys. Res. , 114, B10313, doi:10.1029/2008JB006261. https://doi.org/10.1029/2008JB006261 Google Scholar CrossRef Search ADS   Kieling K., Roessler D., Krueger F., 2011. Receiver function study in northern Sumatra and the Malaysian peninsula, J. Seismol ., 15( 2), 235– 259. https://doi.org/10.1007/s10950-010-9222-7 Google Scholar CrossRef Search ADS   Kikuchi M., Kanamori, H., 1982. Inversion of complex body waves, Bull. seism. Soc. Am. , 72, 491– 506. Kim B. et al.  , 2016a. Subsurface shear wave velocity characterization using P-wave seismograms in central and Eastern North America, Earthq. Spectra , 32( 1), 143– 169. https://doi.org/10.1193/123013EQS299M Google Scholar CrossRef Search ADS   Kim S., Dettmer J., Rhie J., Tkalčić H., 2016b. Highly efficient Bayesian joint inversion for receiver-based data and its application to lithospheric structure beneath the southern Korean Peninsula, Geophys. J. Int. , 206( 1), 328– 344. https://doi.org/10.1093/gji/ggw149 Google Scholar CrossRef Search ADS   Kind R., Kosarev G.L., Petersen N.V., 1995. Receiver functions at the stations of the German Regional Seismic Network (GRSN), Geophys. J. Int. , 121( 1), 191– 202. https://doi.org/10.1111/j.1365-246X.1995.tb03520.x Google Scholar CrossRef Search ADS   Kiselev S., Vinnik L., Oreshin S., Gupta S., Rai S.S., Singh A., Kumar M.R., Mohan G., 2008. Lithosphere of the Dharwar craton by joint inversion of P and S receiver functions, Geophys. J. Int. ,, 173, 1106– 1118. Google Scholar CrossRef Search ADS   Krüger F., 1994. Sediment structure at GRF from polarization analysis of P waves of nuclear explosions, Bull. seism. Soc. Am. , 84( 1), 149– 170. 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 Google Scholar CrossRef Search ADS   Last R.J., Nyblade A.A., Langston C.A., Owens T.J., 1997. Crustal structure of the East African Plateau from receiver functions and Rayleigh wave phase velocities, J. geophys. Res. , 102( B11), 24 469–24 483. https://doi.org/10.1029/97JB02156 Google Scholar CrossRef Search ADS   Lawrence J.F., 2004. Combined receiver-function and surface wave phase-velocity inversion using a niching genetic algorithm: Application to Patagonia, Bull. seism. Soc. Am. , 94( 3), 977– 987. https://doi.org/10.1785/0120030172 Google Scholar CrossRef Search ADS   Li X., Li Z., Hao T., Wang S., Xing J., 2017. A multi-frequency receiver function inversion approach for crustal velocity structure, Comput. Geosci. , 102, 45– 55. https://doi.org/10.1016/j.cageo.2017.02.009 Google Scholar CrossRef Search ADS   Li Z., Ni S., Somerville P., 2014. Resolving shallow shear-wave velocity structure beneath station CBN by waveform modeling of the Mw 5.8 mineral, Virginia, Earthquake Sequence, Bull. seism. Soc. Am , 104, 944– 952. https://doi.org/10.1785/0120130190 Google Scholar CrossRef Search ADS   Ligorría J.P., Ammon C.J., 1999. Iterative deconvolution and receiver-function estimation, Bull. seism. Soc. Am. , 89( 5), 944– 952. Moorkamp M., Jones A.G., Fishwick S., 2010. Joint inversion of receiver functions, surface wave dispersion, and magnetotelluric data, J. geophys. Res. , 115, B04318, doi:10.1029/2009JB006369. https://doi.org/10.1029/2009JB006369 Google Scholar CrossRef Search ADS   Ni S., Li Z., Somerville P., 2014. Estimating subsurface shear velocity with radial to vertical ratio of local P waves, Seism. Res. Lett. , 85( 1), 82– 90. https://doi.org/10.1785/0220130128 Google Scholar CrossRef Search ADS   Nuttli O., Whitmore J.D., 1961. An observational determination of the variation of the angle of incidence of P waves with epicentral distance, Bull. seism. Soc. Am. , 51( 2), 269– 276. Owens T.J., Zandt G., Taylor S.R., 1984. Seismic evidence for an ancient rift beneath the Cumberland Plateau, Tennessee: a detailed analysis of broadband teleseismic P waveforms, J. geophys. Res. , 89( B9), 7783– 7795. https://doi.org/10.1029/JB089iB09p07783 Google Scholar CrossRef Search ADS   Özalaybey S., Savage, M., Sheehan, A., Louie, J., Brune, J., 1997. Shear-wave velocity structure in the northern Basin and Range province from the combined analysis of receiver functions and surface waves, Bull. seism. Soc. Am. , 87( 1), 183– 199. Peng H., Hu J., Yang H., Wen L., 2012. An effective technique to constrain the non-uniqueness of receiver function inversion, Chinese J. Geophys. , 55( 2), 194– 205. https://doi.org/10.1002/cjg2.1714 Google Scholar CrossRef Search ADS   Randall G.E., 1989. Efficient calculation of differential seismograms for lithospheric receiver functions, Geophys. J. Int. , 99( 3), 469– 481. https://doi.org/10.1111/j.1365-246X.1989.tb02033.x Google Scholar CrossRef Search ADS   Rondenay S., 2009. Upper mantle imaging with array recordings of converted and scattered teleseismic waves, Surv. Geophys. , 30( 4–5), 377– 405. https://doi.org/10.1007/s10712-009-9071-5 Google Scholar CrossRef Search ADS   Rothman D.H., 1986. Automatic estimation of large residual statics corrections, Geophysics , 51( 2), 332– 346. https://doi.org/10.1190/1.1442092 Google Scholar CrossRef Search ADS   Sambridge M., 1999. Geophysical inversion with a neighbourhood algorithm-I. Searching a parameter space, Geophys. J. Int. , 138( 2), 479– 494. https://doi.org/10.1046/j.1365-246X.1999.00876.x Google Scholar CrossRef Search ADS   Sarkar D., Kumar M.R., Saul J., Kind R., Raju P.S., Chadha R.K., Shukla A.K., 2003. A receiver function perspective of the Dharwar craton (India) crustal structure, Geophys. J. Int. , 154( 1), 205– 211. https://doi.org/10.1046/j.1365-246X.2003.01970.x Google Scholar CrossRef Search ADS   Schiffer C., Jacobsen B.H., Balling N., Ebbing J., Nielsen S.B., 2015. The east Greenland caledonides-teleseismic signature, gravity and isostasy, Geophys. J. Int. , 203( 2), 1400– 1418. https://doi.org/10.1093/gji/ggv373 Google Scholar CrossRef Search ADS   Schiffer C., Stephenson R., Oakey G.N., Jacobsen B.H., 2016. The crustal structure of Ellesmere Island, Arctic Canada-teleseismic mapping across a remote intraplate orogenic belt, Geophys. J. Int. , 204( 3), 1579– 1600. https://doi.org/10.1093/gji/ggv539 Google Scholar CrossRef Search ADS   Shen W., Ritzwoller M.H., Schulte-Pelkum V., Lin F., 2013. Joint inversion of surface wave dispersion and receiver functions: a Bayesian Monte-Carlo approach, Geophys. J. Int. , 192( 2), 807– 836. https://doi.org/10.1093/gji/ggs050 Google Scholar CrossRef Search ADS   Shibutani T., Sambridge M., Kennett B., 1996. Genetic algorithm inversion for receiver functions with application to crust and uppermost mantle structure beneath eastern Australia, Geophys. Res. Lett. , 23( 14), 1829– 1832. https://doi.org/10.1029/96GL01671 Google Scholar CrossRef Search ADS   Svenningsen L., Jacobsen B.H., 2007. Absolute S -velocity estimation from receiver functions, Geophys. J. Int. , 170( 3), 1089– 1094. https://doi.org/10.1111/j.1365-246X.2006.03505.x Google Scholar CrossRef Search ADS   Tkalčić H., Pasyanos M.E., Rodgers A.J., Gök R., Walter W.R., Al‐Amri A., 2006. A multistep approach for joint modeling of surface wave dispersion and teleseismic receiver functions: implications for lithospheric structure of the Arabian Peninsula, J. geophys. Res. , 111, B11311, doi:10.1029/2005JB004130. Google Scholar CrossRef Search ADS   Vinnik L.P., Reigber C., Aleshin I.M., Kosarev G.L., Kaban M.K., Oreshin S.I., Roecker S.W., 2004. Receiver function tomography of the central Tien Shan, Earth planet. Sci. Lett. , 225( 1–2), 131– 146. https://doi.org/10.1016/j.epsl.2004.05.039 Google Scholar CrossRef Search ADS   Zhao L., Sen M.K., Stoffa P., Frohlich C., 1996. Application of very fast simulated annealing to the determination of the crustal structure beneath Tibet, Geophys. J. Int. , 125( 2), 355– 370. https://doi.org/10.1111/j.1365-246X.1996.tb00004.x Google Scholar CrossRef Search ADS   Zhou L., 2000. Seismic properties of the central Indian shield, Bull. seism. Soc. Am. , 90( 5), 1295– 1304. https://doi.org/10.1785/0119990039 Google Scholar CrossRef Search ADS   Zhu L., Kanamori H., 2000. Moho depth variation in southern California from teleseismic receiver functions, J. geophys. Res. , 105( B2), 2969– 2980. https://doi.org/10.1029/1999JB900322 Google Scholar CrossRef Search ADS   Wu Q., Li Y., Zhang R., Zeng R., 2007. Wavelet modelling of broad-band receiver functions, Geophys. J. Int. , 170( 2), 534– 544. https://doi.org/10.1111/j.1365-246X.2007.03467.x Google Scholar CrossRef Search ADS   © 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

# Receiver function HV ratio: a new measurement for reducing non-uniqueness of receiver function waveform inversion

, Volume 212 (2) – Feb 1, 2018
11 pages

Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx464
Publisher site
See Article on Publisher Site

### Abstract

Summary It is known that a receiver function has relatively weak constraint on absolute seismic wave velocity, and that joint inversion of the receiver function with surface wave dispersion has been widely applied to reduce the trade-off of velocity with interface depth. However, some studies indicate that the receiver function itself is capable for determining the absolute shear-wave velocity. In this study, we propose to measure the receiver function HV ratio which takes advantage of the amplitude information of the receiver function to constrain the shear-wave velocity. Numerical analysis indicates that the receiver function HV ratio is sensitive to the average shear-wave velocity in the depth range it samples, and can help to reduce the non-uniqueness of receiver function waveform inversion. A joint inversion scheme has been developed, and both synthetic tests and real data application proved the feasibility of the joint inversion. Joint inversion, Body waves, Crustal structure 1 INTRODUCTION Receiver function analysis is a widely used technique for resolving crustal and upper-mantle structures. It is sensitive to the local structure beneath a seismic station, and is almost insensitive to the earthquake source and propagation effects through the mantle (Burdick & Langston 1977; Lang 1979). It has been applied to obtain the crustal structure beneath a seismic station by waveform inversion (Owens et al. 1984; Randall 1989; Kind et al. 1995; Zhao et al. 1996; Sambridge 1999; Chang et al. 2004; Lawrence & Wiens 2004; Vinik et al. 2004) or image the seismic interfaces by converted phases or multiples found within the receiver functions (Zhu & Kanamori 2000; Chen et al. 2005; Rondenay 2009; Cheng et al. 2016). It has been generally accepted that a receiver function has a non-uniqueness issue in constraining the absolute velocity because of the strong trader-off between the velocity and interface depth (Ammon et al. 1990). A popular way to release the non-uniqueness is through joint inversion of the receiver function waveform (RFWF) and surface wave dispersion (e.g. Last et al. 1997; Özalaybey et al. 1997; Du & Foulger 1999; Juliá et al. 2000; Chang et al. 2004; Lawrence & Wiens 2004; Tkalčić et al. 2006; Moorkamp et al. 2010; Bodin et al. 2012; Shen et al. 2013; Chen & Niu 2016; Kim et al. 2016b). Recently, Chong et al. (2016) tried to invert receiver function with Rayleigh wave ellipticity, and found that the joint inversion had better constraints on the velocity gradient than the absolute velocity. Some studies, however, indicate that the receiver function itself is capable of resolving the absolute shear-wave velocity (Ammon 1991; Svenningsen & Jacobsen 2007) and the non-uniqueness of the waveform inversion can be reduced by using other aspects of receiver function (Wu et al. 2007; Jacobsen & Svenningsen 2008; Li et al. 2017). Wu et al. (2007) and Li et al. (2017) used multifrequency receiver functions to reduce the non-uniqueness of RFWF inversion. Jacobsen & Svenningsen (2008) stratified the velocity model in quasi-continuous delay time to reduce the non-uniqueness. Svenningsen & Jacobsen (2007) used the apparent shear-wave velocity (Vs,app) measured from RFWF to constrain shear-wave velocity structure. Inspired by this Vs,app approach, Kieling et al. (2011) and Peng et al. (2012) developed empirical relations between the Vs,app and the shear-wave velocity structure to construct starting models for the kernel based RFWF inversion, and they claimed that the non-uniqueness was reduced. Schiffer et al. (2015, 2016) applied the joint inversion of the RFWF and Vs,app to study crustal and upper-mantle structure. Hannemann et al. (2016) studied oceanic lithospheric shear-wave velocity structure with Vs,app measured from the P wave recorded by an OBS (ocean bottom seismometer). All these studies show that P wave particle motion is an important parameter for determining crustal structure, and it can help to improve waveform inversion from receiver functions. Although the Vs,app method provides constraints on absolute Vs, it seems to have weak sensitivity to deep discontinuities. Instead, the time-domain RFWF may resolve discontinuities better. A joint inversion with Vs,app and receiver function might improve inversion of both absolute Vs and interface positions. In this study, we propose to measure the time and frequency dependent amplitude information of receiver function to determine the crustal shear-wave velocity. Since the proposed measurement contains the amplitude ratio between the radial and vertical P wave train, and we name it receiver function HV ratio. The receiver function HV ratio contains not only the information of P wave amplitude ratio (particle motion) but also the amplitude information of other phases such as P-to-S conversion and its multiples. We investigated the sensitivity of the receiver function HV ratio to the velocity structure, and the feasibility of joint inversion of the HV ratio and waveform of the receiver function. Both synthetic tests and real data applications indicate that the receiver function HV ratio can be used to reduce the non-uniqueness of RFWF inversion. 2 RECEIVER FUNCTION HV RATIO (RFHV) Studies indicate that P wave particle motion carries velocity structure information beneath a seismic station (Kruger 1994; Svenningsen & Jacobsen 2007; Ni et al. 2014). For a half space model, P wave particle motion on the surface is determined by the ray parameter and the subsurface shear-wave velocity (Svenningsen & Jacobsen 2007), which can be derived from the relation between the apparent and true incident angle of P wave (Nuttli & Whitmore 1961). The relation can be written as following:   $$\mathop V\nolimits_S = \frac{{\sin [{\textstyle{1 \over 2}}\overline {{i_p}} ]}}{p},$$ (1)where $${\bar{i}_p}$$is the apparent incident angle of P wave, and p is the ray parameter. For the teleseismic P wave, the ray parameter can be calculated based on a reference earth model, thus the apparent incident angle or amplitude ratio can be used to determine the subsurface shear-wave velocity. Direct measuring of the apparent incident angle could be difficult due to the complicated seismic P waveform (Kruger 1994). Svenningsen & Jacobsen (2007) proposed a convenient way to measure the apparent incident angle from receiver function. They measure the frequency dependent apparent incident angle of P wave, and derive the apparent shear-wave velocity Vs, app(T) to determine the crustal velocity structure. Here, we propose to measure a transformation of the RFWF, which contains the amplitude ratio between the radial and vertical P wave train. The amplitude ratio is measured for different time window, with a squared-cosine smooth filter (weighting function) centred at t = 0 s following Svenningsen & Jacobsen (2007). It is named the ‘receiver function HV ratio’ (shortened as RFHV hereafter), which can be measured as following:   $${\rm{RFHV(}}T{\rm{)}} = {{{\int_{{{\rm{ - }}T}}^{T}{{{R_{RF}}(\tau ){{\cos }^2}\left({{{\pi \tau } \over {2T}}}\right){\rm{d}}\tau }}} \over {\int_{{ - T}}^{T}{{{Z_{RF}}(\tau ){{\cos }^2}\left({{{\pi \tau } \over {2T}}}\right){\rm{d}}\tau }}}}},$$ (2)where T is half length of the time window, RRF is the radial RFWF. And ZRF is ‘Z’ receiver function (Svenningsen & Jacobsen 2007) obtained from deconvolution of vertical P wave train by itself, which is ideally the Gaussian filtering function used in calculating receiver function. In this equation, the radial receiver function is integrated in a variable time window, with a time window depended squared-cosine weighting function which smooth the waveform for integration. It is then scaled with a scaling factor which is calculated from the ZRF with a integration in the same way. Since in the calculation of receiver function, a Gaussian filtering function is applied, thus the scaling factor helps to remove the effect of the filtering function. The RFHV contains information on P wave particle motion and also amplitude information for P-to-S converted phases or multiples from Moho and other interfaces. Thus we can invert for the velocity structure beneath the seismic station by fitting the RFHV ratio curve, without constructing a starting velocity model for inversion using empirical relations (Kieling et al. 2011; Peng et al. 2012). Before applying RFHV ratio to invert for the crustal structure, we take a one layered crustal model (Fig. 1a) as an example to analyse the characteristics and sensitivity of this new measurement to the velocity structure. Synthetic receiver functions are calculated from the ray parameter range from 0.04 to 0.08 s km–1 with an increment of 0.005 s km–1. The Gaussian filtering factor is chosen to be 2.5, which is a commonly used value for crustal structure study. Here only radial receiver functions are displayed (Fig. 1b). The RFHV ratio curves are then calculated following eq. (2), and are shown in Fig. 1(c). We found that for this synthetic example, the effects of the Gaussian factor on the RFHV ratio is fairly small in the commonly applied range (1.0–4.0) due to the scaling factor in the eq. (2). As we can see in Fig. 1(c), for each ray parameter, the RFHV ratio curve varies mostly in the time window of T < 20 s and varies slightly after this time window. This is because most of the multiples generated in the crust end before 20 s. For the first few seconds, the measurement is related to P wave HV ratio, and each of amplitude jumps on the RFHV ratio curves correspond to the contribution of a new phase coming into the integral time window. Another important feature is that the RFHV ratio increases almost linearly with the ray parameter for a given T (integral parameter), as shown in Fig. 2(a). Here, the RFHV ratios with integral parameter T of 1.0, 6.0, 10.0 and 20.0 s are taken as an example to demonstrate their linear relationship with the ray parameter. Because of this relationship, the real data measurements can be fitted with a line in the RFHV and ray parameter domain for each filter-parameter to depress the anomalies, and RFHV ratio for given ray parameters can be obtained for further analysis. Since the RFHV maybe unevenly sampled in the ray parameter domain, linear fitting of the measurements in RFHV and ray parameter domain helps to avoid the bias, thus reduce the uncertainty of the observations. For sensitivity analysis, we calculate the sensitivity kernel numerically with finite difference method. We find that the RFHV ratio is more sensitive to Vs, and less sensitive to density and Vp/Vs, while being almost insensitive to Vp, as shown in Fig. 2(b). The longer time window measurement is more sensitive to the deeper structure (Fig. 2c). Also, the sensitivity kernel for Vs is mostly positive, which means the RFHV ratio is sensitive to the average shear-wave velocity in the depth range it samples. Figure 1. View largeDownload slide One layered crustal model, synthetic radial RFWFs and the corresponding RFHV ratio curves. (a) The one layered crustal model, with the Moho located at a depth of 30 km, the Vp/Vs ratio is 1.732. (b) The synthetic P wave RFWFs (radial component) with the Gaussian factor of 2.5, displayed according to the ray parameters. (c) The synthetic RFHV ratio curves with ray parameters marked. Figure 1. View largeDownload slide One layered crustal model, synthetic radial RFWFs and the corresponding RFHV ratio curves. (a) The one layered crustal model, with the Moho located at a depth of 30 km, the Vp/Vs ratio is 1.732. (b) The synthetic P wave RFWFs (radial component) with the Gaussian factor of 2.5, displayed according to the ray parameters. (c) The synthetic RFHV ratio curves with ray parameters marked. Figure 2. View largeDownload slide Features of RFHV ratio. (a) Almost linear relation between the RFHV ratio and the ray parameter for each integral parameter T. Examples of four different integral parameters (T = 1.0, 6.0, 10.0 and 20.0 s) are plotted with different symbols. (b) The sensitivity kernels respected to Vp, Vs, density and Vp/Vs ratio for T = 6.0 s. (c) Variation of the sensitivity kernel (Ks) with integral parameter T. Figure 2. View largeDownload slide Features of RFHV ratio. (a) Almost linear relation between the RFHV ratio and the ray parameter for each integral parameter T. Examples of four different integral parameters (T = 1.0, 6.0, 10.0 and 20.0 s) are plotted with different symbols. (b) The sensitivity kernels respected to Vp, Vs, density and Vp/Vs ratio for T = 6.0 s. (c) Variation of the sensitivity kernel (Ks) with integral parameter T. From the analysis, we see that the RFHV ratio is sensitive to the average shear-wave velocity, thus it may be combined with the RFWF to reduce the non-uniqueness of the inversion. 3 JOINT INVERSION OF RFHV AND RFWF As demonstrated in former section, the RFHV ratio is sensitive to the average shear-wave velocity, which could help to reduce the non-uniqueness of RFWF inversion. Here we developed a tool to invert for crustal velocity structure with these two measurements and investigate the feasibility of the joint inversion method. We will describe the inversion method, and present numerical tests and real data application to demonstrate the feasibility of the joint inversion. 3.1 Inversion method The crustal model is parameterized into sublayers with shear-wave velocity (Vs), ratio of P-wave velocity and shear-wave velocity (Vp/Vs), and layer thickness (hi). The density of each layer is scaled to Vp using the relation by Brocher (2005). A subroutine by Shibutani et al. (1996) is modified to calculate the radial and ‘Z’ receiver functions, and the RFHV ratio is calculated following eq. (2) in Section 2. The objective function for fitting both the RFWF and RFHV ratio is defined as following:   \begin{eqnarray} ObjFun &=& {{{{W_{{\rm{RFWF}}}}} \over {{N_{{\rm{RFWF}}}}}}}\sum\limits_{i = 1}^{{N_{{\rm{RFWF}}}}}\nonumber\\ &&{\sqrt {{{1 \over {NPT{S_i}}}}\sum\limits_{j = 1}^{NPT{S_i}} {{{{{{\left( {RFWF_{\rm obs}^i(j) - RFWF_{\rm syn}^i(j)} \right)}^2}} \over {\sigma _{{\rm{RFWF}}}^2(j)}}}} } }\nonumber\\ && + \,{{{{W_{{\rm{RFHV}}}}} \over {{N_{{\rm{RFHV}}}}}}}\sum\limits_{k = 1}^{{N_{{\rm{RFHV}}}}}\nonumber\\ &&{\sqrt {{{1 \over {N{T_k}}}}\sum\limits_{l = 1}^{N{T_k}} {{{{{{\left( {RFHV_{\rm obs}^k(l) - RFHV_{\rm syn}^k(l)} \right)}^2}} \over {\sigma _{{\rm{RFHV}}}^2(l)}}}} } } ,\nonumber\\ && +\, {W_S}*Smoothness \end{eqnarray} (3)where WRFWF and WRFHV are the weights for RFWF and RFHV, respectively, and WS is the weight for model smoothness factor. RFWFobs and RFWFsyn represent observed and synthetic radial receiver functions, while RFHVobs and RFHVsyn denote observed and synthetic receiver function HV ratios. NPTS and NT are the number of points in the RFWF and RFHV ratio. Uncertainties of the RFWF and the RFHV ratio (σRFWF and σRFHV) are used to weight each measurement. To avoid abnormal velocity jumps in the inversion, a smoothness factor is applied to the objective function, which is defined as following,   $$Smoothness = \sqrt {{{1 \over {{N_{\rm layer}}{\rm{ - 1}}}}}\sum\limits_{n = 1}^{{N_{\rm layer}}{\rm{ - 1}}} {{{{{{\left( {{V_j} - {V_{j + 1}}} \right)}^{\rm{2}}}} \over {V_{\rm mean}^2}}}} }$$ (4)Or   $$Smoothness {=} \sqrt {{{1 \over {{N_{\rm layer}}{\rm{ - 2}}}}}\sum\limits_{n = 1}^{{N_{\rm layer}}{\rm{ - 2}}} {{{{{{\left( {{{1 \over 2}}\left( {{V_j} {+} {V_{j {+} {\rm{2}}}} - {V_{j + {\rm{1}}}}} \right)} \right)}^2}} \over {V_{\rm mean}^2}}}} } ,$$ (5)where Nlayer is the number of layers in the model, and Vj is the velocity of the jth layer. Eqs (4) and (5) are for the first and second order velocity changes, respectively, with the former favouring half-space like models, while the latter favours velocity models with gradients (Chong et al. 2016). A fast simulated annealing algorithm is employed as the driver for the inversion. It is a non-linear inversion method, which searches for the optimal solution in a given model space. Just like other Monte Carlo non-linear inversion methods, it only needs forward modelling of each measurement without calculating the sensitivity kernels (Rothman 1986; Ingber 1989). More description of the fast simulated annealing algorithm can be found in section 4 of Chong et al. (2015). 3.2 Synthetic tests In order to test the feasibility of the joint inversion method, we carried out numerical inversion tests before applying it to real data. Similar to the inversions in fig. 4 of Chong et al. (2016), we use a velocity model composed of smoothly varying velocity layers under a low velocity layer for the synthetic inversion tests. Since RFWF may fail to determine the smoothly varying velocity structure (Chong et al. 2016), this kind of model will be a good choice for investigating the stability of the inversion. RFWF and RFHV are calculated with a Gaussian factor of 2.5 and ray parameter of 0.06 s km–1. The uncertainty of the radial receiver function is set to be 5 per cent of its maximum amplitude, and for the RFHV ratio it is set to be 5 per cent of its mean value. Numerical tests are carried out for inversions with only single data set and combination of the two data sets. During the inversions, the Vp/Vs ratio is fixed at 1.732, the same with the value of the true model, in which way only the shear-wave velocity and layer thickness are inverted. As one may expect, the inversion with RFWF alone is able to resolve features of the sharp velocity jumps, however, it cannot determine the absolute shear-wave velocity (Vs), as in Figs 3(a) and (b). Although RFWF is well fitted, the optimal velocity model from the inversion cannot predict the RFHV ratio curve (Fig. 3c). From the comparison, the optimal model predicted RFHV ratio is larger than the synthetic data, which indicates that the inverted optimal model is faster in shear-wave velocity than the true model. On the contrary, the inversion with RFHV ratio alone is able to determine a model which matches the average velocity of the true model, however, the sharp velocity jumps like the thin low velocity layer and the Moho interface are not captured (Figs 3d and e). Consequently, the optimal model inverted with RFHV ratio alone cannot predict RFWF (Fig. 3f). The single data set inversion tests suggest that there is complementarity information between the RFHV and RFWF, and that the combination should do a better job. Indeed, the joint inversion of the two data sets is able to resolve a better crustal model, and both the velocity jumps and the absolute velocity are well determined, as in Figs 3(g)–(i). Figure 3. View largeDownload slide Numerical inversion tests for a smooth variation model with a low velocity layer. (a–c) Inversion with receiver function only. The input model (blue) and optimal model (red) are displayed in (a), while the grey lines indicate the 100 best models from the inversion. (b) Waveform comparison between the data (thick black curve) and the optimal model prediction (red curve). Thin black curves indicate the uncertainty of the data. (c) Comparison between the RFHV ratio data (black) and the optimal model prediction (red). (d–f) Inversion of RFHV ratio only. (d) Same as (a) but for the result of RFHV ratio inversion. (e) Comparison between the data (thick black curve) and the optimal model predicted RFHV ratio (red curve), uncertainty is shown with the black bars. (f) Comparison between the RFWF data (black) and the optimal model prediction (red). (g–i) Joint inversion of the two data sets. (g) Same as (a) but for joint version. (h) Waveform fitting of the RFWF. (i) Fitting of the RFHV ratio curve. Notice that the data and the optimal model prediction overlap with each other. Figure 3. View largeDownload slide Numerical inversion tests for a smooth variation model with a low velocity layer. (a–c) Inversion with receiver function only. The input model (blue) and optimal model (red) are displayed in (a), while the grey lines indicate the 100 best models from the inversion. (b) Waveform comparison between the data (thick black curve) and the optimal model prediction (red curve). Thin black curves indicate the uncertainty of the data. (c) Comparison between the RFHV ratio data (black) and the optimal model prediction (red). (d–f) Inversion of RFHV ratio only. (d) Same as (a) but for the result of RFHV ratio inversion. (e) Comparison between the data (thick black curve) and the optimal model predicted RFHV ratio (red curve), uncertainty is shown with the black bars. (f) Comparison between the RFWF data (black) and the optimal model prediction (red). (g–i) Joint inversion of the two data sets. (g) Same as (a) but for joint version. (h) Waveform fitting of the RFWF. (i) Fitting of the RFHV ratio curve. Notice that the data and the optimal model prediction overlap with each other. 3.3 Application to the Indian craton In this section, we apply this method to station HYB in the Indian craton to obtain the crustal structure. Since the crustal structure beneath HYB is relatively simple and has been extensively studied (Zhou et al. 2000; Sarkar et al. 2003; Kiselev et al. 2008; Julia et al. 2009; Chong et al. 2016; Li et al. 2017), it is a suitable place to investigate the feasibility of the method. As done for the synthetic tests, single data set inversions and two data set joint inversions are performed. 3.3.1 Data processing High quality teleseismic waveforms from 785 earthquakes that occurred between 2000 and 2007 with epicentral distances between 30° and 90° and magnitudes between M6 and M7, were collected for the station HYB (Fig. 4). The instrument responses, trends and means of each component were removed from the raw data. All the waveforms are then selected according to the signal-to-noise ratio of P waves on the three components (ENZ). After the waveform selection, P wave trains were obtained by applying a 100-s time window which starts from 10 s before the onset of P wave. The E–N components were then rotated to radial and tangential components. Receiver functions are calculated with an iterative deconvolution method (Kikuchi & Kanamori 1982; Ligorría & Ammon 1999). Figure 4. View largeDownload slide Station location of HYB in Indian Craton and distribution of the earthquakes. Black triangle for the station and circles for the earthquakes. Figure 4. View largeDownload slide Station location of HYB in Indian Craton and distribution of the earthquakes. Black triangle for the station and circles for the earthquakes. In total, 403 high quality RFWFs were selected for further processing, Fig. 5(a). We first remove the outliers according to the mean and standard deviation of the RFWF, this process helps to remove waveforms with long period noise. Then the waveforms are further selected according to the cross-correlation coefficients between the individual waveform and the mean RFWF. For the waveform inversion, we need to stack the radial receiver functions to reduce the number of records and improve the signal-to-noise ratio. We carried out two-step stacking of the radial receiver functions. For the first step, the raw radial receiver functions are stacked using a ray parameter bin of 0.002 s km–1 with an overlap of 0.002 s km–1. This procedure generates stacked receiver functions that are nearly evenly spaced in the ray parameter domain, avoiding large contributions from some ranges of the ray parameter. In the second step, the radial receiver functions are further processed with a ray parameter based stacking method proposed by Chen & Niu (2013), which corrects for the moveouts and amplitudes for all ray parameters based on a one layer crustal model obtained with H–K stacking method (Zhu & Kanamori 2000). And the method is modified to calculate the uncertainties for the stacked waveform, which will be used for the data weighting during the inversion. The radial receiver functions are stacked to three reference ray parameters (0.05, 0.06 and 0.07 s km–1) with the ray parameter based stacking method (Fig. 5b). Figure 5. View largeDownload slide RFWF and RFHV ratio curves for station HYB. (a) Radial RFWF measured from the records of earthquakes in Fig. 4. (b) Stacked radial RFWFs with ray parameter correction. Uncertainties of the stacked waveform are shown with thin lines. (c) RFHV ratio curves. (d) RFHV ratio curves for three ray parameters (0.05, 0.06 and 0.07 s km−1) calculated from the linear fitting of the measurements in (c), the uncertainties are shown with vertical bars. Details of the linear fitting can be found in Fig. 6. Figure 5. View largeDownload slide RFWF and RFHV ratio curves for station HYB. (a) Radial RFWF measured from the records of earthquakes in Fig. 4. (b) Stacked radial RFWFs with ray parameter correction. Uncertainties of the stacked waveform are shown with thin lines. (c) RFHV ratio curves. (d) RFHV ratio curves for three ray parameters (0.05, 0.06 and 0.07 s km−1) calculated from the linear fitting of the measurements in (c), the uncertainties are shown with vertical bars. Details of the linear fitting can be found in Fig. 6. RFHV ratios are computed from the raw radial receiver functions shown in Fig. 5(a) according to eq. (2) in Section 2. After eliminating the outliers, we obtained 252 RFHV ratio curves, as shown in Fig. 5(c). The RFHV ratios are then fitted in the ray parameter domain with a line for each integral parameter (T) according to the linear relation between the RFHV ratio and the ray parameter, as shown in Figs 6(a)–(c). And we used the RFHV values on the fitted line of three ray parameters (0.05, 0.06 and 0.07 s km−1) for inversion (Fig. 5d). Figure 6. View largeDownload slide (a–c) Linear fitting of the RFHV ratios measured for station HYB. Linear fitting of RFHV ratios with integral parameters of T = 6, 10 and 20 seconds are shown in (a–c), respectively. Open triangles for data, thick black lines for fitted lines of RFHV ratio of each integral parameter. Thin grey lines indicate the one standard deviation of linear fitting, and half of the standard deviation is used for the uncertainty of measurement. The linear fitting is carried out with the function ROBUSTFIT in MATLAB. It performs for a robust multilinear regression, and in the data processing we used default algorithm with iteratively reweighted least squares with a bi-square weighting function. Figure 6. View largeDownload slide (a–c) Linear fitting of the RFHV ratios measured for station HYB. Linear fitting of RFHV ratios with integral parameters of T = 6, 10 and 20 seconds are shown in (a–c), respectively. Open triangles for data, thick black lines for fitted lines of RFHV ratio of each integral parameter. Thin grey lines indicate the one standard deviation of linear fitting, and half of the standard deviation is used for the uncertainty of measurement. The linear fitting is carried out with the function ROBUSTFIT in MATLAB. It performs for a robust multilinear regression, and in the data processing we used default algorithm with iteratively reweighted least squares with a bi-square weighting function. 3.3.2 Inversion results Inversions with single data sets and the combination of the two data sets are carried out. From the inversion of only RFWF, we get an optimal model with the Moho located at the depth of about 29 km, and the uppermost mantle shear-wave velocity of about 4.0 km s–1 (Fig. 7a). However, the Moho depth and shear-wave velocity of the uppermost mantle of the optimal model are quite different from previous studies (Zhou et al. 2000; Sarkar et al. 2003; Kiselev et al. 2008; Julia et al. 2009; Li et al. 2017). Although the optimal model is able to fit the RFWF, the RFHV ratio curves are not well-predicted (Figs 7b and c). From the comparison, the amplitude of predicted RFHV ratio is smaller than the data, which means that the shear-wave velocity of the optimal model is slower than that of the real model. The inversion with only RFHV ratio curves gives us a model in which velocity increases smoothly with the depth (Figs 7d and e), but the Moho interface has not been captured. Thus the RFWF, especially the P-to-S conversion and multiples from the Moho, are not well predicted (Fig. 7f). The joint inversion of these two data sets, gives us a crustal model with a Moho located at depth of about 32 km, and the shear-wave velocity increasing smoothly from about 3.5 km s–1 near the surface to roughly 3.8 km s–1 above the Moho (Figs 7g–i). From the comparison with other studies (Sarkar et al. 2003; Kiselev et al. 2008; Zhou et al. 2000; Julia et al. 2009; Li et al. 2017), we found the crustal model from the joint inversion of RFWF and RFHV is similar to other models in the lower crust, but slightly slower in the upper and middle crust than most previous velocity models except for the one by Kiselev et al. (2008) (Fig. 8(a)). The shear-wave velocity of the uppermost mantle is about 4.3 km s–1, which is close to the model by Li et al. (2017), but a little bit slower than others. Figure 7. View largeDownload slide Inversions for the station HYB in Indian Craton. (a–c) Inversion with RFWF alone. (d–f) Inversion with RFHV ratio alone. (g–i) Jointly inversion of the two data sets. Figure 7. View largeDownload slide Inversions for the station HYB in Indian Craton. (a–c) Inversion with RFWF alone. (d–f) Inversion with RFHV ratio alone. (g–i) Jointly inversion of the two data sets. Figure 8. View largeDownload slide Figures for comparing the results of this study with previous studies. (a) Crustal model comparison of the station HYB in Indian Craton. (b) Crustal model comparison with the one by Chong et al. (2016) from the joint inversion of RFWF and Rayleigh wave ellipiticity (SWZH). (c) Comparison between the observed Rayleigh wave ellipiticity (black) and prediction (red) by the optimal model of this study. (d) Comparison between the observed RFHV ratio curves (black) and prediction (red) by the model of Chong et al. (2016). Figure 8. View largeDownload slide Figures for comparing the results of this study with previous studies. (a) Crustal model comparison of the station HYB in Indian Craton. (b) Crustal model comparison with the one by Chong et al. (2016) from the joint inversion of RFWF and Rayleigh wave ellipiticity (SWZH). (c) Comparison between the observed Rayleigh wave ellipiticity (black) and prediction (red) by the optimal model of this study. (d) Comparison between the observed RFHV ratio curves (black) and prediction (red) by the model of Chong et al. (2016). 4 DISCUSSION The inversion of RFWF usually has difficulties in constraining the absolute shear-wave velocity. The RFHV technique, which is based on a integral transformation of the RFWFs, provides a way to constrain the absolute shear-wave velocity. Sensitivity analysis and numerical inversions indicate that RFHV provides constraints on the average shear-wave velocity (Figs 2b–c and Figs 3d–e). Since RFWF provides constraints to the sharp interfaces, the RFHV and RFWF can be combined to put tighter constraints on velocity. The numerical inversion tests and real data application indicate that joint inversion is more reliable than inversions with single data set (Figs 3 and 7). Both the RFWF and RFHV ratio are single station measurements with localized sensitivity, and the joint inversion can be applied to regions with either a sparse network or dense array. In real data applications, the joint inversion may be affected by the background noise. We carried out two numerical tests to investigate the effect of noise on the joint inversion of RFHV and RFWF. The crustal model in Fig. 3 is used for the tests. Radial and vertical receiver functions are calculated with a Gaussian factor of 2.5 and ray parameter of 0.06 s km–1. White Gaussian noise is generated with the function AWGN in MATLAB (see data and resources) and added to the radial receiver function but not to the ‘Z’ receiver function, because the noise has been cancelled on the vertical component during the deconvolution. And the RFHV ratio curve is then measured from the radial and ‘Z’ receiver functions. Numerical tests are carried out with the signal-to-noise-ratio of 5 and 10 for radial receiver function. Same inversion parameters are used as in the joint inversion in Fig. 3. From the tests, we find that the joint inversion with a SNR of 10 is quite stable and the result is close to that without noise (Fig. 9a). When the noise level is increased, absolute velocity of the middle crust and the uppermost mantle is affected, as in the inversion with SNR of 5 (Fig. 9d). From the comparison, we find that noise mainly affects the RFHV ratio of larger integral parameter T, while the part of short integral time window is less affected. As a result, the upper crustal velocity is better constrained compared to the deeper part (Fig. 9d). We have to point out that the numerical tests with artificial noise may not be representative of the real noise effect on the Earth. Real noise may be composed by both short and long periods coming from the background. One practical way to suppress the noise is to use more data for stacking of RFWF and linear fitting of RFHV ratio. In addition, we should use the RFHV ratio of short integral time window which is less affected, if the noise level is relatively high. Figure 9. View largeDownload slide Numerical inversion tests with noise. (a–c) Joint inversion of RFHV with RFWF test for the signal to noise ratio of 10. (d–f) Joint inversion of RFHV with RFWF test for the signal-to-noise ratio of 5. Figure 9. View largeDownload slide Numerical inversion tests with noise. (a–c) Joint inversion of RFHV with RFWF test for the signal to noise ratio of 10. (d–f) Joint inversion of RFHV with RFWF test for the signal-to-noise ratio of 5. The joint inversion of RFHV and RFWF may also be improved with surface wave measurements. Surface wave dispersion may help the RFHV ratio to determine the shear-wave velocity in the middle or lower crust, in case the long period RFHV is affected by background noise. And Rayleigh wave ellipiticity may also be added to the joint inversion to help to constrain the shear-wave velocity gradient. In previous study, we found that Rayleigh wave ellipiticity can help determine the smooth long period velocity change with depth, which is complementary to RFWF (Chong et al. 2016). And the joint inversion of Rayleigh wave ellipiticity and RFWF has been applied to obtain a crustal model of the station HYB. From comparing the crustal model of station HYB by Chong et al. (2016) and this study, we found there are some difference in the velocity gradient, as shown in Fig. 8(b). The optimal model of this study does not predict the observed Rayleigh wave ellipiticity very well (Fig. 7c), suggesting that the velocity gradient determination could be improved by Rayleigh wave ellipiticity. The model by Chong et al. (2016) cannot predict the observed RFHV ratio in this study, the predicted RFHV is smaller than data in the range of T < 5 s as in Fig. 8(d). This suggests the absolute velocity of the model by Chong et al. (2016) is slower than the true model, which has been demonstrated in by Chong et al. (2016) and proves that the RFHV ratio can help to determine the absolute shear-wave velocity. 5 CONCLUSIONS We proposed a new seismic measurement: the RFHV ratio. It makes use of the time and frequency dependence of the amplitude information of receiver functions. The sensitivity analysis indicates that it provides integral constraints to the average velocity in the depth range it samples, thus can be jointly inverted with the RFWF. Numerical tests indicate that the RFHV ratio alone can give us a smooth velocity model with the average velocity well determined, but the sharp interfaces are not well resolved. However, the joint inversion of both data sets better constrains both absolute velocity and velocity jumps. Both synthetic and real data applications proves that the RFHV ratio and RFWF can be jointly inverted, and non-uniqueness reduced. The joint inversion can be used for sparse networks and dense arrays, since they are single station measurements and have localized sensitivities. The RFHV ratio might also be useful for studying the sedimentary or subsurface shear-wave velocity structure. From previous studies, the amplitude ratio between the vertical direct P and the radial P wave or sedimentary converted Ps phase can be used to estimate the shallow shear-wave velocity structure (Li et al. 2014; Ni et al. 2014; Kim et al. 2016a). Thus, in a similar way, the receiver function HV ratio technique is possible for determining shallow shear-wave velocity structure. It may also be useful for determining the crustal azimuthal anisotropy by studying azimuthal variation of RFHV ratio (Cristiano et al. 2016). A possible limitation is that, the joint inversion may be affected by the background noise and structural generated noise, as a result, the absolute velocity in deeper part may not be well constrained. So, surface wave measurements, like dispersion or Rayleigh wave ellipiticity, which can sample the deeper structure, may be combined to put tighter constraints to the velocity structure. 5.1 Data and resources The seismic waveform data of the station HYB were obtained from the Incorporated Research Institutions for Seismology Data Management Center (http://ds.iris.edu/ds/nodes/dmc/, last accessed January 2017). The deconvolution code for receiver function is provided by Lupei Zhu of Saint Louis University (http://www.eas.slu.edu/People/LZhu/home.html, last accessed January 2017). The function AWGN and ROBUSTFIT in MATLAB is used to add white Gaussian noise to the receiver function and perform the linear fitting (www.mathworks.com/products/matlab, last accessed January 2017). ACKNOWLEDGEMENTS This study is supported by National Science Foundation of China (Grant Number 41504075, 41461164003, 41474049 and 41661164035). Many thanks to two anonymous reviewers for careful review and constructive suggestions. REFERENCES Ammon C.J., 1991. The isolation of receiver effects from teleseismic P waveforms, Bull. seism. Soc. Am. , 81( 6), 2504– 2510. Ammon C.J., Randall G.E., Zandt G., 1990. On the nonuniqueness of receiver function inversions, J. geophys. Res. , 95( B10), 15 303–15 318. https://doi.org/10.1029/JB095iB10p15303 Google Scholar CrossRef Search ADS   Bodin T., Sambridge M., Tkalčić H., Arroucau P., Gallagher K., Rawlinson N., 2012. Transdimensional inversion of receiver functions and surface wave dispersion, J. geophys. Res. , 117, B02301, doi:10.1029/2011JB008560. https://doi.org/10.1029/2011JB008560 Google Scholar CrossRef Search ADS   Brocher T.M., 2005. Empirical relations between elastic wavespeeds and density in the Earth's crust, Bull. seism. Soc. Am. , 95( 6), 2081– 2092. https://doi.org/10.1785/0120050077 Google Scholar CrossRef Search ADS   Burdick L.J., Langston C.A., 1977. Modeling crustal structure through the use of converted phases in teleseismic body-waveforms, Bull. seism. Soc. Am. , 67, 677– 691. Chang S.-J., 2004. Joint analysis of teleseismic receiver functions and surface wave dispersion using the genetic algorithm, Bull. seism. Soc. Am. , 94( 2), 691– 704. https://doi.org/10.1785/0120030110 Google Scholar CrossRef Search ADS   Chen Y., Niu F., 2013. Ray-parameter based stacking and enhanced pre-conditioning for stable inversion of receiver function data, Geophys. J. Int. , 194, 1682– 1700. Google Scholar CrossRef Search ADS   Chen Y., Niu F., 2016. Joint inversion of receiver functions and surface waves with enhanced preconditioning on densely distributed CNDSN stations: crustal and upper mantle structure beneath China, J. geophys. Res. , 121, 743– 766. Google Scholar CrossRef Search ADS   Chen L., Wen L., Zheng T., 2005. A wave equation migration method for receiver function imaging: 2. Application to the Japan subduction zone, J. geophys. Res. , 110, B11310, doi:10.1029/2005JB003666. Cheng C., Bodin T., Allen R.M., 2016. Three-dimensional pre-stack depth migration of receiver functions with the fast marching method: a Kirchhoff approach, Geophys. J. Int. , 205( 2), 819– 829. https://doi.org/10.1093/gji/ggw062 Google Scholar CrossRef Search ADS   Chong J., Ni S., Zhao L., 2015. Joint inversion of crustal structure with the Rayleigh wave phase velocity dispersion and the ZH ratio, Pure appl. Geophys. , 172( 10), 2585– 2600. https://doi.org/10.1007/s00024-014-0902-z Google Scholar CrossRef Search ADS   Chong J., Ni S., Chu R., Somerville P., 2016. Joint inversion of body-wave receiver function and rayleigh-wave ellipticity, Bull. seism. Soc. Am. , 106( 2), 537– 551. https://doi.org/10.1785/0120150075 Google Scholar CrossRef Search ADS   Cristiano L., Meier T., Krüger F., Keers H., Weidle C., 2016. Teleseismic P-wave polarization analysis at the Gräfenberg array, Geophys. J. Int. , 207( 3), 1456– 1471. https://doi.org/10.1093/gji/ggw339 Google Scholar CrossRef Search ADS   Du Z.J., Foulger G.R., 1999. The crustal structure beneath the northwest fjords, Iceland, from receiver functions and surface waves, Geophys. J. Int. , 139( 2), 419– 432. https://doi.org/10.1046/j.1365-246x.1999.00945.x Google Scholar CrossRef Search ADS   Hannemann K., Krüger F., Dahm T., Lange D., 2016. Oceanic lithospheric S-wave velocities from the analysis of P-wave polarization at the ocean floor, Geophys. J. Int. , 207( 3), 1796– 1817. https://doi.org/10.1093/gji/ggw342 Google Scholar CrossRef Search ADS   Ingber L., 1989. Very fast simulated re-annealing, Math. Comput. Model. , 12( 8), 967– 973. https://doi.org/10.1016/0895-7177(89)90202-1 Google Scholar CrossRef Search ADS   Jacobsen B.H., Svenningsen L., 2008. Enhanced uniqueness and linearity of receiver function inversion, Bull. seism. Soc. Am. , 98( 4), 1756– 1767. https://doi.org/10.1785/0120070180 Google Scholar CrossRef Search ADS   Juliá J., Ammon C.J., Herrmann R.B., Correig A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations, Geophys. J. Int. , 143( 1), 99– 112. https://doi.org/10.1046/j.1365-246x.2000.00217.x Google Scholar CrossRef Search ADS   Juliá J., Jagadeesh S., Rai S.S., Owens T.J., 2009. Deep crustal structure of the Indian shield from joint inversion of P wave receiver functions and Rayleigh wave group velocities: implications for Precambrian crustal evolution, J. geophys. Res. , 114, B10313, doi:10.1029/2008JB006261. https://doi.org/10.1029/2008JB006261 Google Scholar CrossRef Search ADS   Kieling K., Roessler D., Krueger F., 2011. Receiver function study in northern Sumatra and the Malaysian peninsula, J. Seismol ., 15( 2), 235– 259. https://doi.org/10.1007/s10950-010-9222-7 Google Scholar CrossRef Search ADS   Kikuchi M., Kanamori, H., 1982. Inversion of complex body waves, Bull. seism. Soc. Am. , 72, 491– 506. Kim B. et al.  , 2016a. Subsurface shear wave velocity characterization using P-wave seismograms in central and Eastern North America, Earthq. Spectra , 32( 1), 143– 169. https://doi.org/10.1193/123013EQS299M Google Scholar CrossRef Search ADS   Kim S., Dettmer J., Rhie J., Tkalčić H., 2016b. Highly efficient Bayesian joint inversion for receiver-based data and its application to lithospheric structure beneath the southern Korean Peninsula, Geophys. J. Int. , 206( 1), 328– 344. https://doi.org/10.1093/gji/ggw149 Google Scholar CrossRef Search ADS   Kind R., Kosarev G.L., Petersen N.V., 1995. Receiver functions at the stations of the German Regional Seismic Network (GRSN), Geophys. J. Int. , 121( 1), 191– 202. https://doi.org/10.1111/j.1365-246X.1995.tb03520.x Google Scholar CrossRef Search ADS   Kiselev S., Vinnik L., Oreshin S., Gupta S., Rai S.S., Singh A., Kumar M.R., Mohan G., 2008. Lithosphere of the Dharwar craton by joint inversion of P and S receiver functions, Geophys. J. Int. ,, 173, 1106– 1118. Google Scholar CrossRef Search ADS   Krüger F., 1994. Sediment structure at GRF from polarization analysis of P waves of nuclear explosions, Bull. seism. Soc. Am. , 84( 1), 149– 170. 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 Google Scholar CrossRef Search ADS   Last R.J., Nyblade A.A., Langston C.A., Owens T.J., 1997. Crustal structure of the East African Plateau from receiver functions and Rayleigh wave phase velocities, J. geophys. Res. , 102( B11), 24 469–24 483. https://doi.org/10.1029/97JB02156 Google Scholar CrossRef Search ADS   Lawrence J.F., 2004. Combined receiver-function and surface wave phase-velocity inversion using a niching genetic algorithm: Application to Patagonia, Bull. seism. Soc. Am. , 94( 3), 977– 987. https://doi.org/10.1785/0120030172 Google Scholar CrossRef Search ADS   Li X., Li Z., Hao T., Wang S., Xing J., 2017. A multi-frequency receiver function inversion approach for crustal velocity structure, Comput. Geosci. , 102, 45– 55. https://doi.org/10.1016/j.cageo.2017.02.009 Google Scholar CrossRef Search ADS   Li Z., Ni S., Somerville P., 2014. Resolving shallow shear-wave velocity structure beneath station CBN by waveform modeling of the Mw 5.8 mineral, Virginia, Earthquake Sequence, Bull. seism. Soc. Am , 104, 944– 952. https://doi.org/10.1785/0120130190 Google Scholar CrossRef Search ADS   Ligorría J.P., Ammon C.J., 1999. Iterative deconvolution and receiver-function estimation, Bull. seism. Soc. Am. , 89( 5), 944– 952. Moorkamp M., Jones A.G., Fishwick S., 2010. Joint inversion of receiver functions, surface wave dispersion, and magnetotelluric data, J. geophys. Res. , 115, B04318, doi:10.1029/2009JB006369. https://doi.org/10.1029/2009JB006369 Google Scholar CrossRef Search ADS   Ni S., Li Z., Somerville P., 2014. Estimating subsurface shear velocity with radial to vertical ratio of local P waves, Seism. Res. Lett. , 85( 1), 82– 90. https://doi.org/10.1785/0220130128 Google Scholar CrossRef Search ADS   Nuttli O., Whitmore J.D., 1961. An observational determination of the variation of the angle of incidence of P waves with epicentral distance, Bull. seism. Soc. Am. , 51( 2), 269– 276. Owens T.J., Zandt G., Taylor S.R., 1984. Seismic evidence for an ancient rift beneath the Cumberland Plateau, Tennessee: a detailed analysis of broadband teleseismic P waveforms, J. geophys. Res. , 89( B9), 7783– 7795. https://doi.org/10.1029/JB089iB09p07783 Google Scholar CrossRef Search ADS   Özalaybey S., Savage, M., Sheehan, A., Louie, J., Brune, J., 1997. Shear-wave velocity structure in the northern Basin and Range province from the combined analysis of receiver functions and surface waves, Bull. seism. Soc. Am. , 87( 1), 183– 199. Peng H., Hu J., Yang H., Wen L., 2012. An effective technique to constrain the non-uniqueness of receiver function inversion, Chinese J. Geophys. , 55( 2), 194– 205. https://doi.org/10.1002/cjg2.1714 Google Scholar CrossRef Search ADS   Randall G.E., 1989. Efficient calculation of differential seismograms for lithospheric receiver functions, Geophys. J. Int. , 99( 3), 469– 481. https://doi.org/10.1111/j.1365-246X.1989.tb02033.x Google Scholar CrossRef Search ADS   Rondenay S., 2009. Upper mantle imaging with array recordings of converted and scattered teleseismic waves, Surv. Geophys. , 30( 4–5), 377– 405. https://doi.org/10.1007/s10712-009-9071-5 Google Scholar CrossRef Search ADS   Rothman D.H., 1986. Automatic estimation of large residual statics corrections, Geophysics , 51( 2), 332– 346. https://doi.org/10.1190/1.1442092 Google Scholar CrossRef Search ADS   Sambridge M., 1999. Geophysical inversion with a neighbourhood algorithm-I. Searching a parameter space, Geophys. J. Int. , 138( 2), 479– 494. https://doi.org/10.1046/j.1365-246X.1999.00876.x Google Scholar CrossRef Search ADS   Sarkar D., Kumar M.R., Saul J., Kind R., Raju P.S., Chadha R.K., Shukla A.K., 2003. A receiver function perspective of the Dharwar craton (India) crustal structure, Geophys. J. Int. , 154( 1), 205– 211. https://doi.org/10.1046/j.1365-246X.2003.01970.x Google Scholar CrossRef Search ADS   Schiffer C., Jacobsen B.H., Balling N., Ebbing J., Nielsen S.B., 2015. The east Greenland caledonides-teleseismic signature, gravity and isostasy, Geophys. J. Int. , 203( 2), 1400– 1418. https://doi.org/10.1093/gji/ggv373 Google Scholar CrossRef Search ADS   Schiffer C., Stephenson R., Oakey G.N., Jacobsen B.H., 2016. The crustal structure of Ellesmere Island, Arctic Canada-teleseismic mapping across a remote intraplate orogenic belt, Geophys. J. Int. , 204( 3), 1579– 1600. https://doi.org/10.1093/gji/ggv539 Google Scholar CrossRef Search ADS   Shen W., Ritzwoller M.H., Schulte-Pelkum V., Lin F., 2013. Joint inversion of surface wave dispersion and receiver functions: a Bayesian Monte-Carlo approach, Geophys. J. Int. , 192( 2), 807– 836. https://doi.org/10.1093/gji/ggs050 Google Scholar CrossRef Search ADS   Shibutani T., Sambridge M., Kennett B., 1996. Genetic algorithm inversion for receiver functions with application to crust and uppermost mantle structure beneath eastern Australia, Geophys. Res. Lett. , 23( 14), 1829– 1832. https://doi.org/10.1029/96GL01671 Google Scholar CrossRef Search ADS   Svenningsen L., Jacobsen B.H., 2007. Absolute S -velocity estimation from receiver functions, Geophys. J. Int. , 170( 3), 1089– 1094. https://doi.org/10.1111/j.1365-246X.2006.03505.x Google Scholar CrossRef Search ADS   Tkalčić H., Pasyanos M.E., Rodgers A.J., Gök R., Walter W.R., Al‐Amri A., 2006. A multistep approach for joint modeling of surface wave dispersion and teleseismic receiver functions: implications for lithospheric structure of the Arabian Peninsula, J. geophys. Res. , 111, B11311, doi:10.1029/2005JB004130. Google Scholar CrossRef Search ADS   Vinnik L.P., Reigber C., Aleshin I.M., Kosarev G.L., Kaban M.K., Oreshin S.I., Roecker S.W., 2004. Receiver function tomography of the central Tien Shan, Earth planet. Sci. Lett. , 225( 1–2), 131– 146. https://doi.org/10.1016/j.epsl.2004.05.039 Google Scholar CrossRef Search ADS   Zhao L., Sen M.K., Stoffa P., Frohlich C., 1996. Application of very fast simulated annealing to the determination of the crustal structure beneath Tibet, Geophys. J. Int. , 125( 2), 355– 370. https://doi.org/10.1111/j.1365-246X.1996.tb00004.x Google Scholar CrossRef Search ADS   Zhou L., 2000. Seismic properties of the central Indian shield, Bull. seism. Soc. Am. , 90( 5), 1295– 1304. https://doi.org/10.1785/0119990039 Google Scholar CrossRef Search ADS   Zhu L., Kanamori H., 2000. Moho depth variation in southern California from teleseismic receiver functions, J. geophys. Res. , 105( B2), 2969– 2980. https://doi.org/10.1029/1999JB900322 Google Scholar CrossRef Search ADS   Wu Q., Li Y., Zhang R., Zeng R., 2007. Wavelet modelling of broad-band receiver functions, Geophys. J. Int. , 170( 2), 534– 544. https://doi.org/10.1111/j.1365-246X.2007.03467.x Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 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
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### 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. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations