Assessing uncertainties in high-resolution, multifrequency receiver-function inversion: A comparison with borehole data

Estimating the elastic properties of the shallow crust is key in the exploration and location of earth resources such as petroleum or geothermal reservoirs. Active seismic surveys are the primary sources of information about the subsurface distribution of P-wave velocities $VP$ and about the location of discontinuities. Complementary information on S-wave velocities $VS$ can be obtained using passive seismic investigations such as receiver function (RF) analysis. RF analysis is an established technique to map first-order seismic discontinuities in the crust and upper mantle, exploiting P-to-S converted waves generated during the propagation of a teleseismic wavefield (Langston, 1977, 1979; Piana Agostinetti et al., 2008, 2009; Licciardi et al., 2014). Teleseismic P-waves that cross a seismic discontinuity are partially converted to S-waves. Due to the near-vertical incidence of teleseismic wavefronts, converted S-waves are recorded mainly on the horizontal components of a seismogram, which are usually rotated to the radial and transverse directions. The vertical (⁠$z$⁠) component is dominated by P-waves, and the general assumption made to construct a RF is that this $z$-component records instrumental, source, and path effects. Radial and transverse RFs are obtained by deconvolving the $z$-component from the horizontal components and contain P-to-S converted phases due to seismic discontinuities beneath the station. Because this deconvolution operation is unstable due to noise and the low spectral values in the $z$-component signal, it is typically stabilized by applying a low-pass filter with a characteristic cut-off frequency (Langston, 1977).

Several methods have been developed to invert an RF time series measured at a single station and infer a vertical profile of $VS$ in the crust and upper mantle beneath the station (Ammon et al., 1990; Sambridge, 1999). In recent years, RF analysis has been extended to investigations of the shallow crust (Leahy et al., 2012; Licciardi and Piana Agostinetti, 2017). Although these studies demonstrated the potential of using the propagation of teleseismic wavefields to explore the first few kilometers of the crust, the resolving power of RFs at these shallow depths is still not fully understood.

Several factors must be considered to realistically estimate the uncertainties of the inverted $VS$ profiles. First, RF inversion is a nonlinear and nonunique inverse problem (Ammon et al., 1990). Linearized inversion schemes do not provide reliable estimates of the uncertainties on the results due to the subjective choice of the damping/smoothing operator (Snoke and Sambridge, 2002). Transdimensional algorithms, developed in a full Bayesian framework, have been applied to eliminate the need for regularization and obtain more realistic estimation of the uncertainties of the inverted models (Bodin et al., 2012).

Second, uncertainties of the earth-model parameters critically rely on realistic estimates of the errors in the data. Poor assumptions on the error statistics can strongly affect the inversion results. For example, measurement errors in time-series data are often assumed to be uncorrelated white noise. Errors in RFs computed with frequency-domain estimators are formally uncorrelated in the frequency domain (e.g., Park and Levin, 2016b). However, processed seismic waveforms such as RFs are always band-limited and measurement errors are necessarily correlated in the time domain, in which the fit between observed and predicted data is evaluated. Assuming uncorrelated errors overestimates the amount of information provided by the data and leads to an underestimation of the uncertainty in the inversion. For an example of this effect in RF inversion, see Chai et al. (2017).

A third factor that affects the uncertainty of RF inversion is that RF data are widely known to be more sensitive to $VS$ jumps than to absolute velocity values (e.g., Kind et al., 1995; Schlindwein, 2006). This means that there is a trade-off between the inferred depth of a velocity discontinuity and the average $VS$ value above it. To reduce the uncertainties related to this trade-off, the ratio of the amplitude of the horizontal over the vertical component of the P-wave arrival at different cut-off frequencies has been inverted to estimate an absolute $VS$ profile (Hannemann et al. [2016], the so-called $H/Z$ ratio or P-wave polarization). Svenningsen and Jacobsen (2007) define an apparent $VS$ profile computed from RFs obtained for different cut-off frequencies of a low-pass filter used in processing teleseismic data. The simultaneous inversion of several RFs computed for different cut-off frequencies combines the benefits of a classic RF inversion, which can resolve $VS$ discontinuities, and apparent $VS$ profiles, which constrain the absolute values of $VS$⁠. Similar approaches have been adopted to invert for a multiscale seismic model in full-waveform inversions (Hong and Sen, 2009).

