# Anisotropic tomography of the European lithospheric structure from surface wave studies

Blanka Nita, Satish Maurya, Jean-Paul Montagner

► To cite this version:

Blanka Nita, Satish Maurya, Jean-Paul Montagner. Anisotropic tomography of the European lithospheric structure from surface wave studies. Geochemistry, Geophysics, Geosystems, 2016, 17, pp.2015-2033. 10.1002/2015GC006243. insu-03581659

HAL Id: insu-03581659

https://insu.hal.science/insu-03581659

Submitted on 20 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

10.1002/2015GC006243

**Key Points:**

* We present azimuthally anisotropic model of the European upper mantle down to a depth of 300 km with a high lateral resolution of 200 km

* We find traces of the midlithospheric discontinuity within the cratonic lithosphere in Eastern Europe

* Anisotropic pattern suggests block structure of the Fennoscandia

**Correspondence to:**

B. Nita,
nita@ipgp.fr

**Citation:**

Nita, B., S. Maurya, and J.-P. Montagner (2016), Anisotropic tomography of the European lithospheric structure from surface wave studies, *Geochem. Geophys. Geosyst.*, 17, 2015–2033, doi:10.1002/2015GC006243.

Received 30 DEC 2015
Accepted 9 MAY 2016
Accepted article online
Published online 10 JUN 2016

# Anisotropic tomography of the European lithospheric structure from surface wave studies

Blanka Nita¹, Satish Maurya¹, and Jean-Paul Montagner¹

¹Institut de Physique du Globe de Paris, 1 rue Jussieu, Paris, France

**Abstract** 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.

## 1. 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á, 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
windows, the authors proposed efficient fitting quality criteria of simplified waveforms and the synthetics. The AMI has been successfully applied to image mantle structure beneath Europe by Legendre et al. [2012].

