Anisotropic tomography of the European lithospheric structure from surface wave studies

We present continental‐scale seismic isotropic and anisotropic imaging of shear wave upper‐mantle structure of tectonically diversified terranes creating the European continent. Taking into account the 36–200 s period range of surface waves enables us to model the deep subcontinental structure at different vertical scale‐lengths down to 300 km. After very strict quality selection criteria, we have obtained phase wave speeds at different periods for fundamental Rayleigh and Love modes from about 9000 three‐component seismograms. Dispersion measurements are performed by using Fourier‐domain waveform inversion technique named “roller‐coaster‐type” algorithm. We used the reference model with a varying average crustal structure for each source‐station path. That procedure led to significant improvement of the quality and number of phase wave speed dispersion measurements compared to the common approach of using a reference model with one average crustal structure. Surface wave dispersion data are inverted at depth for retrieving isotropy and anisotropy parameters. The fast axis directions related to azimuthal anisotropy at different depths constitute a rich database for geodynamical interpretations. Shear wave anomalies of the horizontal dimension larger than 200 km are imaged in our models. They correlate with tectonic provinces of varying age‐provenance. Different anisotropy patterns are observed along the most distinctive feature on our maps–the bordering zone between the Palaeozoic and Precambrian Europe. We discuss the depth changes of the lithosphere‐asthenosphere boundary along the profiles crossing the chosen tectonic units of different origin and age: Fennoscandia, East European Craton, Anatolia, Mediterranean subduction zones. Within the flat and stable cratonic lithosphere, we find traces of the midlithospheric discontinuity.


Introduction
Surface-wave tomography is based on the frequency-dependent wave speeds in a dispersive medium. A layered structure and the presence of the free surface perfectly fulfill the conditions for the excitation of Rayleigh and Love waves near the surface of the earth. High quality waveforms give us the possibility to retrieve dispersion properties which can be translated into the shear-wave seismic structure of the region under investigation.
A massive collection of seismic waveforms was recorded at a growing number of broadband stations in Europe over the last 15 years. Dense and more complete path coverage of the study volume compared to previous studies motivated us to extract the information contained in radial and transverse components of surface waves in order to revise the sublithospheric structure from the point of view of their anisotropic properties. Highly dispersive characteristics of the medium ensure higher vertical resolution of the mantle models than the ones obtained using body waves. A lot of debate concerns separating the anisotropy effect from the heterogeneity since both of them can have the same order of magnitude [Montagner and Guillot, 2002]. Till the end of nineties, the tomographic studies on seismic anisotropy in Europe were mostly based on body wave P-residuals. Including anisotropy, improved the correlation of isotropic wave speed perturbation models with the thickness of the lithosphere [Plomerov a, 1997].
A wide variety of approaches aimed at exploring the information from surface waves to study European mantle structure have been reported in literature.
In the middle of the last decade, the Automated Multimode Inversion (AMI) method has been elaborated by Lebedev et al. [2005] in order to treat large data sets of surface and S waves. Special attention was focused on removing the scattered waves from the signal. Thanks to filtering in a series of time-frequency Some authors focused their mantle tomographic studies on extracting the information from group wave speed databases. A first example of that case provides transversely isotropic models of the European mantle which evolved toward enhanced resolution models Morelli, 2009, 2011] by incorporating more and more broadband array data.
Other seismologists [Boschi et al., 2009;Schaefer et al., 2011] have undertaken to parameterize the model grid in an adaptive way which takes into account the changing lateral resolution of the data. These changes are mainly dictated by unevenly distributed data coverage.
The dramatic increase of processing capacities in the last 10 years led to modeling of the waveforms in a broad period range going down to 10 s. The results of filling the gap between the travel-time tomography and the inversion of surface wave dispersion properties are presented in the work of Zhu et al. [2012]. By using a spectral-element method (SEM), the authors model three-component seismograms to fit phase anomalies for different frequencies. They obtained P-wave isotropic and S-wave anisotropic model which map tectonic deformations [Zhu and Tromp, 2013]. The last results of the application of adjoint tomography technique to image European upper mantle are presented in the paper by Zhu et al. [2015]. Computational costs of the authors' approach are very demanding since a big parameter space must be inverted for retrieving the elastic and anelastic part of the model.
Interesting remedy has been proposed by [Lekić and Romanowicz, 2011] to reduce high computational costs of the inversion procedure. The tomographers proposed a hybrid approach in which sensitivity kernels are calculated using normal mode summation theory and then spectral-element method is applied to calculate the synthetics. Different approaches have been proposed to retrieve phase wave speeds from the dispersive surface wave trains. In global studies, the seismologists explore the information contained in fundamental modes to image lithospheric structures, as presented in the recent work by Burgos et al. [2014], or they include higher modes to reach into the mantle transition zone. The second case example is a work by Debayle and Ricard [2012] in which automated waveform inversion, reinforced with cross-correlation techniques, has been applied to Rayleigh fundamental and five higher order harmonics to image the mantle SV structure.
In this study, we use the information extracted from dispersive nature of Rayleigh and Love waveforms. Differently polarized types of surface waves give complementary information about model parameters. The starting database of around 16,000 three-component seismograms was thoroughly quality checked to ensure the time separation of the events being examined. Strict quality checking criteria applied after the calculation of phase velocity dispersion curves constrained our database to about 9000 seismograms. To measure surface-wave fundamental mode phase wave speeds, we apply a ''roller-coaster'' technique [Beucler et al., 2003] which ensures that fundamental mode dispersion calculations are not biased by higher modes.
The reference model for calculating the synthetics takes into account the average crustal thickness for each source-receiver path. We use the CRUST1.0 [Laske et al., 2013] model with a data resolution of 18 by 18. After applying the regionalization method by Montagner [1986] and Debayle and Sambridge [2004], we obtain local phase wave speeds at the grid points in the study region. Then we follow the isotropic depth inversion procedure in which we take into account the local crustal structure from CRUST1.0. Finally, the first order perturbation theory is applied to retrieve azimuthal anisotropy parameters from Rayleigh wave database. We derive a 3D model of the European mantle down to a depth of 300 km.