A fourth factor that affects the uncertainties of RF inversion is the amount of teleseismic data available. Permanent seismic stations with a decade or more of continuous recording allow for collecting hundreds of high signal-to-noise ratio (S/N) teleseismic events. Conversely, RF analysis for data recorded during short temporary deployments must rely on a limited number of teleseismic events.

In this paper, we investigate in detail the resolving power of the RF methodology by comparing inversion results with geologic and geophysical data from a deep borehole near the seismic station. We compute four stacks of RFs processed with different cut-off frequencies from four independent sets of teleseismic data, and we estimate the full covariance matrix of the data errors using data-derived autocorrelation functions. We then compare the results of two inversions. First, we invert only the stack of RFs obtained at the highest cut-off frequency (4 Hz) computed from a limited number of teleseismic events. We then use four RF stacks (for cut-off frequencies of 0.5, 1, 2, and 4 Hz) in a joint multifrequency inversion, almost doubling the total number of teleseismic events used. The RF inversion is done with a transdimensional Markov chain Monte Carlo sampling algorithm, which allows the data to determine the number of seismic interfaces in the earth model. Results from the two inversions are compared with the sonic log and the lithostratigraphy in the nearby borehole. Finally, we discuss the estimated uncertainties in the retrieved velocities and depths of seismic discontinuities as a function of the number of inverted RFs.

In this study, we reconstruct $VS$ and $VP/VS$ profiles beneath a seismic station and compare the results with geologic and geophysical data acquired in a nearby borehole. The PUGLIA-1 deep borehole was drilled to 7070 m below the surface in the 1980s, and sonic log and lithostratigraphy data are now open access. The drillsite is on the well-studied Apulian carbonate platform in Southern Italy. A permanent seismic station (MRVN, located approximately 900 m away from the PUGLIA-1 borehole) has been deployed in the past few decades, providing a large amount of teleseismic records.

The Apulian platform is a multilayer carbonate platform that developed in the Mesozoic-early Cenozoic on the southern margin of the Tethyan Ocean (Channell et al., 1979). In the last 15 Ma, the Apulian platform has been partially involved in the Apenninic compressive tectonics. The PUGLIA-1 borehole was drilled in the undeformed portion of the Apulian platform, which is the foreland of the Apenninic mountain belt (Figure ¹).

The lithostratigraphy of the PUGLIA-1 borehole consists of five main units (Figure ¹). On top, a wackestone carbonate formation (Calcari di Cupello) outcrops at the surface and is found down to 890 m depth. Below, a limestone formation with different degrees of dolomitization and some mudstone layers can be followed to 3535 m depth. Below this depth, the borehole was drilled in a dolomitized formation and the lithostratigraphy lacks the resolution found at shallower depth. The evaporites of the Burano Formation (anhydrites) are found at 5065 m depth. Finally, at 6112 m depth, the borehole reached a red argillite formation above a thick sandstone layer. The borehole ended at 7070 m in sandstones cemented by silica and carbonate.

Sonic-log measurements were acquired between 910 and 7070 m depth (Figure ¹). The higher values of $VP$ in the sonic log are the more reliable measurements of the in situ velocity; low values are likely due to borehole irregularities. The average $VP$ measured by the sonic log is approximately $6km/s$ in the top 3.5 km of the borehole. Then, $VP$ increases to approximately $6.5km/s$ and remains constant down to 6 km depth. Below this depth, $VP$ drops to approximately $5km/s$ in the siliciclastic formation and remains approximately constant down to the total depth.

The regional seismic structure in the crust and upper mantle in the southern Italian peninsula has been investigated using RF analysis in several previous publications (Steckler et al., 2008; Piana Agostinetti and Amato, 2009; Miller and Piana Agostinetti, 2011). Beneath the Apulian platform, these studies found a very smooth topography of the crust-mantle boundary with an average Moho depth of approximately 30 km. In recent work focused on the entire Apulian foreland, Amato et al. (2014) find variations in the total thickness of the carbonate platform (4–12 km thick) and in the thickness of single layers (as also noticed in Improta et al., 2000). Also, a 3–4 km thick siliciclastic layer at the base of the Apulian carbonate platform is present at the northern end, but not at the southern end. Given the low vertical resolution of their RF data set, Amato et al. (2014) could only make a qualitative comparison of their $VS$ profile with the sonic log of the PUGLIA-1 borehole.

