Resolution analysis of joint inversion of seismic receiver function and surface wave dispersion curves in the “13 BB Star” experiment

Joint inversion of Rayleigh wave phase velocity dispersion and P receiver function has been applied to study the structure of the upper mantle beneath the south-western margin of the East European Craton. The data were gathered in the passive seismic experiment “13 BB Star” (2013–2016) in the area of the crust recognized from previous borehole and refraction surveys. Several fundamental issues inherent in the linearised inversion were addressed in this work, including exploitation of a priori knowledge, choice of model’s depth, trapping by local minima associated with non-uniqueness of the 5 misfit-function optimization problem, proper weighting of data sets characterized by different uncertainties, and credibility of the final models. The last was investigated with novel 1D checkerboard tests juxtaposed with resolution matrix analysis. We advocate the usefulness of linearised approach when handled with proper care, and show that the resolution analysis is an indispensable step when choosing the inversion parameters. It allowed us to obtain reliable S-wave velocity models down to 200 km depth beneath the “13 BB Star” array, indicating the presence of a Paleozoic asthenosphere and the ceiling of the 10 deeper, Precambrian, lithosphere-asthenosphere transition zone.


Introduction
Limited resolvable depth and sensitivity are the two major problems inherent in receiver function (RF) and surface wave dispersion (SWD) analysis.The first issue results mainly from the presence of noise and imperfection of the measurements.
The second constitutes the intrinsic feature of the data.Combined, they make the credible models of subsurface difficult to obtain.Although it is the inversion of RF which is usually regarded as exceptionally non-unique (Ammon and Randall, 1990;Bodin et al., 2014), the SWD lacks uniqueness either, not being able to discriminate between the fine structures, especially at greater depths (Tsuboi and Saito, 1983;Romanowicz, 2002).
The simultaneous inversion of RF and SWD, having been applied to study deep lithosphere for 20 years (Özalaybey et al., 1997;Du and Foulger, 1999;Julià et al., 2000), mitigates these issues thanks to the complementary sensitivity of the data (Shen et al., 2013).Even in this case however, there remains a certain ambiguity, showed in the following sections, that should be taken into careful consideration throughout inversion.A single model obtained by minimization of the data-misfit function cannot be regarded as a full solution of the inverse problem (Gubbins, 2004), which is often a case of linearised inversion (Julià et al., 2003;Horspool et al., 2006;Wang et al., 2014;Sosa et al., 2014;Bao et al., 2015;Li et al., 2016).One needs to probe 1 Solid Earth Discuss., https://doi.org/10.5194/se-2017-58Manuscript under review for journal Solid Earth Discussion started: 14 July 2017 c Author(s) 2017.CC BY 4.0 License. the whole space of physically plausible solutions, by drawing inferences from an ensemble (Sambridge and Mosegaard, 2002).
Bayesian Monte Carlo methods (Green and Hastie, 2009;Bodin et al., 2012;Shen et al., 2013;Deng et al., 2015;Fontaine et al., 2015) appear to be a natural approach to achieve this goal, not only performing importance sampling (Sambridge, 1999), but also properly exploiting a priori knowledge of the problem (Malinverno, 2002).They are not flawless though, but suffer from their own inherent nuisance, e.g. the lack of objective criterion of convergence (Roberts et al., 1996(Roberts et al., , 1997)), not to mention the computational cost.They also require nontrivial tuning to prevent solutions from following the shape of prior velocity distribution other than homogeneous (Minato et al., 2008;Wathelet, 2008).
Here we propose the compromise between the simplicity and transparency of the method on the one hand, and a full solution of the inverse problem on the other, by performing linearised inversion with ensemble of starting models covering entire space of acceptable solutions.We demonstrate this approach with teleseismic data collected in the area of well-recognized crust -the a priori knowledge which we introduce into inversion.We show a remnant non-uniqueness of the joint RF and SWD inversion, which justifies the careful selection of the values of initial parameters.To this end, we perform the resolution analysis based on novel 1D checkerboard tests and resolution matrices before the actual inversion.