Data
To create our surface waveform database for European tomographic studies, we have selected and checked the signal quality of the records of events with body wave magnitude higher than 5.5 with the epicenters Geochemistry, Geophysics, Geosystems The instrumental response was removed from the raw data to give bandpassed (0.003-0.1 Hz) ground velocity seismograms rotated into Radial/Transverse/Vertical components.
We have obtained about 16,000 seismograms for which signal-to-noise ratio was high enough to retrieve fundamental mode amplitude from the background noise. This database has been subsequently constrained to around 9000 after applying additional quality criteria while retrieving phase wave speeds as presented in section 2.2.1. Table 1 presents the numbers of the good quality selected Rayleigh and Love wave waveforms recorded at the networks. The ray coverage of the region by Rayleigh and Love waves is presented in Figure 1.

Modeling
The workflow is divided into the following steps: 1. Calculation of synthetic reference seismograms using normal mode summation (nms) method and taking into account the variations in crustal thickness for each event-station path (period range 36-500 s); 2. Inversion of real data to obtain phase wave speed dispersion curves for Rayleigh and Love waves; 3. Regionalization of the geographically unevenly distributed phase wave speed measurements for the chosen periods; 4. Isotropic depth inversion of Rayleigh wave speeds to obtain a starting model for anisotropic inversion; 5. Anisotropic depth inversion of Rayleigh wave data set to retrieve azimuthal anisotropy parameters.
Following the arguments presented by Ritzwoller and Levshin [1998], we assume the propagation of the rays along the great circles as fully justified by the range of periods being treated (higher than 36 s) and the length-scale of the paths.

Phase Wave Speed Measurements
Prior to proceeding to calculation of the fundamental mode dispersion curves, all the seismograms were windowed in order to separate the X-phases (including the higher modes) and the fundamental mode.
The synthetic seismograms were calculated in the PREM model [Dziewonski and Anderson, 1981] for each event-station path using the normal mode summation method [Gilbert, 1971]. As surface waves are sensitive to an average structure along the path, we found justified to improve our starting model by assigning average crustal thickness, rather than composition, between each source-receiver path. The average thickness of the crust was calculated from the CRUST1.0 model [Laske et al., 2013]   Geochemistry, Geophysics, Geosystems