To compute a RF data set, we adopt the approach developed by Di Bona (1998), which is based on the frequency-domain deconvolution scheme of Oldenburg (1981) applied to RFs by Ammon (1992). To stabilize the ill-posed deconvolution, this method uses a time-domain Gaussian averaging filter that acts as a low-pass filter with a cut-off frequency that depends on the width of the filter. Frequency-domain deconvolution can give stable results for high-frequency cutoffs, e.g., using multiple-taper methods (Park and Levin, 2000, 2016b). The scheme of Di Bona (1998) also estimates a variance of the error associated to the RF derived from each event that takes into account the background noise and signal-generated noise. We use this variance to weight the stacking (see below).

We computed RFs from teleseismic events recorded at seismic station MRVN, which belongs to the Italian National Seismic Network (network code IV) and has been deployed in 2007. We first selected few hundreds of high S/N events after visual inspection of the P-wave arrivals. In a successive step, RFs were computed excluding noisy seismograms and traces that were contaminated by ringing effects. Given the vast amount of data, the back-azimuthal coverage is nearly complete. We first show the RF data as a function of the back azimuth to test the assumption of a horizontally stratified earth model. In Figure ², we plot the RF data set computed using all available events (approximately 100 RF with high S/N) using a Gaussian filter that cuts-off frequencies greater than approximately 2 Hz. The RFs in Figure ² were binned in $20°×40°$ back-azimuthal and epicentral distance bins. The bins use a 50% overlapping scheme, so that each RF contributes to two bins. Blue (positive) and red (negative) amplitudes correspond to converted P-to-S phases, and the radial and transverse RFs are shown. Positive pulses on the radial RF are generated by velocity contrasts, where $VS$ increases with increasing depth. The presence of only a few low-amplitude pulses on the transverse RFs suggests a horizontally layered structure beneath this seismic station.

To further investigate this point, we also computed the angular harmonics of the RF data set (e.g., Piana Agostinetti and Miller, 2014). The first harmonic (⁠$k=0$⁠) is the RF signal that shows no variation with back azimuth; the second harmonic (⁠$k=1$⁠) quantifies the variation of the RFs with a $2π$ periodicity in back azimuth, and contains energy converted at dipping discontinuities or at the top or bottom of anisotropic layers with a dipping symmetry axis; the third harmonic (⁠$k=2$⁠) quantifies the variation of the RFs with a $π$ periodicity in back azimuth, and it contains energy converted at the top or bottom of anisotropic layers with an horizontal axis of symmetry. In principle, the $k>0$ harmonics can be used to obtain additional information about dipping interfaces and anisotropy beneath the seismic station (Park and Levin, 2016a). The $k=0,1,2$ harmonics in our data set are shown in Figure ². Although the amplitudes for $k=1$ and 2 are not negligible, they are much smaller than those in the $k=0$ trace. In this paper, we neglect the possible effects of anisotropy and dip and concentrate on the 1D isotropic, horizontally stratified structure beneath the seismic station.

The $k=0$ harmonics were calculated from the full data set of Figure ². In our inversion, we assembled four subsets from the full data set and computed four RFs for different cut-off frequencies by stacking (rather than taking the $k=0$ harmonic). We only used RFs from events within a small range in epicentral distance (70°–90°, ray parameter range $Δp=0.013s/km$⁠) to avoid having converted waves generated at the same interface having different time delays. To obtain stacks of RFs at different cut-off frequencies with statistically independent errors, we subdivided the selected 40 teleseismic events in four nonoverlapping subsets and computed RFs for cut-off frequencies of 0.5, 1, 2, and 4 Hz. Given the evidence for a horizontally stratified structure beneath the station, we expect to have statistically consistent RFs at each frequency, even if the back-azimuthal coverage is less complete compared with the full RF data set in Figure ². A list of all the teleseisms used for each cut-off frequency is shown in Table S1 (Supplementary information can be accessed through the following link: S1.).