Data
Passive seismic experiment "13 BB Star" was dedicated to study the deep structure of the Earth's interior in the marginal zone of the East European Craton (EEC) in northern Poland (Grad et al., 2015).The seismic network consisted of 13 broadband stations on the area of ca.120 km in diameter (Fig. 1).The network was located in the area of well-known sedimentary cover and crustal structure (Grad et al., 2009;Polkowski and Grad, 2015;Grad et al., 2016).The "13 BB Star" seismic stations were equipped with Reftek 151-120 Observer seismometers and Reftek 130 data logger.The stations were operated from June 2013 to October 2016.The distance and azimuthal epicentral distribution of analysed earthquakes are shown in Fig. 2.

P receiver function
The receiver functions techniques (Langston, 1977;Vinnik, 1977) have been used to investigate the structure of the lithosphereasthenosphere system.Receiver functions were calculated by slightly modified method presented by Wilde-Piórko (2015); Wilde-Piórko et al. (2017).Seismograms used in receiver function analysis has been selected manually.Additionally, the second manual selection was done after the calculation of receiver functions to choose the traces with the highest signal to noise ratio.Ultimately, the total number of 99 events was taken into consideration.Seismograms were cut 300 s before and 300 s after the theoretical P-onset calculated for the iasp91 model (Kennett and Engdahl, 1991) for earthquakes with epicentral distance of 30-100 • and viewed to leave for further analysis only that ones with visible energy on vertical components.Onepass low-pass filtering with Butterworth filter of corner frequency 5 Hz was applied before resampling seismograms to 20 Hz.
Calculation of backazimuth and polarization angle were performed in two steps: (1) seismograms were filtered with two-pass band-pass Butterworth filter with corner periods of 2 and 10 s for searching the backazimuth angle for teleseismic P-wave; N and E components of seismograms were rotated with backazimuth angle from 0 to 360 • every 1 • and radial receiver function (RFR) were calculated by time-domain Wiener deconvolution; next RFR with maximal sum of amplitudes between time 0 and 1 s were selected (equivalent of rotation ZN E components to ZRT); (2) seismograms were filtered with two-pass band-pass Butterworth filter with corner periods of 2 and 10 s for searching the polarization angle for teleseismic P-wave; seismograms were rotated from ZNE to ZRT components by backazimuth angle found in step (1) and then rotated from ZRT to LQT with polarization angle from 0 to 45 • every 1 • and Q receiver function (RFQ) were calculated; next RFQ were cut 5 s before and 5 s after direct P-wave, a mean value and trend were removed and root square mean (rms) values were calculated for each RFQ for time window between -2 and 0 s; then as best were chosen RFQ for which the rms value was lower than rms calculated for RFQ of next polarization angle (equivalent of rotation ZRT components to LQT).The final RF were filtered with Gaussian filter with parameter 4.
An example of presented procedure is shown in Fig. 3 together with RFs calculated in traditional way -ZNE components were rotated to LQT ones for theoretical backazimuth and polarization angles.The final RFR used in linearised joint inversion were move-out corrected for mean slowness of each station and stacked (Fig. 4).

Fundamental mode Rayleigh wave phase velocity dispersion
To calculate dispersion curves from surface waves the ASWMS package was used (Jin and Gaherty, 2015).It allows measuring phase and amplitude of surface waves from raw seismic waveforms and use these measurements to calculate phase velocity maps for different wave periods using Eikonal and Helmholtz equation.
The ASWMS package bases calculation on surface waves incoming to array of stations form natural teleseismic events.
Before performing phase velocity calculation data obtained during "13 BB Star" experiment required preparation.Firstly, list of all recorded events was created and manually checked for quality.Total of 706 events from period from 2013-07-10 to 2016-09-23 were selected for further study.All 706 events were manually verified on station basis.For each event list of stations with good quality recording was prepared.Stations with no data, partial data or bad data quality were omitted.135 events were selected as good quality on all 13 stations, 478 events were recorded with good quality on over 50 % of stations 10 and 103 events rejected on all stations.From all 9178 station-event pairs 5697 were selected as good quality (over 62 %).In average 8 stations were selected as good quality per event.For each station-event pair raw seismic data was converted form daily mini-seed files to single SAC files with embedded event parameters (latitude, longitude, depth and magnitude).SAC data    phase velocity maps for different wave periods.In this case maps were calculated for periods from 10 to 300 s every 5 s.Results were then prepared as dispersion curves, with values taken from corresponding grid cells on all maps.Final dispersion curve was obtained by averaging all dispersion curves for each period and is shown as green dots and green line in Fig. 5, with ±3σ deviation.

