Waveform Fitting of Receiver Functions for Enhanced Retrieval of Crustal Structure in the Presence of Sediments

The receiver function technique is widely used to image crustal structure using P‐to‐S converted phases at the Moho discontinuity. However, the presence of sedimentary layer generates additional P‐to‐S conversions and reverberations, which can overprint the Moho phases and pose problems in imaging crustal structure with standard receiver function techniques. We introduce a robust two‐step method that uses H‐κ stacking to determine average thickness and Vp/Vs of the sedimentary layer, followed by waveform‐fitting of the observed receiver function to constrain the average crustal thickness and sub‐sediment Vp/Vs. We tested the method using both synthetic data and real‐data from stations located on sedimentary layers in the Netherlands and USA. We show that the new method outperforms other common approaches in retrieving accurate Moho depth and sub‐sediment Vp/Vs estimates, even in cases where the Moho phases are completely overprinted by large‐amplitude phases related to sedimentary layers.


Introduction
Receiver function (RF) analysis is a well-established technique to image impedance contrasts in the crust by using isolated P-to-S wave conversions and the reverberations generated at such discontinuities.The velocity contrast at the crust-mantle boundary (Moho) represents a major discontinuity that RFs are highly sensitive to, making RFs useful to constrain average crustal thickness and Vp/Vs (e.g., Ammon, 1991;Langston, 1979;Owens et al., 1984).However, the presence of low-velocity sedimentary layers can obscure Moho signals due to the additional P-to-S wave conversions (PbS) and associated reverberated phases created at the sediment-basement discontinuity.These additional intra-crustal phases can have large amplitudes and similar arrival times to Moho signals, which makes it difficult to retrieve Moho information using standard RF techniques (e.g., van der Meijde et al., 2003;Yeck et al., 2013;Yu et al., 2015;Zelt & Ellis, 1999;Zhu & Kanamori, 2000).
Several methods have been developed to deal with the complex RFs obtained in sedimentary basins.Iterative trial-and-error fitting of observed RFs with synthetic ones (i.e., forward modeling) is a common approach (e.g., Sheehan et al., 1995;Zelt & Ellis, 1999).However, the presence of near-surface reverberations can mask deeper intra-crustal and Moho signals.This expands the range of possible models and aggravates the intrinsic ambiguity of finding a unique model via forward modeling.The wavefield continuation and decomposition technique of Langston (2011), which was improved by Tao et al. (2014) retrieves reliable crustal structures provided that good a priori information on Vp and density is available for the sediments, the sub-sediment units and the uppermost mantle.Since access to reliable a priori information is rarely possible, other studies proposed the application of band-pass or resonance filters directly to the RFs to remove the sediment signals (e.g., Leahy et al., 2012;Reed et al., 2014;Yu et al., 2015).Despite the advantage of working with "filtered" RFs, these methods tend to be casedependent and it is not always obvious what frequency band to remove and/or how to guarantee that the sediment functions for stations overlying a sedimentary layer • The approach is effective even when Moho phases are completely overprinted by positive or negative amplitudes from the sedimentary layer • Synthetic tests and real-data applications show the advantages of the developed approach over common methods for addressing sediment impact Supporting Information: Supporting Information may be found in the online version of this article.
signals have been completely eliminated from the RFs.For instance, Yu et al. (2015) showed that resonance filtering performs well in cases where the reverberations due to the sediment are strong, contain a single dominant frequency, and have a couple of cycles.Furthermore, the effectiveness of the resonance filter degrades when the delay in the PbS arrival is large (e.g., due to a thick or very low-velocity sediment layer), when the impedance contrast at the base the sediment layer is low or when sediment layer is very thin (<one fourth of the wavelength of a PbS phase) or very thick (∼5 km or thicker) (Yu et al., 2015).Zhang and Olugboji (2023) further extended the reverberation filtering technique in a data-driven approach for cases where two or more types of reverberant layers are present (sediment, ocean and glaciers) using homomorphic analysis.A slightly different, two-step approach was proposed by Yeck et al. (2013) based on the original work of Zhu and Kanamori (2000).The average thickness and Vp/Vs of the sediment layer is determined in the first step from H-κ stacking of highfrequency RFs.These are then used in a second step to correct for the time delay in deeper conversions before stacking the Moho phases.The simplicity and effectiveness of this method is attractive.However, its original implementation corrects for timing delays of phases due to slow(er) sediments and is limited to the case where the Moho phases are not overprinted by positive or negative amplitudes of sediment phases as shown by Yeck et al. (2013).
Here, we present an alternative approach for imaging the average crustal thickness and Vp/Vs beneath stations overlying sedimentary layer using RFs.It is a two-step method that integrates H-κ stacking for sediment detection and characterization with waveform-fitting via grid search to characterize sub-sediment crustal layer.The technique relies on the coherence between the Moho and sediment converted phases and the respective reverberations (van der Meijde et al., 2003) and is independent of the thickness, velocity, and strength of the reverberation of the sediment.It adopts a combined misfit criteria integrating the correlation coefficient, rootmean-squared error and phase goodness-of-fit measures to determine two-layer synthetic RF that best fits the observed RF.Because the solution of the sedimentary layer (Vp/Vs and thickness) is determined in the first step, fewer free parameters are solved for in the waveform-fitting step, which helps in constraining the non-uniqueness of the solution (Ammon et al., 1990).It overcomes the problem of negative interference of sediment phases with Moho phases, and does not depend on subjective choices of specific filtering to be applied.In the following sections, we describe the method and demonstrate its effectiveness and limitations using both synthetic tests and real-data.The latter comes from stations located on top of sedimentary layers with different thicknesses deployed as part of the Network of Autonomously Recording Seismographs-Netherlands (NARS-NL), the Netherlands Seismic and Acoustic Network operated by Royal Netherlands Meteorological Institute (NSAN-KNMI) and the USArray Transportable Array (EarthScope).