In Figure ³, we show the RFs computed for the cut-off frequency of 4 Hz and the resulting stacked RF. Stacking is performed by weighting each RF with the inverse of its noise variance, estimated with the method of Di Bona (1998). The stack residuals (differences between each RF and the stacked value) allow for characterizing the errors in the stacked RF. We use the stack residuals and the stacking weights to compute a time-dependent standard deviation of the weighted average in the RF stack. We also calculate the average autocorrelation of the stack residuals (again, applying the stacking weights). In summary, for each frequency $f$⁠, we obtain from the stacking procedure three vectors: a stacked RF $d[f](t)$⁠, a time-dependent standard deviation of the errors $σ[f](t)$⁠, and an autocorrelation function of the errors $r[f](t)$⁠. These vectors are illustrated in Figure ⁴.

Following this approach, we obtain the covariance matrix of the RF data errors from a straightforward analysis of the stacking process, ensuring a realistic characterization of the data errors. We do not need to assume that the data errors are white noise or to add to the inversion unknown parameters of the data-error autocorrelation (Malinverno and Briggs, 2004) or error variances (Malinverno and Parker, 2006). To better illustrate this point, in the following, we also report results of an inversion, in which such parameters are used. Finally, in contrast to other applications of similar algorithms (e.g., Bodin et al., 2012), in our implementation, the standard deviation of the data errors is not constant with time, so that we can account for different variabilities in the stacked RFs at different times (Figure ⁴).

The RF inverse problem, i.e., the reconstruction of 1D profiles of $VS$ and $VP/VS$ ratio from RF data is nonlinear and nonunique. Several solution strategies have been implemented from linearized inversion (Ammon et al., 1990) to Monte Carlo sampling (Sambridge, 1999). In this study, we apply the algorithm developed by Piana Agostinetti and Malinverno (2010), which used transdimensional Markov chain Monte Carlo sampling. The algorithm is based on a Bayesian framework, in which the objective of the inversion is a posterior probability distribution (PPD) of the earth-model parameters. Our earth model is a 1D sequence of layers with constant values of $VS$ and $VP/VS$ in each layer. A key advantage of a transdimensional approach is that the number of parameters (in our case, the number of layers in the earth model) is not fixed a priori but is determined by the data. In this study, the investigated depth range extends from the free surface to 60 km depth. Although we are mainly interested in the seismic structure of the shallow crust (0–10 km depth), we model the lithospheric structure down to 60 km depth to avoid misinterpretations of the arrivals. In fact, we expect free-surface multiples from shallow crustal interfaces to arrive within the same time window as P-S converted phases from deeper interfaces (e.g., the Moho). Prior constraints for $VS$ and $VP/VS$ are the same as those described in Amato et al. (2014), except for the near-surface $VS$⁠. To account for surface-geology information, the prior probability distribution of $VS$ in the 0–250 m depth range has a prior mean of $2km/s$ and a prior standard deviation of $0.1km/s$⁠.

As an example of a standard RF inversion for the shallow crust, we first invert the stacked RF computed with a cut-off frequency of 4 Hz (Figure ⁴), which is the stack of 27 RFs with a high S/N. A 4 Hz cut-off frequency has been widely used to image the fine-scale structure of the crust (e.g., Leahy and Park, 2005). We run 100 independent sampling chains and sample $6×105$ earth models in each random walk. We discard the first $105$ models sampled in each chain (the so-called burn-in phase). Because only one earth-model parameter is modified at each step, the sampled models make up a highly correlated sequence. To estimate properties of the PPD, we use one model out of each thousand, so that the final sample consists of $5×104$ earth models. The computation time is approximately 2 h per chain on a 100 CPU cluster where each chain is run independently on 1 CPU.

In Figure ⁵, we illustrate the full PPD as posterior marginal distributions of interface depths, $VS$⁠, and $VP/VS$ as a function of depth. This is a first-order, conservative representation of the PPD that does not account for posterior correlations between earth-model parameters (see the “Discussion” section). In the distribution of the depths to interfaces (Figure ⁵), we recognize four modes at 0.2, 1.0, 2.8, and 5.5 km depth. The PPD for $VS$ (Figure ⁵) shows relatively well-defined values, whereas it is quite broad for the $VP/VS$ ratio (Figure ⁵). The PPD of $VS$ confirms the presence of velocity discontinuities at the depths of model interfaces (Figure ⁵).