10.1002/2015GC006243
The eigenfunctions calculated separately for spheroidal and toroidal modes include the higher modes with radial order up to 10 and angular order up to 1500. The period range of the eigenfunctions (10-1000 s) is wider than the frequency range of sources used in our modeling (36-500 s). While calculating the reference dispersion curves based on our eigenfunctions the Stoneley modes are excluded. In our long-period waveform modeling, the source parameters are provided by the moment tensor solution from the CMT Harvard database [Ekstr€ om et al., 2012]. In order to eliminate the influence of the recorder characteristic from the seismogram, we deconvolve the real data with the instrumental response function to get ground velocity.
The so-called roller-coaster-type method [Beucler et al., 2003] is applied to calculate phase wave speed dispersion curves. The spectra of the real data and of the synthetics are calculated for the time window of the seismograms corresponding to fundamental and higher modes, separately. For the aim of this study, we focus on our fundmental mode database. In the inversion procedure the misfit function between the real and synthetic spectra is minimized after applying a gradient least-squares algorithm to compare the spectra. Using ''roller-coaster'' method, we start exploring the space of model parameters looking for large-scale variations and we finish with retrieving short-wavelength solution in local variations of the misfit function. The effect of phase indeterminacy of 2p is minimized in the approach we use. The reasons for applying the ''rollercoaster'' method were the following: 1. It assures efficient separation of fundamental modes from the higher ones on the seismograms; 2. It provided us the disperion curves for single stations which means that our choice of source distributions was larger than in a case of interstation measurements; 3. We modified the method by applying a starting model that takes into account the average crustal thickness between the source and the receiver in order to improve the quality of dispersion data measurements.  show that in the case of Love waves the discussed limits are lower. The different seismograms presented in Figure 2 (black-real data, green-synthetic calculated by nms method, red-synthetic after applying a phase correction calculated from the retrieved dispersion curves) and the phase matching of real and synthetic data are a key argument of the accuracy of the chosen approach. That fitting was significantly improved by considering the varying crustal thickness in the chosen reference model.
In Figure 3, we present an example of a synthetic test made in order to evaluate the influence of crustal thickness on phase wave speed calculations. The left panel of the figure refers to modeling with the starting reference PREM model with average 27 km-thick crust and the right one to the starting model with the crust thickness of 35 km used to obtain synthetic (real) data. The closest the initial crustal thickness to the real model, the best resolution in recovering of a dispersion curve is obtained.
After calculation of phase wave speeds, we have correlated phase-corrected synthetics with real data in a time window corresponding to fundamental mode (marked as traingles in Figure 2). Our criteria for accepting the seismogram was the value of correlation coefficient higher than 90% and the dispersion measurements within a range ðl12rÞ and ðl22rÞ, where l denotes data mean and r standard deviation of the wave speed data. As a result of the quality checking procedure, we have constrained our final database to approximately 9000 seismograms.