Receiver Functions
The teleseismic coda contains different mode conversions of body waves (P-to-S wave and S-to-P wave) related to velocity discontinuities beneath the station.In this study, the focus is on P-to-S wave conversions.P-receiver functions (hereafter referred to as RFs) are isolated P-to-S wave mode conversions obtained by deconvolving the vertical component from the radial component of 3-component seismic records of teleseismic earthquakes (Langston, 1977).The deconvolution process removes the source and propagation imprints, as well as the instrument response, from the raw seismogram.It returns signals directly related to the incoming P-phase, converted P-to-S phases (Ps) and their reverberations phases (PpPs and PpSs + PsPs) from discontinuities beneath the station (Ammon et al., 1990).

Effect of Near-Surface Strong Velocity Contrast
The effects of sedimentary layers on RFs have been discussed in detail by several authors (e.g., van der Meijde et al., 2003;Yeck et al., 2013;Yu et al., 2015;Zelt & Ellis, 1999).Here we only give a succinct account of how the presence of sediment affect the RFs.The two major effects of the presence of sediment on the RFs are overprinting and time shift of the Moho phases.The three main factors to consider are: (a) the thickness of the sediment layer, (b) the magnitude of the velocity contrast between the sediment and sub-sediment layers and (c) the complexity of the sediment layer for example, anisotropy and impedance contrast within the basin.Regarding the first factor, synthetic RF studies show that a thinner sediment layer results in a broadened first arrival, which is a composite of the direct P wave and the PbS ( Figures 1a, 1b, 1e, and 1f).With increasing sediment thickness, the direct P and the PbS signals become increasingly distinguishable.However, the arrival of the PbS reverberations may coincide with Moho phases and can result in negative signals where positive Moho phases are expected (Figures 1c and 1d).Moreover, sediment layer delays the arrival of conversion phases from deeper interfaces as they pass through the low-velocity layer and the apparent delay increases with thickness of the sediment layer (e.g., Figures 1a, 1b, 1e, and 1f).In the presence of sedimentary layer, existing RF techniques such as one-crustal layer H-κ stacking, sequential H-κ stacking and resonance removal methods may fail to recover the correct crustal structure due to reversed or delayed Moho phases, very thin or very thick sediment layer, or weak or nonmonotonic sediment reverberations.Regarding the second factor, large sediment to sub-sediment velocity contrast (much slower sediment layer) results in larger amplitudes of sediment phases that may lead to stronger interference if they arrive at similar times as the Moho phases (compare Figures 1a-1d with Figures 1e-1h).Smaller velocity contrast may decrease the interference with Moho signals, although a general delay in the Moho phases is observed (Figure S1 in Supporting Information S1).Lastly, the presence of intra-basin velocity contrast introduces additional conversion phases and increases the complexity of the RFs as well as its interpretation/ inversion (Figures 1i and 1j).