The posterior mean value of $VS$ is approximately $2km/s$ near the surface and increases to $3km/s$ at 1 km depth. The posterior mean $VS$ remains nearly constant down to 3 km depth, where it increases to a value of $3.2km/s$⁠. Approximately 5.5 km depth, we observe a clear decrease of $VS$ to approximately $2.7km/s$⁠. Below this discontinuity, $VS$ increases smoothly to approximately $3.4km/s$ at 10 km depth. We measure posterior uncertainties as the half-width of the interval containing the central 90% of the sampled values. $VS$ uncertainties range between $±0.5$ and $±1km/s$ along the 1–10 km depth interval. The PPD for the $VP/VS$ ratio is less informative, and a constant $VP/VS$ ratio between 1.7 and 1.85 would be within the posterior uncertainties. The depths of interfaces in Figure ⁵ can be directly related to prominent features in the RF (Figure ³). For example, the 5.5 km depth interface, a negative velocity jump where $VS$ decreases with depth, corresponds to the negative arrival at approximately 0.8 s. The positive velocity jump at 2.8 km depth is mapped as a positive arrival at approximately 0.5 s. Finally, the interface at 1 km depth is responsible for the broadening of the direct P pulse at approximately 0.2 s because the direct P pulse and the P-S phase generated at this shallow depth cannot be fully separated at 4 Hz.

The PPD of the number of interfaces in the sampled earth models (Figure ⁵) shows a broad distribution with a mode at approximately 40 interfaces. It is important to note that the algorithm uses several closely spaced interfaces to reproduce gradational changes of $VS$ with depth. Figure ⁵ shows about five major layers in the top 10 km. Looking at the same distribution as in Figure ⁵ for the full depth range modeled (0–60 km, not shown), there are about eight major layers in the sampled models. The intrinsic parsimony of the RjMcMC algorithm ensures that the sampled models only contain interfaces, and thus layers that are needed to fit the data (e.g., Malinverno, 2002). Finally, Figure ⁶ displays the distribution of the RF data predicted by the sampled earth models. The predicted values generally fall within the $±2σ$ uncertainty envelope of the observed RF stack.

In our second inversion, we use as input data the four stacked RFs computed for different cut-off frequencies (Figure ⁴). We follow the same procedure used in the first inversion, obtaining a final sample of $5×104$ earth models to evaluate the PPD. The computation time is again approximately 2 h per chain on a 400 CPU cluster (here, we use 4 CPUs for each Markov chain, one for each frequency).

The PPD of the depths to seismic interfaces (Figure ⁷) shows four clear modes at approximately 0.2, 1.0, 3.2, and 5.8 km depth. These modes are slightly better defined over the background than in the single-frequency inversion (Figure ⁵). The posterior mean of $VS$ (Figure ⁷) shows well-defined velocity discontinuities corresponding to the modes in Figure ⁷. The posterior mean $VS$ starts at approximately $2.5km/s$ near the surface and increases to $3km/s$ approximately 1 km depth. $VS$ retains this value down to 3.2 km depth, where it increases to $3.6km/s$⁠. Between 3.2 and 5.8 km depth, $VS$ increases slightly and then drops to approximately $2.9km/s$⁠, maintaining this value down to 7.5 km depth. Below this depth, $VS$ smoothly increases again to approximately $3.3km/s$ at 10 km depth. The posterior uncertainties in $VS$ are approximately $±0.3km/s$ down to 6 km depth, and then increase to approximately $±0.5km/s$⁠. The $VS$ uncertainties in the multifrequency inversion are clearly less than those obtained with a single-frequency RF (Figure ⁵).

As in the single-frequency inversion, the PPD of the $VP/VS$ ratio has a relatively large uncertainty (Figure ⁷). In the multifrequency inversion, however, the variation with depth of $VP/VS$ is more pronounced, being approximately 1.8 in the 1–6 km depth interval and then decreasing to 1.7 between 6 and 9 km depth. Comparing Figures ⁵ and ⁷, we note that the PPD of the number of layers is similar, suggesting that the number of layers is controlled by the RF with the highest frequency content.