Ambient noise Rayleigh wave tomography for the European continent [Yang et al., 2007] shows correlation with tectonic units but the depth of analyses is limited by the noise range periods (10–50 s) to sedimentary and crustal depths.

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 [Schivardi and 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 1° by 1°. 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.

## 2. Data and Modeling Technique

### 2.1. 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
<table><caption>**Table 1.** Rayleigh and Love Wave Statistics for the Networks Considered in the Databaseª</caption><thead><tr><th>NETWORK</th><th>NS</th><th>NR</th><th>NL</th></tr></thead><tbody><tr><td>BE</td><td>1</td><td>108</td><td>52</td></tr><tr><td>BN</td><td>4</td><td>30</td><td>2</td></tr><tr><td>CH</td><td>9</td><td>778</td><td>413</td></tr><tr><td>CZ</td><td>2</td><td>438</td><td>194</td></tr><tr><td>DK</td><td>1</td><td>125</td><td>19</td></tr><tr><td>G</td><td>4</td><td>1192</td><td>452</td></tr><tr><td>GB</td><td>7</td><td>81</td><td>31</td></tr><tr><td>GE</td><td>50</td><td>2469</td><td>486</td></tr><tr><td>GO</td><td>5</td><td>35</td><td>9</td></tr><tr><td>GR</td><td>2</td><td>288</td><td>32</td></tr><tr><td>HL</td><td>2</td><td>288</td><td>32</td></tr><tr><td>HT</td><td>4</td><td>52</td><td>1</td></tr><tr><td>II</td><td>5</td><td>2454</td><td>906</td></tr><tr><td>IU</td><td>7</td><td>1467</td><td>1372</td></tr><tr><td>KO</td><td>10</td><td>294</td><td>18</td></tr><tr><td>LX</td><td>4</td><td>171</td><td>47</td></tr><tr><td>MN</td><td>17</td><td>1407</td><td>321</td></tr><tr><td>NI</td><td>4</td><td>253</td><td>55</td></tr><tr><td>NO</td><td>1</td><td>53</td><td>0</td></tr><tr><td>OE</td><td>3</td><td>259</td><td>176</td></tr><tr><td>PL</td><td>4</td><td>30</td><td>8</td></tr><tr><td>PM</td><td>3</td><td>209</td><td>39</td></tr><tr><td>SS</td><td>1</td><td>137</td><td>26</td></tr><tr><td>TT</td><td>2</td><td>83</td><td>41</td></tr><tr><td>XG</td><td>28</td><td>254</td><td>90</td></tr><tr><td>XK</td><td>23</td><td>372</td><td>99</td></tr><tr><td>XS</td><td>37</td><td>70</td><td>44</td></tr><tr><td>XY</td><td>1</td><td>18</td><td>5</td></tr><tr><td>Y1</td><td>4</td><td>70</td><td>5</td></tr><tr><td>Y4</td><td>1</td><td>7</td><td>0</td></tr><tr><td>YL</td><td>7</td><td>33</td><td>27</td></tr><tr><td>YY</td><td>6</td><td>46</td><td>1</td></tr><tr><td>ZU</td><td>10</td><td>104</td><td>1</td></tr><tr><td>3C</td><td>5</td><td>17</td><td>9</td></tr><tr><td>AC</td><td>1</td><td>3</td><td>1</td></tr></tbody></table>

within the rectangle of the longitudinal extension –80W to 100E and latitudinal –25S to 75N. The number of earthquakes fulfilling those criteria is 1292. All seismic stations (275 in total) are located on the European plate.

We have retrieved the surface wave parameters for the waveforms recorded at the following 35 permanent and temporary networks in Europe: BE, BN, CH, CZ, DK, G, GB, GE, GO, GR, HL, HT, II, IU, KO, LX, MN, NI, NO, OE, PL, PM, SS, TT, XG, XK, XS, XY, Y1, Y4, YL, YY, ZU, 3C, AC.

The period of data request span a time period of 13 years (1 January 2000 to 30 June 2014).

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.

## 2.2. 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.

### 2.2.1. 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] which is the updated version of CRUST2.0 [Laske et al., 2001] with a data resolution of 1° by 1°.

<footnote>ªNS—number of stations; NR—number of seismograms with good quality Rayleigh waves; NL—number of seismograms with good quality Love waves.</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Ray coverage of the study region by Rayleigh (upper part) and Love (lower part) waves. Yellow stars-epicenters; red circles-stations.</figcaption></figure>

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öm 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 fundamental 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 $2\pi$ is minimized in the approach we use. The reasons for applying the “roller-coaster” method were the following:

1. It assures efficient separation of fundamental modes from the higher ones on the seismograms;

2. It provided us the dispersion 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.

Figure 2 displays examples to show how the reference model synthetics fit the real data. The longer the epicentral distance the longer the maximum period retrieved by the applied technique is. The lower limit of distance, of the order of 2000 km, enables us to retrieve the phase dispersion curves up to about 150 s (wavelength of the order of 600 km), whereas the seismograms recorded at epicentral distances of 6000 km enable us to retrieve the phase velocity up to the period of 250 s (wavelength of the order of 1000 km). The crucial parameter determining the retrieved period range is the event magnitude. The presented examples
<figure><img src="image_5.png" /></figure>

Mb=5.8

<figure><img src="image_6.png" /></figure>

delta=6159 km

Mb=6.8

**Figure 2.** Examples of phase wave speed dispersion curves (red curves) calculated using the roller-coaster method for two different epicentral distances. Black dashed line seismograms—real data, green ones—synthetics, red ones—phase-corrected synthetics. The curves below the dispersion ones show the phase wave speed residua and the modeling errors. Dashed line dispersion curve refers to the PREM model with average crustal structure calculated for each raypath separately. Left column—Rayleigh waves, right column—Love waves.
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 (μ+2σ) and (μ-2σ), where μ denotes data mean and σ 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.

### 2.2.2. Regionalization and Depth Inversion

The procedure of regionalization resulted in evenly spaced phase wave speed measurements 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(\omega, \Theta, \Phi, \Psi)$ can be expanded as a Fourier series of azimuth $\Psi$:

$$C(\omega, \Psi) = A_0(\omega) + A_1(\omega)\cos(2\Psi) + A_2(\omega)\sin(2\Psi) + A_3(\omega)\cos(4\Psi) + A_4(\omega)\sin(4\Psi), \text{ where } A_i(\omega) \text{ are depth integral functions related by Montagner and Nataf [1986] to Cartesian elastic coefficients } c_{ij}. A_0 \text{ term refers to the equivalent tranversely isotropic medium with vertical symmetry axis (no azimuthal anisotropy). For Rayleigh waves } 2\Psi \text{ terms are dominant whereas for Love waves } 4\Psi \text{ 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 1° by 1°). 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 assess 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° by 1° cells in the presented geographical projection. Central and Western Europe is densely
<figure><img src="image_7.png" /></figure>

synthetic test, event depth=15 km, Mb=5.5, delta=6250 km

<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

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 3 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><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Variance reduction for different periods and correlation lengths. Top: Rayleigh waves. Bottom: Love waves.</figcaption></figure>

36–150 s enables to retrieve the depth structure for Rayleigh waves (V<sub>SV</sub>–phase wave speed of vertically polarized shear waves) and for Love waves (V<sub>SH</sub>–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 parame- ters [Takeuchi and Saito, 1972] defined as L = ρV<sup>2</sup><sub>SV</sub> and N = ρV<sup>2</sup><sub>SH</sub>. L determines the wave speed of SV-wave whereas N the wave speed of the SH-wave. The radial anisotropy can be retrieved from the formula ψ=N/L.

Figures 8 and 9 show the chosen period sections of phase wave speed perturbations of Rayleigh and Love waves, respectively. The perturbations are calculated with respect to the median values of the data sets. The azimuthal anisotropy values and their directions are superimposed on the S-wave speed maps.

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
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 5.** Path density maps for the selected periods of Rayleigh and Love waves presenting number of rays crossing 1° by 1° cells. The color scale is saturated at a value of 100 rays.</figcaption></figure>

posteriori model with the error values. Diagonal terms of the a posteriori covariance matrix are used to estimate a final error on the model parameters. The inversion procedure enables to give Gaussian vertical correlation lengths changing with depth.

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 (Ψ parameter) in the study region.
<figure><img src="image_13.png" /><figcaption>**Figure 6.** Rose diagrams of azimuthal coverage of the regionalized data for Rayleigh and Love waves.</figcaption></figure>

### 3. 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.
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 7.** Kernels presenting the depth sensitivity of (top) Rayleigh and (bottom) Love wave phase wave speed to L and N parameters. Different colors correspond to the following periods: red=36 s, blue=50 s, black=100 s, green=150 s.</figcaption></figure>

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 order 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 than 40 s which translates to subcrustal structure.

The architecture of the lithosphere in Eastern Europe is dominated by the old cold Precambrian 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 explanations of the origin of that discontinuity
<figure><img src="image_16.png" /><figcaption>**Figure 8.** Phase wave speed maps for the selected periods of Rayleigh waves. Plate boundaries are marked as brown curves.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 9.** Phase wave speed maps for the selected periods of Love waves. Plate boundaries are marked as brown curves.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 10.** Errors calculated using the bootstrapping method for the selected periods of Rayleigh waves.</figcaption></figure>

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 [Sel-way 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á
<figure><img src="image_19.png" /><figcaption>**Figure 11.** Depth maps of Rayleigh wave isotropic Vsv wave speeds. The bars indicate the fast axes of azimuthal anisotropy.</figcaption></figure>
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 12.** Cross sections showing S-wave wave speed field along the chosen profiles. The color palette was adjusted to show a contrast between the lithosphere (blue) and asthenosphere (red). The depth axes start at 40 km. The tectonic sketch map was modified after Encyclopedia Britannica.</figcaption></figure>

et al. [2011]. The authors model the region as the assemblage of blocks with different anisotropy pattern changes. Their model rather excludes the formation of the cratonic lithospheric mantle as progressive thickening of the lithospheric roots. In our depth sections, the Fennoscandian anisotropy pattern suggests block structure imaged down to deep lithospheric depths. The region of the extreme positive wave speed anomaly in Fennoscandia (see depth sections in Figure 11) down to 250 km correlates with the oldest Archean part of the Baltic Shield with the thickened lithospheric roots. Profile AB in Figure 12
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öğüş 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.

## 4. 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.

## References

* Beucler, E., É. Stutzmann, and J.-P. Montagner (2003), Surface wave higher-mode phase velocity measurements using a roller-coaster-type algorithm, *Geophys. J. Int.*, **155**(1), 289–307.

* Beyreuther, M., R. Barsch, L. Krischer, T. Megies, Y. Behr, and J. Wassermann (2010), Obspy: A Python toolbox for seismology, *Seismol. Res. Lett.*, **81**(3), 530–533.

* Boschi, L., B. Fry, G. Ekström, and D. Giardini (2009), The European upper mantle as seen by surface waves, *Surv. Geophys.*, **30**(4-5), 463–501.

## Acknowledgments

The calculation of the synthetic data and phase wave speed dispersion curves was done with the use of the modeslib0.2 package modified by Yann Capdeville after *Woodhouse and Dziewonski* [1984]. The massive collection of seismograms was downloaded using a python toolbox: obspyDMT [Scheingraber et al., 2013]. The data processing was made using the SAC analysis tools [Goldstein et al., 2003] and the ObsPy toolbox [Beyreuther et al., 2010]. The maps and waveforms were visualized with the use GNU licence GMT software [Wessel and Smith, 1998] and the gsac application [Herrmann et al., 2004]. We highly appreciate the discussion with Eleonore Stutzmann and Gaël Burgos concerning the methodological aspect of our work. The work was supported by the Marie Curie grant F108 P7-PEOPLE-2011-IEF-JIBSA (grant agreement 300034).
Burgos, G., J.-P. Montagner, E. Beucler, Y. Capdeville, A. Mocquet, and M. Drilleau (2014), Oceanic lithosphere-asthenosphere boundary from surface wave dispersion data, *J. Geophys. Res. Solid Earth*, 119, 1079–1093, doi:10.1002/2013JB010528.

Carminati, E., and C. Doglioni (2012), Alps vs. Apennines: The paradigm of a tectonically asymmetric Earth, *Earth Sci. Rev.*, 112(1), 67–96.

Debayle, E., and Y. Ricard (2012), A global shear velocity model of the upper mantle from fundamental and higher Rayleigh mode measurements, *J. Geophys. Res.*, 117, B10308, doi:10.1029/2012JB009288.

Debayle, E., and M. Sambridge (2004), Inversion of massive surface wave data sets: Model construction and resolution assessment, *J. Geophys. Res.*, 109, B02316, doi:10.1029/2003JB002652.

Dziewonski, A., and D. Anderson (1981), Preliminary reference Earth model, *Phys. Earth Planet. Inter.*, 25(362), 297–356.

Ekström, G., M. Nettles, and A. Dziewoński (2012), The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes, *Phys. Earth Planet. Inter.*, 200, 1–9.

Endrun, B., S. Lebedev, T. Meier, C. Tirel, and W. Friederich (2011), Complex layered deformation within the Aegean crust and mantle revealed by seismic anisotropy, *Nat. Geosci.*, 4(3), 203–207.

Faccenna, C., and T. W. Becker (2010), Shaping mobile belts by small-scale convection, *Nature*, 465(7298), 602–605.

Fichtner, A., E. Saygin, T. Taymaz, P. Cupillard, Y. Capdeville, and J. Trampert (2013), The deep structure of the North Anatolian fault zone, *Earth Planet. Sci. Lett.*, 373, 109–117.

Gilbert, F. (1971), Excitation of the normal modes of the earth by earthquake sources, *Geophys. J. Int.*, 22(2), 223–226.

Göğüş, O. H., and R. N. Pysklywec (2008), Mantle lithosphere delamination driving plateau uplift and synconvergent extension in eastern anatolia, *Geology*, 36(9), 723–726.

Goldstein, P., D. Dodge, M. Firpo, and L. Minner (2003), Sac2000: Signal processing and analysis tools for seismologists and engineers, in *The IASPEI International Handbook of Earthquake and Engineering Seismology*, vol. 81, pp. 1613–1620, Academic Press, London, U. K.

Griffin, W., S. Y. Oreilly, J. C. Afonso, and G. Begg (2009), The composition and evolution of lithospheric mantle: A re-evaluation and its tec- tonic implications, *J. Petrol.*, 50(7), 1185–1204.

Griot, D.-A., J.-P. Montagner, and P. Tapponnier (1998), Phase velocity structure from rayleigh and love waves in tibet and its neighboring regions, *J. Geophys. Res.*, 103(B9), 21,25–21,232.

Haned, A., E. Stutzmann, M. Schimmel, S. Kiselev, A. Davaille, and A. Yelles-Chaouche (2016), Global tomography using seismic hum, *Geo-phys. J. Int.*, 204(2), 1222–1236.

Herrmann, R., C. Ammon, and K. Koper (2004), GSAC-Generic Seismic Application Computing, *Eos Trans. AGU*, 85(47), Fall Meet. Suppl., Abstract #S21A-0273.

Janutyte, I., M. Majdanski, P. Voss, E. Kozlovskaya, and PASSEQ Working Group (2014), Upper mantle structure around the Trans-European Suture Zone obtained by teleseismic tomography, *Solid Earth Discuss.*, 6, 1723–1763.

Kind, R., T. Eken, F. Tilmann, F. Sodoudi, T. Taymaz, F. Bulut, X. Yuan, B. Can, and F. Schneider (2015), Thickness of the lithosphere beneath Turkey and surroundings from S-receiver functions, *Solid Earth Discuss.*, 7, 1315–1346.

Laske, G., G. Masters, and C. Reif (2001), CRUST2.0: A new global crustal model at 2 × 2 degrees, Inst. of Geophys. and Planet. Phys., Univ.y of Calif., San Diego, Calif. [Available at http://mahi.ucsd.edu/Gabi/rem.dir/crust/crust2.html.]

Laske, G., G. Masters, Z. Ma, and M. Pasyanos (2013), Update on CRUST1. 0-A 1-degree global model of Earths crust, Geophys. Res. Abstracts, 15, 2658.

Lebedev, S., G. Nolet, T. Meier, and R. D. Van Der Hilst (2005), Automated multimode inversion of surface and S waveforms, *Geophys. J. Int.*, 162(3), 951–964.

Legendre, C., T. Meier, S. Lebedev, W. Friederich, and L. Viereck-Götte (2012), A shear wave velocity model of the European upper mantle from automated inversion of seismic shear and surface waveforms, *Geophys. J. Int.*, 191(1), 282–304.

Lekić, V., and B. Romanowicz (2011), Inferring upper-mantle structure by full waveform tomography with the spectral element method, *Geophys. J. Int.*, 185(2), 799–831.

Montagner, J. (1986), Regional three-dimensional structures using long-period surface waves, *Ann. Geophys., Ser. B*, 4(3), 283–294.

Montagner, J.-P., and L. Guillot (2002), Seismic anisotropy and global geodynamics, *Rev. Mineral. Geochem.*, 51(1), 353–385.

Montagner, J.-P., and N. Jobert (1988), Vectorial tomography: II. Application to the Indian Ocean, *Geophys. J. Int.*, 94(2), 309–344.

Montagner, J.-P., and H.-C. Nataf (1986), A simple method for inverting the azimuthal anisotropy of surface waves, *J. Geophys. Res.*, 91(B1), 511–520.

Plomerová, J. (1997), Seismic anisotropy in tomographic studies of the upper mantle beneath Southern Europe, *40(1).*

Plomerová, J., L. Vecsey, V. Babuška, and LAPNET Working Group (2011), Domains of Archean mantle lithosphere deciphered by seismic anisotropy—initial results from the LAPNET array in northern Fennoscandia, *Solid Earth Discuss.*, 3(2), 655–678.

Ritzwoller, M. H., and A. L. Levshin (1998), Eurasian surface wave tomography: Group velocities, *J. Geophys. Res.*, 103(B3), 4839–4878.

Schaefer, J., L. Boschi, and E. Kissling (2011), Adaptively parametrized surface wave tomography: Methodology and a new model of the European upper mantle, *Geophys. J. Int.*, 186(3), 1431–1453.

Scheingraber, C., K. Hosseini, R. Barsch, and K. Sigloch (2013), Obs-load: A tool for fully automated retrieval of seismological waveform data, *Seismol. Res. Lett.*, 84(3), 525–531.

Schivardi, R., and A. Morelli (2009), Surface wave tomography in the European and Mediterranean region, *Geophys. J. Int.*, 177(3), 1050–1066.

Schivardi, R., and A. Morelli (2011), EPmantle: A 3-D transversely isotropic model of the upper mantle under the European Plate, *Geophys. J. Int.*, 185(1), 469–484.

Selway, K., H. Ford, and P. Kelemen (2015), The seismic mid-lithosphere discontinuity, *Earth Planet. Sci. Lett.*, 414, 45–57.

Soomro, R., C. Weidle, L. Cristiano, S. Lebedev, T. Meier, and P. W. Group (2016), Phase velocities of Rayleigh and Love waves in central and northern Europe from automated, broad-band, interstation measurements, *Geophys. J. Int.*, 204(1), 517–534, doi:10.1093/gji/ggv462.

Takeuchi, H., and M. Saito (1972), Seismic surface waves, *Methods Comput. Phys.*, 11, 217–295.

Vecsey, L., J. Plomerová, and V. Babuška (2014), Mantle lithosphere transition from the East European Craton to the Variscan Bohemian Massif imaged by shear-wave splitting, *Solid Earth*, 5(2), 779–792.

Wessel, P., and W. H. Smith (1998), New, improved version of Generic Mapping Tools released, *Eos Trans. AGU*, 79(47), 579–579.

Woodhouse, J. H., and A. M. Dziewonski (1984), Mapping the upper mantle: Three-dimensional modeling of earth structure by inversion of seismic waveforms, *J. Geophys. Res.*, 89(B7), 5953–5986.

Yang, Y., M. H. Ritzwoller, A. L. Levshin, and N. M. Shapiro (2007), Ambient noise Rayleigh wave tomography across Europe, *Geophys. J. Int.*, 168(1), 259–274, doi:10.1111/j.1365-246X.2006.03203.x.

Yuan, H., B. Romanowicz, K. M. Fischer, and D. Abt (2011), 3-D shear wave radially and azimuthally anisotropic velocity model of the North American upper mantle, *Geophys. J. Int.*, 184(3), 1237–1260.
<figure><img src="image_22.png" /></figure>

Zhu, H., and J. Tromp (2013), Mapping tectonic deformation in the crust and upper mantle beneath Europe and the North Atlantic Ocean, *Science*, **341**(6148), 871–875, doi:10.1126/science.1241335.

Zhu, H., E. Bozdağ, D. Peter, and J. Tromp (2012), Structure of the European upper mantle revealed by adjoint tomography, *Nat. Geosci.*, **5**(7), 493–498.

Zhu, H., E. Bozdağ, and J. Tromp (2015), Seismic structure of the European upper mantle based on adjoint tomography, *Geophys. J. Int.*, **201**(1), 18–52.