Determination of Sediment and Sub-Sediment Crustal Structures
In our two-step method, we calculate separate RFs for each step.Initially, the teleseismic earthquake data is filtered using a Butterworth band-pass filter with corner frequencies between 0.05 and 2.5 Hz for the highfrequency RFs and between 0.05 and 1.25 for the low-frequency RFs.The horizontal components from both groups are rotated to radial and transverse directions and their vertical components are deconvolved from each, using a time domain iterative deconvolution (Ligorria & Ammon, 1999).In the deconvolution process, we use a Gaussian filter of 0.75 s pulse width, corresponding to Gaussian filter parameter, a = 5, for the high-frequency RFs.For the low-frequency RFs, we used a low-pass Gaussian filter with pulse width of 1.5 s (a = 1.25).In all examples, the Moho is fixed at 35 km with crustal Vp of 6.9 km/s.The Vp models are shown for all.For all layers, Vs and density are derived directly from Vp using Brocher's regression fit (Brocher, 2005).The gray lines show RFs with no sedimentary layer but only Moho.The top panel shows the effect of single sediment layer on the RFs.(a-d) We used a sediment layer of 3.0 km/s Vp and thicknesses of 0.5, 1, 2, and 3 km.(e-h) We used a sediment layer of 2.0 km/s Vp and thicknesses of 0.5, 1, 2, and 3 km.The lower panel shows the effect of two sediment layers.We used a 3-crustal layers model with Moho parameters as above and two sediment layers of (i) 2 and 3 km/s Vp at 2 and 3 km depths, respectively.(j) 2 and 4 km/s Vp at 2 and 3 km depths, respectively.

First Step: Determination of Sediment Thickness and Vp/Vs
In the first step, the average sediment thickness and Vp/Vs are determined using the H-κ stacking, which is based on a grid search of all combinations of sediment thickness (H) and Vp/Vs ratio (κ) within a defined range, using an assumed Vp of the sediment layer (for more details, see Yeck et al. (2013) and Zhu and Kanamori (2000)).We use high-frequency RFs, and although they are typically noisier, they have larger amplitudes of basin conversion phases due to the narrower pulse width and can resolve shallow structures (Yeck et al., 2013).For RFs with relatively thin sedimentary layers (up to 2 km using the assumed velocities), the Vp/Vs may not be accurately constrained due to the arrival of the sediment PbS phase almost simultaneously as its reverberated phases (Yeck et al., 2013).In this case, we propose to retrieve the Vs of the sediment layer (needed as a priori information for the second step) from available Vs estimates or directly from the Vp values using Brocher's regression fit (Brocher, 2005).