Figure ⁸ displays the distribution of the multifrequency RF data predicted by the sampled earth models. As in the single-frequency inversion, the predicted values generally fall within the $±2σ$ uncertainty envelope of the observed RF stacks. The only exception is the direct pulse (0 s time delay) for the 0.5 Hz RF. In this case, the mean predicted RF follows the $−2σ$ lower boundary.

As expected, the results from the two inversions for a single RF and four RFs at different cut-off frequencies are quite similar in many aspects, but at the same time, they present interesting differences. We use these differences and the comparison of the retrieved $VS$ and $VP/VS$ profiles with the PUGLIA-1 borehole data to evaluate the resolving power of the RF methodology and discuss its potential as a tool to explore the shallow crust.

Figure ⁹ compares the PPD for the depth of the interfaces in the two inversions with the interpreted lithostratigraphy in the PUGLIA-1 borehole. Both inversions show four main discontinuities in the borehole interval, confirming that the maximum frequency content (4 Hz in both cases) controls the degree of resolution in depth. Comparing the position of these main discontinuities with the lithostratigraphy, we interpret them as the contact between wackestone and limestones (at approximately 1 km depth), the contact between limestones and dolomites (approximately 3 km depth) and the contact between evaporites and siliciclastics (5.5–5.8 km depth). This comparison indicates that the RF methodology can recover the main lithologic discontinuities, but finer scale layering cannot be resolved because different layers are too thin and have relatively small differences in elastic parameters.

The depth of the first interface at approximately 1 km is about the same in both inversions. Our results indicate that there is no substantial difference in using one or more frequencies (and also doubling the number of teleseismic events) to resolve the very shallow structure of the crust. We argue that the relatively small uncertainty in the inverted $VS$ in the top kilometer or so results in a relatively small error in interface depths. Our results have been obtained for a station on a thick carbonate sedimentary sequence. Similar tests should be carried out with stations deployed in different geologic setting (e.g., over thick unconsolidated sequences in a clastic sedimentary basin) to extend these considerations.

In the deeper portion of the profile, the multifrequency inversion gives an estimated depth of the velocity inversion between evaporites and siliciclastics of 5.8 km, which is approximately 300 m shallower than the depth observed in the borehole (6.1 km, the “true” value hereafter). We observe that the histogram of the depth of the interfaces for the multifrequency inversion is quite broad and clearly encompasses the true depth of the discontinuity (Figure ⁹). On the other hand, the estimated depth for this interface from the single RF inversion is approximately 600 m shallower than the true value. In this case, inverting a single 4 Hz RF could bias substantially the retrieved depth of the interface, indicating that uncertainties in $VS$ decrease the accuracy of the inversion at this depth.

In Figure ¹⁰, we compare $VP$ profiles reconstructed from the two inversions with $VP$ measurements from the PUGLIA-1 borehole sonic log. The reconstructed $VP$ profiles are computed as the product of the posterior mean $VS$ and the posterior mean $VP/VS$⁠. The uncertainties in the reconstructed $VP$ are computed from the posterior uncertainty in $VS$ times the posterior mean $VP/VS$⁠. Although this definition of the uncertainties on the reconstructed $VP$ profile is not entirely consistent (it does not account for uncertainties in $VP/VS$ and for possible correlations), it gives a general idea of the relative variation of $VP$ uncertainties with depth.

The results of the two inversions broadly agree with the sonic log, showing two intervals of increasing velocity with depth followed by a velocity inversion approximately 5.5–6 km. The reconstructed $VP$ values are generally lower in the inversion of a single RF. This translates into a shallower estimate of the depth of the interfaces, as shown in Figure ⁹. In particular, in the 3.5–6 km depth range where dolomites and evaporites are found, the $VP$ values retrieved from the multifrequency inversion are closer to the values measured in the sonic log. These higher $VP$ values are also more consistent with laboratory measurements on evaporitic samples from the Apennines (approximately $7km/s$ for dolostones and approximately $6.4km/s$ for anhydrites; Trippetta et al., 2010). Both inversions give similar magnitudes of the two main velocity contrasts, namely the wackestone/limestone contact and the velocity inversion at the evaporites/siliciclastic contact. Both inversions estimate a $ΔVP=+1.8$ and $−1.5km/s$ for the shallower and the deeper interfaces respectively, confirming that RF analysis of teleseismic data can constrain the polarity and size of velocity jumps beneath a seismic station.