5
We used a linearised damped least-squares inversion scheme implemented in CPS package (Julià et al., 2000;Herrmann, 2013), appropriate for over-determined problems such as inversion of RF.The main step of the procedure is to find a general inverse of the forward operator matrix (forward function linearised in current iteration), which is usually done with singular value decomposition.

Theoretical prerequisites 10
The joint data vector is a simple concatenation of RF and SWD components, therefore joint inversion of two or more data sets doesn't differ much from inversion of a single data type.For such a data vector, one can propose a simple misfit function ( et al., 2000): where r i stands for i th of R number of RF samples, s j is a j th of S number of phase velocity values, and |∆d| 2 expresses a modified Mahalanobis distance (Mahalanobis, 1936) between random vectors in case of no correlation (we assume diagonal covariance matrix, that in general is not true, especially for RF).The parameter p is added to allow weighting the relative influence of each data set.For p = 0 we invert only RF and for p = 1 only SWD.The thing worth noticing is that the uncertainties (denoted with σ) are squared, and thus contribute to weighting of each datum more than p.
In order to stabilize the solution, it is desirable to add an extra term to the misfit function resulting in a trade-off between final data fit and smoothness of the model: where θ 2 is a positive parameter called damping, ∆m is a difference between the model in current and previous iteration, and T is a Toeplitz matrix of the following form: . (3)

Exploitation of a priori information
The "13 BB Star" experiment was conducted in the area covered by high-resolution 3D model of the crust (Grad et al., 2016).
To use this prior knowledge (Fig. 6a) in analogy to to Bayesian inference (Sambridge and Mosegaard, 2002), we include it in all starting models (Fig. 6b).We assign a lower weight to ("freeze") those layers, preventing the inversion from modifying them too much.We test two options of freezing the crust: the first (referred to as "option 1") involves the layers characterized by velocities and thicknesses exactly from 3D model (Grad et al., 2016).In the second case ("option 2") we consider thinner layers: 0.1-0.2 and 2 km for upper-middle and lower crust respectively.The velocity-depth dependence remains the same as in option 1. Densities are calculated with combined formulas from (Berteussen, 1977) and (Gardner et al., 1974).We assume 1.8-1.9km/s in sediments, 1.67, 1.73, 1.77 km/s in upper, middle, lower crust respectively, i.e. values similar to those reported by Środa et al. (1999).All mantle P-wave velocities are calculated assuming v P /v S = 1.73.The thickness of mantle layers is 5 km in every model.

Ensemble inference
Among downsides of the linearised approach, trapping by local minima is listed as one of the most profound (Bodin et al., 2012), nonetheless Monte Carlo methods are not completely free of it either (Minato et al., 2008;Wathelet, 2008).There were  several solutions to this issue proposed, e.g. in FWI a functional that pulls the model out of the local minimum (Bharadwaj et al., 2016).As the RF and SWD inversion is not as computationally expensive as FWI (Virieux and Operto, 2009;Warner et al., 2013), we take a different approach based on an ensemble of starting models, which densely cover the space of physically plausible solutions.In our case it consists of 20 homogeneous mantle structures evenly distributed between 4 and 5 km/s (Fig. 6b), within which falls the great majority of S-wave velocities reported for cratonic mantle in other studies (Fischer et al.,5 2010; Jones et al., 2010;Vinnik et al., 2015;Kind et al., 2015).Contrary to the prior distribution used in Bayesian analysis, the solution of linearised inversion can converge to the values beyond this range, which is a slightly weaker initial assumption.
Homogeneity ensures that we do not make any guess about the mantle which would drive the inversion towards the most "reasonable" structure.We rather let the data do it for us.We emphasize the fact that a high density of distribution of the starting models within given interval (here every 0.05 km/s) enables to cover the whole space of physically plausible solutions.to fit into RF tail that represents structure below the bottom of the model.For dispersion curves the reasoning behind choice of optimal maximum period T max is slightly different due to substantial widths of their sensitivity kernels for longer periods (e.g., Darbyshire et al., 2013;Romanowicz, 2002;Tsuboi and Saito, 1983).However, it leads to similar conclusions -too shallow models may be affected by sampling the structure below the model depth.