Second Step: Determination of Crustal Thickness and Vp/Vs
In the second step, we stack the low-frequency RFs to increase its signal-to-noise ratio (SNR) and obtain the average crustal structure.We use waveform-fitting of the RF data with synthetic RFs to estimate the sub-sediment structure by varying the crustal thickness and sub-sediment Vp/Vs.The synthetic RFs are calculated using the hrftn96 program (Herrmann, 2013) with model parameters including an assumed Vp value for the base of the subsediment, average ray parameter, and the thickness and velocities of the sediment layer from the first step.The densities of the sediment and sub-sediment layers used in the synthetic RFs are calculated from the assumed Vp using Brocher's regression fit (Brocher, 2005).We determine the optimal waveform fit by grid search in the crustal thickness and sub-sediment Vp/Vs parameter space.We use correlation coefficient (CC), root-meansquared error (RMSE) and phase goodness-of-fit (Ph) misfit criteria (Kristeková et al., 2009) to establish the optimal fit between the observed RF and synthetic RFs (see Equation 2 below).The Ph criteria evaluates the phase arrival time match between the observed RF and synthetic RFs.For the Ph evaluation, we normalize both the observed and synthetic RFs using Equation 1 to highlight the phase information and minimize ambiguities in the grid search.
where RF is the normalized RF; RF t is the amplitude of the RF at time t.
The overall goodness of fit (GoF) is evaluated as where ĈC is the normalized CC, R MSE the normalized RMSE, Ph the normalized Ph and wc, wr, and wp the weighting factors applied to the CC, RMSE, and Ph criteria, respectively.These factors satisfy wc + wr + wp = 1.0.The ranges of H and κ that are evaluated represent all possible combinations of Moho depth and sub-sediment Vp/Vs in the grid search.
We use the shape of the GoF function as a proxy to estimate the relative uncertainty in the best model.The relative uncertainty is determined by a 1% decrease in the GoF value of the Moho depth and sub-sediment Vp/Vs that correspond to best model fit.Hence, the uncertainties due to assumed Vp and thickness of the sediment, the assumed sub-sediment Vp and trade-off between Moho depth and sub-sediment Vp/Vs are not taken into account.