It is important to stress that the posterior uncertainties displayed in Figures ⁵ and ⁷ are marginal distributions of the earth-model parameters. This is a conservative estimate because ignores posterior correlations between parameters, which clearly exist in RF inversion results. In Figure ¹¹, we illustrate this point with the posterior marginal distribution of $VS$ near the Moho. In this case, the depth of the Moho interface and $VS$ are correlated in that a greater interface depth can be compensated by an increase in $VS$ above (see, for example, the formulas in Zhu and Kanamori, 2000). A simple display of the posterior marginal distribution of $VS$ does not show this correlation and could be misinterpreted. For example, looking at Figure ¹¹, a constant $3.7km/s$ velocity profile between 22 and 32 km depth falls within the posterior uncertainties, apparently suggesting that the inversion does not require a velocity discontinuity at the Moho. However, this is not the case because the sampled values of $VS$ above and below the Moho interface are correlated. Displaying posterior uncertainties with marginal distributions is clearly a first-order approximation; posterior correlations between relevant parameters (e.g., $VS$ at two different depths) can be always computed from the sampled earth models.

Another advantage of the approach used here is the possibility of including additional prior information on the velocity variation with depth, e.g., using $VP$ values from a sonic log in a nearby borehole. Constraining the prior PPD with this additional information could be very useful for exploration geophysics. In this study, we presented a more general case that does not include extra constraints because they are not always available. The posterior uncertainties could be considerably reduced by considering correlations between parameters and by imposing additional prior constraints.

The multifrequency inversion slightly increases the resolution of the $VP/VS$ ratio. Comparing Figures ⁵ and ⁷, the difference in the $VP/VS$ ratio between evaporite and siliciclastic layers is clearly more pronounced in the multifrequency inversion. In both inversions, the upper limestone and evaporite layers displays a $VP/VS=1.8−1.85$⁠, as expected for these rock types (e.g., Trippetta et al., 2010). The $VP/VS$ in the siliciclastic layer at the base of the carbonate platform is constrained between 1.7 and 1.75, as expected for continental clastic deposits (e.g., Christensen, 1996). It is worth noting that the discontinuity in the $VP/VS$ ratio retrieved from the multifrequency inversion results in a $ΔVP$ at the interface that is closer to the value measured in the sonic log. The inverted $ΔVS$ at this interface is approximately $0.5km/s$⁠. Thus, by adopting an average $VP/VS$ for the continental crust (1.732; Christensen and Mooney, 1995), we should have a $ΔVP$ drop equal to $ΔVP=ΔVS·1.732=0.86km/s$⁠, which is less than half the $ΔVP=1.8km/s$ drop measured in the sonic log. Conversely, using the $VP/VS$ ratio in the mean posterior profile, we obtain a predicted velocity drop of approximately $ΔVP=1.5km/s$⁠.

As illustrated above, we inspected the data uncertainties and error correlation to obtain realistic uncertainties on model parameters. However, we recognize that our detailed analysis on the error statistics is based on the large amount of data available. A hierarchical Bayes approach can be adopted to infer information on the absolute value of the data uncertainties, if such information is missing (Malinverno and Briggs, 2004). Moreover, hierarchical Bayes can also be adopted in case, in which assumptions on the forward modeling (e.g., 1D flat layering) could be violated (Bodin et al., 2012). In such cases, a factor $σRF$ that multiplies the standard deviation of the RF errors can be used as a “hyperparameter,” and its PPD can be retrieved at the end of the McMC sampling. This hyperparameter can account for “that part of the data that we do not expect the model to explain” (Scales and Sneider, 1998). We test our estimates of the error statistics by inverting the single 4 Hz RF stack within a hierarchical Bayes scheme. In this way, we directly derive from the data the amplitude of the errors expressed as a factor $σRF$ that scales all the elements in the covariance matrix of the data errors. Figure ¹² shows the results of this test. In our case, the posterior values of the hyperparameter $σRF$ the close to 1.0, indicating that our assumptions about error statistics and 1D horizontal layering are generally correct.