Relative weighting of data sets of different uncertainty
It may seem that the influence parameter p (Julià et al., 2000(Julià et al., , 2003) ) is solely responsible for the relative weighting of different data sets with p = 0.5 implying their equal contribution to the solution (Horspool et al., 2006).However, even more important weighting factor is represented by a squared uncertainty of a given data point (Eq.2), and the proper weighting may require value of p other than 0.5.Stringent assessment of seismic data errors in general is regarded as highly challenging (Julià et al., 2000(Julià et al., , 2003;;Bodin et al., 2012), which makes the problem of finding the optimum value of p essential.We address it by performing resolution tests for different values of p, which we describe in Sect. 4.

Tuning the initial parameters through synthetic tests
As the problem of calculating RF is highly non-linear (Bodin et al., 2012), synthetic tests of the inversion considered are not fully plausible.Even for quite similar structures there may exist a difference among the sets of inversion parameters illuminating them optimally.Here synthetic tests served to get the rough idea of the values of parameters suitable for the considered data types.In such a way we estimated the value of damping (θ 2 ≈ 1.0) and influence parameter (p ≈ 0.1), as well as a number of iterations (after 20, no significant increase in fit was observed).
4 Resolution analysis

Checkerboard test for 1D problems
We adopted a checkerboard test, (e.g., Janutyte et al., 2015) to the 1D problem.The algorithm consists of the following steps: 1. Add the pattern of alternating positive/negative anomalies of constant size (thickness) and amplitude (percentage of model velocity at given depth) to the inversion result (arithmetic mean of the ensemble of final models); Despite the simplicity of its idea, checkerboard test may be misleading as has been shown for 2D tomography problem (Leveque et al., 1993), therefore some care has to be taken also in 1D, non-tomographic case.For instance, usually only one of two possible polarities of anomaly-pattern enables its correct recovery (Fig. 7a and 7b).Moreover, for different sizes, different polarities are favourable (compare Fig. 7a and 7c).The higher amplitudes of anomalies tend to be worse recovered (Fig. 7d), which may result from excessive distortion of synthetic RF time series.

Tuning of the inversion parameters
Each parameter required in inversion was adjusted by performing a dedicated checkerboard test.The optimum set has been shown in Fig. 7. Figure 8 in turn presents the effect of choosing the wrong values, such as too short dispersion curve (Fig. 10  8a), or too low damping (Fig. 8b) and influence parameter (Fig. 8c).Crust frozen with option 2 (Fig. 8d) exhibits the presence fictitious bump and slightly worse resolving power at about 150 km depth.

Resolution matrix
The model resolution matrix R is defined by the relation: 5 and in case of damping is expressed by (Gubbins, 2004): The presence of non-zero off-diagonal elements is caused by the smoothing and the correlation of RF data.The poor crustal resolution of model for option 1 (Fig. 9a) is the effect of freezing.For option 2 there is an apparent dependence of the model  at about 140 km on the deeper structure (Fig. 9b).It seems to agree with the results of checkerboard test, which was unable to recover the anomaly at similar depth (Fig. 8d).The final loss of the resolution begins slightly lower/higher than 200 km depth for option 1 and 2 respectively.

Results
Below we present the results of the inversion (Fig. 10) including the data fit (Fig. 11 and 12) for the central station of the "13 BB Star" array (A0).The following values of inversion parameters were used: p = 0.1 (for joint inversion), θ 2 = 1.0, T max = 180 s, option 1.As shown in the previous section, the interpretation should not be made below 200 km depth.Introducing the SWD data helps to diminish the inherent non-uniqueness of the RF inversion (Fig. 10).The result for option 2 exhibits distinctive crust variations compared to option 1 (fig.13a), which, according to the checkerboard test (Fig. 8d), may not be fully plausible.
Results for all stations are presented for option 1 only (Fig. 13b).

Discussion
The common feature of mantle S-wave velocity models (see Fig. 10 and 13) is lowering velocity in the 80-120 km depth range.However, we do not interpret this zone as EEC asthenosphere, but rather, because of the neighbouring TESZ, a wedge of younger Paleozoic asthenosphere (Hoernle et al., 1995).The asthenosphere could be identified with smooth and deeper velocity lowering starting from the depth 180-200 km.Similar depth of the lithosphere-asthenosphere transition/boundary (LAB) has been obtained along profile P4 located about 200 km south-east of "13 BB Star" array.It was derived from a study of relative P-wave residuals by Świeczak (2007), using the method of Babuška and Plomerová (1992).The average depth of the LAB for the EEC was estimated at 193 km, significantly in contrast to TESZ, where it lies at 119 km depth (Wilde-Piórko et al., 2010).Numerous seismic studies have shown that the LAB beneath Precambrian cratons are not easily detected by seismic methods, e.g.receiver function, because converted phases from smooth transition zones are not well pronounced (Kind et al., 2012).Geissler et al. (2010), using S receiver functions, show that beneath the stations located in the Precambrian platform of Eastern Europe LAB deepens to approximately 200 km.This value is consistent with LAB depths of other old cratons of the Earth (Eaton et al., 2009;Jones et al., 2010).

Conclusions
Using the data from recent "13 BB Star" experiment, we studied the limitations of a linearised joint inversion of receiver function and surface-wave dispersion in terms of its ability to resolve upper mantle structure.We argue for the utility of linearised inversion whenever it is handled properly.having carefully studied the influence of inversion parameters on resolution of the final models, this type of modelling can give insights on the credibility of the solution, including approximate estimation of uncertainty.We advocate that a resolution analysis, e.g. a novel 1-D checkerboard test proposed here, ought to be an integral part of the parameter-tuning workflow, not the afterthought.This approach is versatile and can be applied to both linearised and Monte Carlo inversion scheme.Here, it allowed us to obtain reliable S-wave velocity models down to 200 km depth beneath the "13 BB Star" array.A feature that all 5 final models have in common is a presumable beginning of low-velocity zone at 180-200 km depth, which may be interpreted as an upper part of the cratonic asthenosphere.

Figure 2 .
Figure 2. Epicentral distribution of earthquakes used for calculation of receiver functions and dispersion curves of surface waves plotted relatively to a position of A0 central station of network.
files were then imported by ASWMS package and used for phase velocity calculations.The primary result of calculations are 4 Solid Earth Discuss., https://doi.org/10.5194/se-2017-58Manuscript under review for journal Solid Earth Discussion started: 14 July 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 3 .
Figure 3. Example of rotation of teleseismic P-waves and RFs for the a) C6 and b) B6 seismic stations calculated by the RF-rotation procedure (red line) and theoretical backazimuth and polarization angles (black line).Seismograms and RF are filtered with band-pass Butterworth filter of corner frequencies 0.01 and 0.8 Hz.Time zero refers to the direct P-wave.