Synthetic Experiments
We use a synthetic test to illustrate the performance of our new approach in retrieving crustal structure from RFs from stations overlying low-velocity sedimentary layers.We compare our results with those from the most common standard techniques; the one-crustal layer H-κ stacking, sequential H-κ stacking and resonance filtering methods (Yeck et al., 2013;Yu et al., 2015;Zhu & Kanamori, 2000, respectively).To do this, two suites of synthetic RFs were generated (high and low-frequency) using the hrftn96 program.We added random noise of 10% of the maximum amplitude to each RF trace.We used the following model parameters: for the sedimentary layer, the Vp, Vp/Vs, and density are 3.0 km/s, 2.0 and 2.2 kg/m 3 , respectively; for the sub-sediment layer, the corresponding values are 6.9 km/s, 1.75, and 2.9 kg/m 3 ; and for the mantle, the values are 8.0 km/s, 1.78 and 3.3 kg/m 3 .We tested the above methods with sediment layers of variable thickness (0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, and 9.0 km), while other parameters were kept fixed (including the depth to the Moho, which is at 35 km in all cases; Figure S2 in Supporting Information S1).
We applied the one-crustal layer H-κ stacking, sequential H-κ stacking and the resonance filtering (combined with modified H-κ stacking) techniques to the synthetic RFs (for details of these techniques, the reader is referred to Zhu and Kanamori (2000), Yeck et al. (2013), andYu et al. (2015)).For our two-step method, we first performed H-κ stacking for the sediment layer for a depth range of 0-10 km (interval step of 0.1 km) and a Vp/Vs range of 1.65-2.5 (interval step of 0.0025).In the second step, the waveform-fitting grid search used a depth range of 20-50 km (interval step of 0.2 km) and a Vp/Vs range of 1.65-1.95(interval step of 0.0025).To make the results from all the methods comparable, the same weighting distribution of 0.6, 0.3, and 0.1 was applied the Ps, PpPs, and PpSs + PsPs phases, respectively, to all H-κ stacking analysis in the methods tested.A summary of the results is presented in Figures 2 and 3 and additional information is included in Tables S1-S4 (in Supporting Information S1).We consider derived synthetic crustal values with uncertainties of ±2.0 km and ±0.025 in the Moho depth and sub-sediment Vp/Vs, respectively, to be acceptable results.From the results of the first step of our new approach, the H-κ stacking could recover the sediment thickness and Vp/Vs in the synthetic RFs cases, besides the RF with 0.5 km sediment thickness, where the Vp/Vs is poorly constrained (Figures 2a-2c, Table S1 in Supporting Information S1).For the waveform-fitting for Moho structure, the GoF accurately recovered the input Moho depth and sub-sediment Vp/Vs in all tested cases (Figures 2d  and 2e, Table S2 in Supporting Information S1).We tested the sensitivity of the derived sediment thickness to the assumed sediment Vp and how uncertainties in the H-κ stacking derived sediment thickness and the assumed Vp of sub-sediment units affect the retrieved Moho depth in our new approach.We observed an average of 0.2 km uncertainty in sediment thickness for 0.1 km/s uncertainty in the assumed sediment Vp (Table S3 in Supporting Information S1).In the sub-sediment layer, we observed average offsets of 0.6 and 0.2 km in Moho depth for 0.1 km uncertainty in sedimentary layer and 0.1 km/s uncertainty in the assumed Vp of the sub-sediment layer, respectively (Table S3 in Supporting Information S1).
For the other methods tested, the one-crustal layer H-κ stacking method could recover the Moho depth and Vp/Vs only for the RF with 0.0 and 0.5 km thick sediment layers (Figure 3b).In the other cases with thicker sediments, one-crustal layer H-κ stacking wrongly estimates the Moho depth with offsets of 15 km to +15 km, with total variation of 30 km and offsets of 0.1 to +0.2 for Vp/Vs (Figure 3b).Similar observations were found by other studies (e.g., Yeck et al., 2013).The sequential stacking method could resolve the Moho depth (±1 km) and subsediment Vp/Vs for RFs with thin sediment thicknesses (max 1 km) (Figure 3c), while it provided inaccurate estimates for all other tests.Moreover, the sub-sediment Vp/Vs was poorly resolved for the RFs with 1 km thick sediment with 0.1 offset.In the other synthetic cases, the sequential stacking estimated the Moho wrongly with variations of 14 km and +14 km km, with a total variation of 28 km and 0.1 to +0.2 offsets for the subsediment Vp/Vs (Figure 3c).The results of the interpreted crustal structure from resonance filtering of the RFs and subsequent modified H-κ stacking shows that it wrongly estimates the Moho depth and sub-sediment Vp/ Vs when the resonance filter is poorly defined.The method wrongly estimates the Moho depth by up to 6 km and offset of 0.1 to +0.2 for sub-sediment Vp/Vs for synthetic cases with sediment thickness ≤2 km, due to nonmonotonic sediment reverberations or very thin sediment layers (Figure 3d).For cases with thicker sediment (>2 km), the method is not able to retrieve the Moho and shows offsets of 12 to +24 km with a total variation of 36 km in the Moho depths and offsets of 0.1 to +0.2 for the sub-sediment Vp/Vs.