Regionalization and Depth Inversion
The procedure of regionalization resulted in evenly spaced phase wave speed mesurements presented as phase wave speed maps. We use the technique of regionalization developed and presented in detail in Montagner [1986] and improved by Debayle and Sambridge [2004]. The least square inversion formalism provides local phase wave speed at each point of the chosen grid. The local phase velocity Cðx; H; U; WÞ can be expanded as a Fourier series of azimuth W: Cðx; WÞ5A 0 ðxÞ1A 1 ðxÞcosð2WÞ1A 2 ðxÞsinð2WÞ1A 3 ðxÞcosð4WÞ1A 4 ðxÞsinð4WÞ, where A i ðxÞ are depth integral functions related by Montagner and Nataf [1986] to Cartesian elastic coefficients c ij . A 0 term refers to the equivalent tranversely isotropic medium with vertical symmetry axis (no azimuthal anisotropy). For Rayleigh waves 2W terms are dominant whereas for Love waves 4W terms are more significant.
We applied spatial correlation functions to our data set which gave us the model parameters and their errors defined on each point of the grid (in our case 18 by 18). At each point of the grid we obtained the value of the amplitude as well as the direction of the azimuthal anisotropy. The chosen spatial correlation length of 200 km, defining the smoothness of the final model, was a compromise between: 1. The extension of the region being modeled; 2. Geographical distribution of ray paths; 3. Tradeoff between the wavelengths and the scale of heterogeneity; 4. The size of a grid.
The value of a correlation length is below the seismic wavelength at long period but as far as the ballistic wave is dominant, the geometrical optics is still a good approximation [Montagner, 1986].
We have made tests to assesss how variance reduction changes with the horizontal correlation length. Following the arguments for estimating the optimal correlation length in Montagner and Jobert [1988] and Griot et al. [1998] we calculated the changes in variance reduction for a set of correlation length (100 km, 150 km, 200 km, 250 km, 300 km, 350 km, 400 km) referring to the different periods of Rayleigh and Love waves being regionalized. As shown in Figure 4, the minimum acceptable correlation length for our models is 200 km. The lower values do not improve the variance reduction. The horizontal axis unit is a rescaled value of L due to the fact that the number of parameters in modeling is inversely proportional to the square of the correlation length.
As a result of regionalization, we have obtained phase wave speed maps for Rayleigh waves for the period range 36-200 s as well as for Love waves for the periods 36-150 s with 5 s step interval between the subsequent periods. The choice of this period range was justified by the sufficient path density distribution in the study region. The examples of path density maps for the chosen periods of Rayleigh and Love waves are shown in Figure 5. The blue scale intensity corresponds to the number of ray segments crossing 1 o by 1 o cells in the presented geographical projection. Central and Western Europe is densely Geochemistry, Geophysics, Geosystems 10.1002/2015GC006243 covered with rays coming from a wide range of azimuths (more than 200 counts in a cell) which justifies retrieving azimuthal anisotropy parameters during modeling. Figure 6 presents the azimuthal distribution of Rayleigh and Love data which is crucial in retrieving the azimuthal anisotropy parameters. The azimuthal distribution of the events is almost uniform in Central Europe and satisfying for the adjacent parts in the case of Rayleigh wave database. For Love wave database, the records from northern and southern directions are scarce. The sensitivity kernels presented in Figure 7 show that the chosen period range of Figure 3. Synthetic test to show differences in phase wave speed dispersion fitting to synthetic (real) data for a starting model with (left) average PREM crust thickness and (right) the crust 8 km thicker than global average. Black dashed line seismograms-synthetic data to be recovered in modeling, green ones-seismograms for initial model, red ones-final seismograms. The curves below the seismograms show the improvement in the accuracy of phase wave speeds recovery for longer periods in a case the initial crustal thickness is closer to the real one.

10.1002/2015GC006243
36-150 s enables to retrieve the depth structure for Rayleigh waves (V SV -phase wave speed of vertically polarized shear waves) and for Love waves (V SH -phase wave speed of horizontally polarized S waves) down to 250 km. The sensitivity kernels are calculated with respect to two of the so-called Love parameters [Takeuchi and Saito, 1972] defined as L5qV 2 SV and N5qV 2 SH . L determines the wave speed of SV-wave whereas N the wave speed of the SH-wave. The radial anisotropy can be retrieved form the formula w5N=L. The high computational cost of calculating the a posteriori covariance matrix does not enable us to calculate errors assigned to phase wave speed measurements. Our remedy was to generate statistical errors using the bootstrapping method. In case of Rayleigh wave database we randomly sample 80% of the data 30 times whereas 90% of Love wave data is sampled 40 times. The standard deviations of the stacked samples are treated as statistical errors. Figure 10 presents the estimated errors for the chosen periods of Rayleigh waves.
For the inversion at depth, we start with an isotropic inversion and then perform the anisotropic one. The isotropic inversion was based on the so-called ''trans-dimensional inversion'' approach described in Haned et al. [2016] and implemented in the SWAVETRIN program [Haned et al., 2016]. This approach consists in automatic adapting of model parameterization to phase wave speed uncertainty. At each grid point, an a priori 1-D model is defined. The crustal part of that model is a local CRUST1.0 model and the deeper part is PREM reference model in which 220 km discontinuity has been replaced by a smooth gradual change in wave speed. There are two loops in the procedure: the inner one (computing the optimum model weighting coefficients for a given spline basis by minimizing the misfit between between real and synthetic velocity measurements) and the outer one (determining the optimum spline basis to minimize the data misfit and the S-wave speed model uncertainty).
The anisotropic depth inversion was based on the first order perturbation theory as presented by Montagner and Nataf [1986] and on the least-squares inverse algorithm. The output of the inversion is an a The results of inversion procedure are shown in Figure 11. The depth maps give also the information about the direction of the fast axes of azimuthal anisotropy direction (W parameter) in the study region.