Figure 4 .
Figure 4. Stacked RFR for all "13BB Star" broadband seismic stations.RFs were calculated based on the RF-rotation procedure.Time zero refers to the direct P-wave.

Figure 6 .
Figure 6.A priori information exploited in the inversion: (a) crustal structure beneath all stations of "13 BB Star" array computed from 3D model (Grad et al., 2016); (b) ensemble of 20 starting models with crust from central station (A0), and different homogeneous mantle velocities within physically plausible 4-5 km/s S-velocity range.
We used a 30 s wide RF window (3 s before -27 s after the direct P-wave) finding longer times not reliable enough.By analysing the synthetic RFs calculated from 2-layer models, we conclude that this would encompass converted phases down to 260-280 km depth.We made the inverted model considerably deeper (300 km) to avoid artefacts likely to appear while trying 8 Solid Earth Discuss., https://doi.org/10.5194/se-2017-58Manuscript under review for journal Solid Earth Discussion started: 14 July 2017 c Author(s) 2017.CC BY 4.0 License.

2.
Compute synthetic data from the model perturbed in above way; 3. Invert the synthetics in the exactly same way as the unperturbed model was obtained; 4. Subtract the original model from the result of previous point to recover the anomaly; 5. Compare the recovered and original anomaly.Solid Earth Discuss., https://doi.org/10.5194/se-2017-58Manuscript under review for journal Solid Earth Discussion started: 14 July 2017 c Author(s) 2017.CC BY 4.0 License.4.1.1Size, amplitude and polarity of the anomalies

Figure 10 .
Figure 10.Results of the ensemble-inversion of: (a) SWD only; (b) RF only; and (c) both data sets.

Figure 11 .
Figure 11.RF fit for all final models in the inversion of (a) both data sets; and (b) RF only.

Figure 12 .
Figure 12.SWD fit for all final models in the inversion of: (a) SWD only; and (b) both data sets.