Real-Data Examples
We used data from 5 broadband stations in the Netherlands and 1 in the USA to test the performance of our new approach, one-crustal layer H-κ stacking, sequential H-κ stacking, and resonance filtering (combined with modified H-κ stacking) (Table S5 and Figure S3 in Supporting Information S1; Zhu and Kanamori (2000); Yeck et al. (2013); Yu et al. (2015)).In the Netherlands, stations HGN and GUR1 from the NSAN-KNMI network (KNMI, 1993) and stations NE013, NE05, and NE009 from the NARS-NL network (Utrecht University-UU Netherlands, 1983) were used.In the USA, station G22A from the USArray Transportable Array located in Powder River Basin was used.For all the stations, we used records of teleseismic earthquakes with magnitude ≥5.5 and epicentral distances between 30°and 90°.The seismograms were time-windowed between 20 s before and 50 s after the theoretical P-wave arrival time, calculated using the IASP91 reference model (Kennett & Engdahl, 1991).We calculated two suites of RFs from the data as described in Section 2.3 above.We qualitycontrolled the RFs suites using three criteria: (a) least-square misfit between the original radial component signal and the product of the convolution of the radial RF with the vertical component; (b) SNR defined by the ratio of maximum amplitude in the radial RFs to the maximum amplitude in the transverse RFs; and (c) visual inspection to exclude RFs with spurious signals in both the vertical and radial components.For the highfrequency RFs suite, we used minimum RFs fit with original signal and SNR of 75% and 2.5, respectively, because they are typically noisier.Corresponding values of 90% and 5 were set for the low-frequency RFs suite.We used the Digital Geologic Model-deep (DGM-V4; van Dalfsen et al. (2006)) to constrain the Vp of the sediment layer per station in the Netherlands (Table S6 in Supporting Information S1).We used a fixed subsediment Vp of 6.9 km/s for the stations in the Netherlands, based on the shear velocity model of Yudistira et al. (2017) (Vp derived using Brocher's regression fit of Brocher (2005)).For station G22A, Vp of 3.6 and 6.7 km/s were used for the sediment and sub-sediment layer, respectively (Yeck et al., 2013).For the upper mantle properties, we used the IASP91 model (Kennett & Engdahl, 1991).In the waveform-fitting step of our new approach, we allowed for ±1 km (with 0.1 km step) freedom to the thickness of the sediment derived from the H-κ stacking to account for the uncertainties in the estimation of the Vp of the sediment layer (Section 3).Stacking weight for the sediment phases were chosen based on the clarity of the arrival amplitudes (Table S6 in Supporting Information S1).For the one-crustal layer H-κ stacking and sequential H-κ stacking, we allowed a minimum of 0.6 weight for the Moho Ps phase and details of the weighting are given in Table S8 in Supporting Information S1.For the modified H-κ stacking after resonance filtering, we used weighting of 0.5, 0.4, and 0.1 for the Moho Ps, PpPs and PpSs + PsPs phases, respectively, after Yu et al. (2015).
The results for the analysis with real-data are shown in Figures 4 and 5.Additional details are shown in Tables S6-S8 and Figures S4-S11 in Supporting Information S1.From the first step of our new approach, sediment thicknesses of 0.1 km, 0.6 km, 0.8 km, 1.4 km, and 1.9 km are found for HGN, NE013, NE05, NE009, and GUR1, respectively.These sediment thicknesses are consistent with the depth to the base of the Chalk Group (or North Sea SuperGroup) from the DGM-V4 (Figure 4a; (Doornenbal et al., 2022)).For G22A, we found sediment thickness of 2.8 km, similar to 2.9 km reported by Yeck et al. (2013).
The interpreted Moho structures from the four methods differ substantially.Crustal thickness estimates for the same station, vary at least 10 km, and can be up to 30 km difference in extreme cases.The exception is station NE013 where all methods give consistent results (Figure 4b; Figure S11 in Supporting Information S1).In general, the waveform from the solution of our new approach shows better fit with the RF data than the solutions from one-crustal layer H-κ stacking, sequential H-κ stacking, and resonance filtering methods 5b and Figure S11 in Supporting Information S1).The one-crustal layer H-κ stacking underestimates the Moho depth by 9-23 km relative to results from our new approach for stations with >0.5 km sediment thicknesses (GUR1, NE05, NE009, G22A).Similar underestimation was seen for one-crustal layer H-κ stacking in the synthetic experiment (Section 3).The derived Moho depths from sequential H-κ stacking offset the results from our new approach by 11 to +18 km for stations with >0.5 km sediment thicknesses, similar to observations from the synthetic experiments.The resonance filtering method shows similar Moho depths (±3 km) to our new approach for all the stations, except for HGN (which has very thin sedimentary layer), where an offset of 9 km was observed.In the following, we focus on the result from our new approach and link them to previously published crustal thickness estimates for these stations.Below HGN, we found Moho depth of 37 ± 0.8 km and Vp/Vs of 1.69 ± 0.14.Our Moho depth estimate beneath HGN is consistent with values from other structural highs in the Netherlands (37 km) from deep seismic reflection survey of Duin et al. (1995).Geissler et al. (2008) found a Moho depth of 32 km beneath HGN from their RF analysis using a lower average Vp of 6.5 km/s for the crust.In a study using Moho reflected phases (PmP and SmS), Sichien et al. (2012) reported a Moho depth of 31.7 ± 2.0 km from PmP and 32.7 ± 1.3 km from SmS data, using an assumed Vp velocity of 6.5 km/s for the crust.Our result for station NE013 (Moho at 38 ± 1.0 km and Vp/Vs of 1.66 ± 0.04) is consistent with a previous deep seismic reflection survey near NE013, which found the Moho discontinuity at 37 km depth (Duin et al., 1995).Similarly, our estimates for station GUR1 (Moho depth of 29 ± 1.4 km and Vp/Vs of 1.94 ± 0.07) agree with the estimates reported by Duin et al. (1995) and Yudistira et al. (2017), who found Moho depths of 28 and 29 km, respectively.Beneath NE05, we found the Moho at 32 ± 1.4 km and a Vp/Vs of 1.87 ± 0.08.Paulssen et al. (1993) found a shallower Moho (27 km) beneath NE05 using a non-linear waveform inversion method.However, these authors also reported that the lower crustal structure was not well resolved due to the dependence of the results on the starting model.We estimate Moho depth of 31 ± 1.8 km and Vp/Vs of 1.93 ± 0.02 beneath NE009.No previously published results are available for this location.For G22A, we found Moho depth at 43 ± 0.3 km and Vp/Vs of 1.68 ± 0.02, which is comparable with previous estimate of 41 km reported by Yeck et al. (2013).
The derived RFs in the real-data examples show varying degrees of sediment effects per station.For station HGN, minimal sediment interference is observed in the RFs (Figure 4c), while for the other stations, the sediment arrivals overprint the Moho phases (Figures 4d-4f, 5b and Figure S11 in Supporting Information S1).For station NE009, the Moho PmS phase was not visible in the RFs due to overprinting by negative amplitudes of sediment phases (Figure 4f).