In our study, we set the maximum cut-off frequency to 4 Hz. Our results show that this frequency content allows for a resolution of approximately 0.35 km for the depth of a seismic interface within the 0–5 km depth range. However, recent studies indicate that teleseismic waveforms contain information on the shallow crust at even higher frequencies (e.g., Licciardi and Piana Agostinetti, 2017) to a maximum of approximately 15 Hz (Leahy et al., 2012). In such cases, the vertical resolution of the RF methodology could be even higher, especially in the first few kilometers below the surface. Uncertainties on the interfaces depths also depend on whether the seismic discontinuity is a sharp or a gradational boundary. We followed here a conservative approach and did not consider the frequency band between 4 and 15 Hz, but investigating higher frequencies could be useful to recover the details of seismic interfaces.

Given the large teleseismic database available, we computed four stacked RFs for different cut-off frequencies from four completely independent sets of teleseisms. This allows us to estimate the covariance matrix of the data errors, needed for retrieving realistic uncertainties of the earth-model parameters. If a small number of teleseisms were available, it is possible to compute several stacked RFs for different cut-off frequencies using the same set of teleseisms. This approach could still provide information on the crustal structure given the different frequency content in the source function of different teleseisms, the different dispersion along the source-receiver path for different frequency bands, and the frequency-dependent response of the instrumentation. In this case, however, the errors in the RF data obtained for different cut-off frequencies cannot be assumed to be uncorrelated as is done here. Assuming uncorrelated data errors when they are correlated will result in an underestimate of the uncertainty in the inversion.

Considering an average of about one teleseismic event per day (with magnitude $MW≥5.0$ in the 30°–100° epicentral distance range), a period of three months to a year of continuous recording would be necessary for an RF analysis similar to what is reported in this study. Multichannel (multistation) RF analysis can decrease the amount of data required (Mostafanejad and Langston, 2017). From the point of view of geophysical exploration, the duration of the deployment implies that instrumentation for passive seismic measurements should be installed ahead of active seismic surveys. Moreover, due to the different frequencies in the signals involved and the low probability of simultaneous recording of an explosion and a teleseismic wave, active seismic operations, i.e., explosions, do not interfere with passive seismic data acquisition and the two surveys can be planned to overlap.

We investigated uncertainties in $VS$ and $VP/VS$ in the shallow crust retrieved via RF inversion. We performed two inversions, one using a single RF computed with a high cut-off frequency (4 Hz) from 27 teleseismic events and another using four RFs computed for different cut-off frequencies (0.5, 1, 2, and 4 Hz) from a total of 50 events. Our goal is to test how inversion uncertainties depend on the number of teleseismic events and on the frequency content of the RFs. The 1D velocity profiles obtained with these RF inversions have been compared to independent measurements, namely, the lithostratigraphy and the sonic log from a borehole near the seismic station.

We quantified uncertainties in the results from the spread of PPDs defined in a Bayesian formulation of the inverse problem. These inversion uncertainties depend critically on a realistic quantification of the errors in the data. To this end, we estimated the covariance matrix of the errors in the RF data from an analysis of the residuals in the stacking of RFs obtained for different teleseismic events.

We found that the multifrequency RF inversion results match more closely the vertical variation of elastic properties compared to the inversion of a single RF and that it can be a useful approach for investigating the shallow crust in an interval of interest to exploration geophysicists. In the multifrequency inversion results, the posterior uncertainties of $VS$ are smaller (⁠$±0.3$ to $±0.5km/s$ versus $±0.5$ to $±1km/s$⁠), the reconstructed depth to a velocity inversion approximately 6 km below the surface is more accurate, and the $VP$ sonic-log measurements are matched more closely. We estimate that a minimum period of approximately 3–12 months is necessary to acquire enough teleseismic data to obtain RF inversion results similar to those presented here.

N. Piana Agostinetti’s research is funded by the Austrian Science Fund (FWF): M2218-N29. This publication resulted from research partially conducted with the financial support of Science Foundation Ireland under grant numbers 11/SIRG/E2174 and 14/IFB/2742. N. Piana Agostinetti wishes to acknowledge D. Melini at INGV-Rome for the provision of computational facilities and support. The original well log for PUGLIA-1 can be downloaded from the VIDEPI website (http://unmig.sviluppoeconomico.gov.it/videpi/videpi.asp) and the seismic data for station MRVN from the EIDA website (http://www.orfeus-eu.org/eida/eida.html). Figures were prepared using the Generic Mapping Tools of Wessel and Smith (1998).