Modeling Results and Interpretation
In this section, we focus on distinctive tectonic features imaged in tomographic sections (see Figure 11). Shear-wave wave speed cross sections along different profiles cutting age-diversified tectonic units within the European plate are shown in Figure 12. The color palette facilitates quick distinction between the perturbations ascribed to the lithosphere positive (blue) and the asthenospheric negative (red) anomalies. The lower lid of blue anomalies corresponds to the lithosphere-asthenosphere boundary (LAB) modeled using our surface wave database. It corresponds to a sudden drop in S-wave speed of the order of a few percent.
The prominent feature seen in phase wave speed maps as well as on depth inverted S-wave velocity maps is a strong contrast between Precambrian and Palaeozoic Europe. The huge suture tectonic unit referred to as TESZ (Trans European Suture Zone) is clearly imaged for Rayleigh periods 50-150 s which correspond to a depth range 50-200 km. Such deep penetration of the feature is a strong argument that TESZ has subcrustal prolongation. Deep seismic profiling in the region reaches a depth of 80 km but recent teleseismic tomography studies of Janutyte et al. [2014], based on passive seismic experiment data, resulted in P-wave perturbation distributions down to 300 km. Amplitude and depth range of positive and negative S-wave anomalies observed in the depth section crossing the TESZ are in good agreement with the presented study. The horizontal width of the suture area in our modeling is of the order of 150-200 km. It is clearly imaged in the profiles CD, EF and GH (Figure 12) at the corresponding distances of 2500 km, 1800 km and 3000 km. Their directions follow the SW direction of docking of the suspect terranes to the East European Craton. The recent studies of the anisotropy in the TESZ based on SKS measurements [Vecsey et al., 2014] show that across the suture zone, the changes in fast axes directions of anisotropy are smooth whereas different anisotropy patterns are observed along the feature. The authors postulate the possibility of a penetration of the lithospheric cratonic mantle beneath the Palaeozoic platform units. In our depth sections, two different anisotropy patterns (W-E in the northern part and SW-NW in the southern one) can to be easily identified for a depth range of 70-120 km. At a depth of 150 km that pattern vanishes and the fast axis direction of azimuthal anisotropy is parallel to the suture. The changes across the zone cannot be assessed due to the limits of horizontal data resolution comparable with the width of the suture. The lithosphere depth of the two sides of the zone changes dramatically-from 80 km below the tectonically younger units to 200-220 km beneath the cratonic part. The Central European region around the TTZ has been imaged with the resolution of 150 km by Soomro et al. [2016] with the use of broad-band interstation measurements of Rayleigh and Love waves. The phase velocity maps are presented by authors up to a period of 60 s. A prominent negative shear wave anomaly of the oder of 3% in the lithospheric mantle of the Bohemian Massif is present in that study as well as in our model. In both models, the velocity contrast on the two sides of the TTZ zone is the most sharp for the periods higher that 40 s which translates to subcrustal structure.
The architecture of the lithosphere in Eastern Europe is dominated by the old cold Precambian craton referred to as East European Craton (EEC) which formed a core of the Baltica continent. The oldest Precambrian rocks of c. 3700 Ma occur in Ukraine and the Baltic Shield. The East European Platform is far from being a homogeneous structure. At a depth of around 80-100 km (see profiles: AB, CD, EF, GH) we observe a zone of decreased wave speed which the most probably traces the midlithospheric discontinuity (MLD). According to the results published by Selway et al. [2015], long period data, which is our case, are more suitable for tracing in a coherent way the wave speed contrast at the MLD than the body wave data. The feature being discussed refers to the thick cratonic lithosphere (200-220 km). The studies done by Yuan et al. [2011] in a thick North American craton give arguments for the MLD having been resolved at a depth of about 100 km. For the time being, one of the most plausible explanation of the origin of that discontinuity Geochemistry, Geophysics, Geosystems 10.1002/2015GC006243 seems to the presence of amphibole-rich layers in metasomatized rocks. Nor the thermal effect neither partial melting itself can explain the observed amplitude of seismic wave speed change at the MLD depth [Selway et al., 2015]. Gradual changes of temperature corresponding to cratonic depths could not result is observed discontinuity in wave speed distribution. As for a partial melt, the only solidus which crosses cratonic geothermal at the MLD depth is not representative of volatile-poor cratons [Griffin et al., 2009].
Our model supports the interpretation by Legendre et al. [2012] that there is no indication for cratonic roots below 280 km.
The results of our anisotropic modeling in the Baltic Shield seem to give arguments in favor of the hypothesis concerning the creation of the cratonic lithospheric mantle in the Fennoscandia presented by Plomerov a  shows deep lithospheric cratonic roots beneath the middle part of the Baltic Shield. The part of the profile crossing the southern Caledonides shows no considerable thickening of the lithosphere beneath the orogen.
The complexities in the tectonics of the Mediterranean region result from the still ongoing collision of Africa with Euroasia and the closure of the Tethyan ocean. Complicated subduction geometry pattern typical of the Mediterranean Basin [Faccenna and Becker, 2010] is represented in depth section as adjacent areas of negative and positive anomalies. Subducted bodies are clearly imaged on profile IJ and KL.
Interesting studies have been dedicated to describe the anisotropic seismic structure of the Aegean crust and mantle. Endrun et al. [2011] propose three-layer deformation pattern in the crustal and mantle area bordering the Hellenic slab from the northern direction. Their Rayleigh wave data span the periods significantly lower than ours but the lowest part of their model, the asthenospheric mantle, can be directly compared to the top of our model. The azimuthal anisotropy direction between 70 and 100 km shows SW-NE direction which is consisted not only with the results of Endrun et al. [2011] but also of Zhu and Tromp [2013]. Below the asthenosphere the pattern being considered changes to NW-SE. The layered pattern of azimuthal anisotropy would not have been possible to be obtained using SKS splitting measurements which are less sensitive to depth than the studies using dispersive feature of surface waves.
Three negative anomalies at a depth of 120 km correspond to Neogene extensional basin. A representative example of mantle upwelling in the Mediterranean is presented in profile OP starting in the Calabrian Arc subduction zone (positive wave speed anomalies) and then going along the east side of the Appenines arc to the Massif Central. We observe a lithosphere thinning along the arc structure referring to the upwelling of the mantle as discussed in Carminati and Doglioni [2012].
Profile MN transects the border between the Anatolian and Aegean Plate. In the left corner of the discussed profile, a fragment of the Hellenic slab is clearly imaged. The modeling results beneath the area of Turkey (profile MN) show no indication of the lithosphere thickening beneath the Anatolian Plate. Relatively shallow lithosphere-asthenosphere boundary, seen in our model at a depth of 80-100 km, corresponds to the recently published results of S-receiver functions studies beneath Turkey [Kind et al., 2015]. Seismic evidence of mantle lithosphere absence in that region is in agreement with geodynamical modeling by G€ o g€ uş and Pysklywec [2008]. The authors support a view that lithospheric thinning, plateau uplift and synconvergent extension are caused by the delamination of the mantle lithosphere. Recent results of multiscale full seismic waveform inversion in the region, published by Fichtner et al. [2013], show shallow Anatolian asthenosphere connected to surface instances of volcanic provinces.
Our perspective is to enrich our fundamental mode database by higher mode dispersion measurements in order to reach in the joint-inversion modeling the mantle transition zone depths.

Concluding Remarks
The output of our modeling presents a set of seismic parameters related to the European lithosphere which can be directly compared with the results of regional-scale modeling done using different methods of surface as well as body wave treatment. We have obtained high resolution image, with lateral resolution of about 200 km and vertical one of the order of 50 km, which is not typical of other continental-scale models. We present azimuthal anisotropy direction changes with depth which can enrich the discussions concerning last and present deformation pattern of the European lithosphere. The obtained resolution led to identify regional-scale features reflecting complex geodynamical pattern of Europe, especially in the Mediterranean region. By presenting the seismic evidence of midlithospheric discontinuity within the stable East-European craton, we join a discussion concerning the processes associated with continent formation.