Discussion and Conclusion
We showed how the presence of near-surface sediment layers beneath a seismic station complicates the estimation and modeling of crustal thickness and Vp/Vs using standard RF techniques.In particular, converted phases created at the sediment-basement boundary and their associated reverberations can overprint the Moho phases in the RF and, in some cases, result in negative amplitudes where positive Moho phases are expected.
Synthetic experiments show that one-crustal layer H-κ stacking, sequential H-κ stacking, and resonance filtering (combined with modified H-κ stacking) methods result in up to 30 km variations in estimates of Moho depths for the same structural model.To address these issues, we presented a new method capable of providing reliable estimates of sediment thickness, Moho depth and sub-sediment Vp/Vs structure for the crust.The method is based on a combination (i.e., two-step approach) of H-κ stacking and waveform-fitting.The above results and comparisons using both synthetic and real-data demonstrate that the proposed method can retrieve reliable estimates of Moho depth from RFs, even in the presence of low-velocity sedimentary layers of arbitrary thickness.The method is also independent of the sediment layer velocity and the strength of the reverberations generated in the sediment.However, in the presence of strong intra-basin impedance contrasts, enhanced amplitudes are expected due to the sediments and this interference, particularly at the more sensitive higher frequencies, can influence the reliability of our method in retrieving accurate Moho structure.
Additionally, high-frequency RFs needed in the first step are usually very unstable and noisy and can be affected by amplification of the ground motion due to near-surface sediment layer.Consequently, extended station deployment may be necessary in regions with substantial sediment thickness to collect enough data with high SNR for this step.As previously discussed (Section 2.3.1),Vp/Vs derived from H-κ stacking for thin sedimentary layers up to 2 km may be poorly constrained, potentially introducing errors in the results of the waveform-fitting step, due to wrong estimation of the Vs of the sediment layer used in the synthetic RFs calculation.We observed an average of 0.2 km uncertainty in Moho depth for 0.05 uncertainty in the sediment layer Vp/Vs (Table S4 in Supporting Information S1).In cases where sediment Vp/Vs is poorly constrained, the Vs can be derived from existing Vs estimates when available, or directly from the Vp velocity using Brocher's regression fit (Brocher, 2005).Moreover, the weighting applied to the sediment phases in the H-κ stacking step for the real-data examples are based on visual inspection of the phase arrivals and the clarity of these phases.
While the first step of the new approach is same as the sediment H-κ stacking of Yeck et al. (2013), the waveformfitting step for the Moho depth and sub-sediment Vp/Vs does not involve picking amplitude and arrival time (or corrected time) of the Moho phases.Using synthetic and real-data, we showed that the new method performs well, even in cases where the Moho phases are completely overprinted by large-amplitude phases (positive or negative) from the base of the sediment layer.For the synthetic experiment, we used a single RF trace and its ray parameter as input into the waveform fitting step.We observed an uncertainty of 0.4 km in the Moho depth for 0.005 uncertainty in the ray parameter.In the real-data examples, the observed RFs are stacked to obtain an average crustal structure and the average ray parameter of the traces is used in the waveform fitting step.The average ray parameter does not have significant effects on the Moho depth estimates.The weighting applied to the objective functions of the waveform-fitting step (CC, RMSE, and PG) are user-selected, making the method to be semiautomatic.A limitation of the new method is that it models only two crustal layers with sharp boundaries between them: sediment and sub-sediment.Thus, the focus of the new approach is on modeling the strongest nearsurface impedance layer and the Moho interface.Future work will address these limitations as well as the efficiency and effectiveness of the waveform modeling step.

Figure 1 .
Figure 1.Synthetic examples of the effects of sediment layer(s) on RFs.In all examples, the Moho is fixed at 35 km with crustal Vp of 6.9 km/s.The Vp models are shown for all.For all layers, Vs and density are derived directly from Vp using Brocher's regression fit(Brocher, 2005).The gray lines show RFs with no sedimentary layer but only Moho.The top panel shows the effect of single sediment layer on the RFs.(a-d) We used a sediment layer of 3.0 km/s Vp and thicknesses of 0.5, 1, 2, and 3 km.(e-h) We used a sediment layer of 2.0 km/s Vp and thicknesses of 0.5, 1, 2, and 3 km.The lower panel shows the effect of two sediment layers.We used a 3-crustal layers model with Moho parameters as above and two sediment layers of (i) 2 and 3 km/s Vp at 2 and 3 km depths, respectively.(j) 2 and 4 km/s Vp at 2 and 3 km depths, respectively.

Figure 2 .
Figure 2. (a) Interpreted sediment thickness and Vp/Vs from H-κ stacking of the synthetic RFs.H-κ stacking grid search are shown for synthetic RFs with sediment thickness of (b) 0.5 km and (c) 3 km.Waveform fitting for Moho depth and sub-sediment Vp/Vs showing the best model from Correlation Coefficient (CC), Root-Mean-Square Error (RMSE), phase Goodness-of-Fit (Ph), and Overall Goodness-of-Fit (GoF) grid search for (d) Synthetic RF with 0.5 km sediment thickness; and (e) Synthetic RF with 3 km sediment thickness.

Figure 3 .
Figure 3. Interpreted Moho depths and sub-sediment Vp/Vs using different methods: (a) Waveform-fitting grid search-new approach (b) One-crustal layer H-κ stacking; (c) Sequential H-κ stacking; and (d) Resonance filtering.The gray line in the main plots and black circle in color bars show the reference Moho depth and sub-sediment Vp/Vs in the synthetic RFs, respectively.

Figure 4 .
Figure 4.The interpreted crustal structures for real-data examples in the Netherlands.(a) Sediment layer thickness derived from H-κ stacking compared with depth to Chalk Group and the North Sea SuperGroup (DGM-V4; (Doornenbal et al., 2022)).(b) Interpreted crustal thickness and Vp/Vs from multi-methods.Forward modeled RFs of the solutions from the different methods are plotted against the stacked RF data and its double standard deviation for station (c) HGN (d) GUR1 (e) NE05, and (f) NE009.

Figure 5 .
Figure 5.The interpreted crustal structures for Station G22A.(a) Interpreted crustal thickness and Vp/Vs from multi-methods compared with results from Yeck et al. (2013) (b) Forward modeled RFs of the solutions from the different methods are plotted against the stacked RF data and its double standard deviation.