Traces of the crustal units and the upper-mantle structure in the southwestern part of the East European Craton

The presented study is a part of the passive seismic experiment PASSEQ 2006–2008, which took place around the Trans-European Suture Zone (TESZ) from May 2006 to June 2008. The data set of 4195 manually picked arrivals of teleseismicP waves of 101 earthquakes (EQs) recorded in the seismic stations deployed to the east of the TESZ was inverted using the non-linear teleseismic tomography algorithm TELINV. Two 3-D crustal models were used to estimate the crustal travel time (TT) corrections. As a result, we obtain a model of P -wave velocity variations in the upper mantle beneath the TESZ and the East European Craton (EEC). In the study area beneath the craton, we observe up to 3 % higher and beneath the TESZ about 2–3 % lower seismic velocities compared to the IASP91 velocity model. We find the seismic lithosphere–asthenosphere boundary (LAB) beneath the TESZ at a depth of about 180 km, while we observe no seismic LAB beneath the EEC. The inversion results obtained with the real and the synthetic data sets indicate a ramp shape of the LAB in the northern TESZ, where we observe values of seismic velocities close to those of the craton down to about 150 km. The lithosphere thickness in the EEC increases going from the TESZ to the NE from about 180 km beneath Poland to 300 km or more beneath Lithuania. Moreover, in western Lithuania we find an indication of an uppermantle dome. In our results, the crustal units are not well resolved. There are no clear indications of the features in the upper mantle which could be related to the crustal units in the study area. On the other hand, at a depth of 120–150 km we indicate a trace of a boundary of proposed palaeosubduction zone between the East Lithuanian Domain (EL) and the West Lithuanian Granulite Domain (WLG). Also, in our results, we may have identified two anorogenic granitoid plutons.


Introduction
The East European Craton (EEC) (Fig. 1), the palaeocontinent Baltica, has not been tectonically reworked for at least 1.45 Ga (Bogdanova et al., 2006).The EEC includes a mosaic of tectonic structures.It has formed during the collision of three palaeocontinents: Sarmatia, Volgo-Uralia and Fennoscandia 2-1.7 Ga (Bogdanova et al., 2001;Artemieva, 2007).The EEC in the east is bordered by the Uralides orogen and the Timan Ridge, and in the west by the Trans-European Suture Zone (TESZ), the boundary between Proterozoic Eastern Europe and Phanerozoic western-central Europe (Nolet and Zielhuis, 1994).The inner major sutures in the EEC are the Central Russia Rift System and the Pachelma Rift, which mark amalgation of Baltica in the north, Sarmatia in the west and Volgo-Uralia in the east during the Proterozoic Period (Gorbatschev and Bogdanova, 1993).During a long evolution, the EEC resulted in a complex structure of the crust and the upper mantle, which were intensively investigated during a number of studies (e.g.Artemieva et al., 2006).Guterch et al., 1999Guterch et al., , 2004;;Grad et al., 2006;EUROBRIDGE Seismic Working Group, 1999;Pharaoh et al., 2000;Wilde-Piórko et al., 2010;Knapmeyer-Endrun et al., 2013).
Our study is focused on the SW part of the EEC.The study area covers the NW part of the territory of the passive seismic experiment PASSEQ 2006-2008(Wilde-Piórko et al., 2008), which was carried out around the TESZ in order to study the lithosphere and asthenosphere beneath the area.The aims of our study are to define (1) whether there is a correlation between the crustal units and the upper mantle, and (2) to estimate the seismic P -wave velocity structure of the upper mantle and the lithosphere thickness beneath the study area using the data acquired during the PASSEQ 2006-2008 project and the method of non-linear teleseismic tomography.

Crust and lithosphere structure
The deep seismic sounding (DSS) projects -such as EURO-BRIDGE (EUROBRIDGE Seismic Working Group, 1999), POLONAISE'97 (Guterch et al., 1999), CELEBRATION 2000 (Malinowski et al., 2008), etc. (Fig. 2), carried out around the TESZ in the SW part of the EEC -provided crucial information about the crustal and upper-mantle structure in the area to the depth of about 80 km.The structure of the upper mantle extending to several hundreds of kilometres has been modelled during other studies (e.g.Artemieva et al., 2006;Majorowicz et al., 2003).

Crustal units in Lithuania
The NE part of the EEC is composed of several Svecofennian crustal domains (Figs. 2,3).Grad et al. (2006) and Mo-  Motuza (2005) also interpreted the rocks of the crystalline crust of the WLG as a back-arc complex, rocks of the Mid-dle Lithuania Suture Zone as a volcanic island arc complex, and the rocks of the EL as an accretionary complex.The contact between the EL and the BLG further to the NE is not so prominent (Motuza, 2005).In the WLG the seismic velocities in the uppermost mantle vary from 8.65 to 8.9 km/s and increase along the EUROBRIDGE profile from the west to the east (Motuza et al., 2000).The crustal features of the EL show lineaments extending the NE-SW direction, which coincide with the direction of collision with Sarmatian palaeocontinent (Motuza, 2005;Bogdanova et al., 2001).
The anorogenic magmatism took place around Lithuania and the adjacent areas 1.6-1.5 Ga, resulting in a number of granitoid intrusions (Bogdanova et al., 2006).Two large granitoid bodies of rapakivi-type are present in our study area: the Riga Pluton in western Latvia, and the Mazury Complex in the Kaliningrad District of Russia and NE Poland (Rämö et al., 1996;Dörr et al., 2002).

Crustal units in Belarus
The junction between Fennoscandia and Sarmatia is significant in Belarus (e.g.Bogdanova et al., 1996) (Figs. 2, 3).The crustal pattern in the area shows crustal units with alternating granulite and amphibolite facies which vary in age and origin.The structural features suggest that the accretion was driven by several events of subduction and collision, and the accretionary tectonics prevailed 2.0-1.8Ga (Bogdanova, 1999;Claesson at al., 2001).
At the edge of Sarmatia, there are the Central Belarus Belt (CB) and the Vitebsk Granulite Domain (VG) of the Palaeoproterozoic age (about 2.0 Ga).The VG adjoins the CB in the east and NE.Bogdanova et al. (1996) and Stephenson et al. (1996) indicated the complex crustal structures along the Fennoscandia-Sarmatia junction with the VG and the CB slightly dipping to the SE direction beneath the edge of Sarmatia.The CB consists of bodies of amphibolite and granulite facies (Bogdanova et al., 2001) with significant tectonic faults separating the units of different composition.The study of Claesson et al. (2001) showed that the subcrustal rocks of the VG are similar to those of the southeastern CB.

Crustal units in Poland
The aforementioned (see Sects.2.1 and 2.3) crustal units of prolonged shape (or "belts") from Lithuania and Belarus con-tinue in the SW direction into Polish territory (Bogdanova et al., 2006) and terminate at the TESZ (Fig. 2).The results obtained during the POLONAISE'97 (Guterch et al., 1999) and CELEBRATION 2000 DSS projects provided detailed models of the crust and the upper-mantle structure in Poland (e.g.Czuba et al., 2001;Malinowski et al., 2008).The westernmost part of the EEC adjoining the TESZ has thick continental crust of average thickness of 40-50 km (Grad et al., 2006;Guterch et al., 2004).Dadlez et al. (2005) and Grad et al. (2006) discussed in details the structure of the crust and the uppermost mantle in the SW part of the EEC, obtained by different DSS projects (Fig. 3).Some steep changes in the Moho depths and the seismic velocities along some of the profiles were reported.For example, a "step" with the increase in the Moho depth from 42 to 44 km (which is comparable to the resolution of the method) was found at the P5 profile between the Mazury Complex and the Mazowsze Massif (Czuba et al., 2001).Dadlez et al. (2005) summarized that not all Moho "steps" occur exactly at the places of the proposed terrain boundaries.Moreover, no clear boundaries are visible in the crust between Precambrian terrains postulated by Bogdanova et al. (1996).

Upper mantle structure in the study area
The cratonic lithosphere has been shown to extend to depths of about 200-250 km (Plomerova et al., 2002;Eaton et al., 2009), which is deeper than that of the younger continental regions (e.g.Shomali et al., 2006;Gregersen et al., 2010).Artemieva (2003) found thickness of the thermal lithosphere of about 250-275 km in the EEC for the Archean Kola-Karelian province and some parts of Volgo-Uralia.However, the seismic lithosphere is systematically thicker by about 50 km than the thermal lithosphere (Artemieva, 2007).Sandoval et al. (2004) indicated the high-velocity anomaly extending to a depth of at least 250 km beneath the central part of the Fennoscandian Shield, using the method of bodywave tomography.Hjelt et al. (2006) also reported that, in the Fennoscandian Shield the seismic velocity anomalies extend to the depths of at least 250-300 km.The study of Artemieva et al. (2006) showed the thickness of the thermal lithosphere at about 180 km for the EEC, while the results of geothermal modelling obtained by Majorowicz et al. (2003) indicated thermal lithosphere thickness of 200 km for the EEC.The study of Artemieva et al. (2006) showed thickness of the seismic lithosphere more than 250 km, while Koulakov et al. (2009) observed high P -wave velocities down to 300 km, but Geissler et al. (2010) found no clear seismic LAB beneath the SW part of the EEC.
The reflectors in the upper mantle just beneath the Moho boundary in Fennoscandia were found by Czuba et al. (2001), Yliniemi et al. (2004) and Grad et al. (2002).A major southwards dipping reflector was found beneath the EUROBRIDGE'97 profile, extending from the Moho boundary down to the depth of about 75 km (Thybo et al.,  2003), while a steep SW-dipping mantle reflector reported was below the OMIB, and the VB correlates with a subhorizontal reflector in the EUROBRIDGE'96 profile.Similar subhorizontal lithospheric reflectors were observed beneath the TESZ (Grad et al., 2002;Guterch et al., 2004) and the Baltic Sea (Hansen and Balling, 2004).Beneath the WLGD in the upper mantle reflectors at a depth of 73-82 km were reported, which possibly originated due to delamination processes (Motuza et al., 2000;Motuza, 2005).Moreover, a locally increased heat flow ranging between 55 and 100 mWm −2 was found in the WLGD (Kepezinskas et al., 1996;Rasteniene et al., 1998).

Data set
We used some of the data recorded during the PASSEQ 2006-2008project (Wilde-Piórko et al., 2008), which took place around the TESZ from June 2006 to July 2008 (Fig. 4).Using seismological bulletins of the USGS (http: //earthquake.usgs.gov/)and the ISC (http://www.isc.ac.uk/), we prepared a list of 101 earthquakes (EQs) with epicentral distances from 30 to 92 degrees (Artlitt, 1999;Sandoval, 2002), with respect to the central point at the Lithuanian-Polish border (coordinates 23 • E and 54 • N) and the magnitude range from 5.5 to 7.2 (Fig. 5).The higher and lower values of the epicentral distance ensure that the first-observed arrivals are the direct P waves, and that they hit the target area steeply enough from below.The relatively large magnitudes ensure better quality of the observed seismic signals.
On the other hand, the magnitudes should not be too large, because it is difficult to interpret the seismic signals generated by large-scale seismic sources.We used the Seismic Handler Motif (SHM) program package (http://www.seismic-handler.org/) to perform the analysis and manual picking of the teleseismic P -wave arrivals.During the data analysis, we applied the World Wide Standardized Seismographic Network short period (WWSS-SP) filter, which includes both simulation filtering and instrument response, and picked the P -wave arrivals on seismograms of vertical components (Fig. 6).Every P -wave arrival was assigned with a quality (weighting) factor depending on the time (or picking) error (Table 1).The weighting factor was taken into account during the inversion.We compiled a data set of 4195 P -wave arrivals from the data of 94 seismic stations deployed to the east of the TESZ.
We used Seismic Handler (SH) program package and location information of the listed 101 EQs from the ISC seismological bulletins to calculate the theoretical travel times (TTs) of the first teleseismic P -wave arrivals.Then we applied a subtraction procedure in order to obtain the TT residuals for Out of all the seismograms we picked the best trace (the reference station) with relatively high signal-to-noise ratio (SNR) and picked the absolute P -wave arrival (P_abs) (the onset of the P wave) and the relative P -wave arrival (P_ref) of some well-expressed minima or maxima of the seismic signal on the same trace (i.e.station PP81).Then we compared the waveform of the reference seismogram with the waveforms of other seismograms, and picked the relative P -wave arrivals there.For some EQs, we observed more than one type of the waveform.Thus, we grouped the events with similar waveforms and picked absolute and relative P -wave arrivals for each group separately.Every pick was assigned with a quality factor according to the picking error from 1 (best quality) to 3 (poor quality).The purple picks (stations PP81, PA81, PB50, PD83 and PJ42) indicate P -wave arrivals of quality factor 1, while the red ones (station PF47) indicate P -wave arrivals of lower quality (either 2 or 3).In the data of some stations, we indicated inverted polarities (station PD83).every picked arrival: where T picked is the observed TT, T theoretical is the theoretical TT calculated with SH and T residual is the TT residual.

Teleseismic tomography inversion
We used TELINV code (Voss et al., 2006) to perform inversion of the compiled data set.The program utilizes a nonlinear inversion method and can either (1) calculate propagation of rays through a 3-D velocity model and output TT, raypaths and synthetic relative TT, or (2) invert teleseismic relative P -wave residuals for 3-D velocity structure.
The ray tracing is performed by computing the 3-D minimum TT raypaths, assuming a constant slowness in each cell (Steck and Prothero, 1991).The ray coverage of the cell blocks is affected by horizontal and vertical grid spacing (Arlitt, 1999).For the full description of the inversion procedure, see Thomson andGubbins (1982), Thuber (1983), Menke (1984), Koch (1985) and Aki et al. (1997).

Crustal travel time corrections
In teleseismic tomography, it is very important to use reliable crustal TT corrections in order to eliminate the effects which are created by the earth's crust, while the crust is much more heterogeneous compared to the deeper layers of the earth.
The variation of thickness of the sedimentary cover is significant in the study area ranging from several tens of metres in the Belarus-Mazurian High to almost 20 km in the Polish Basin, while the Moho variation is from 35 km beneath the TESZ to almost 60 km beneath NW Belarus.The teleseismic tomography inversions performed without (Fig. 7a) and with (Fig. 7b) crustal TT corrections show relatively similar distribution of the high and low velocity areas, but the seismic velocity contrast in total is up to about 4 % higher in the results obtained without crustal corrections.The largest differences are observed on the eastern edge of the TESZ beneath the Polish Basin and beneath western Lithuania, where significant sedimentary covers up to 20 and 2 km thick, respectively, are present.
The crustal TT corrections which we use in our study have been compiled using two 3-D crustal models by Majdański (2012) for Poland (Fig. 8a left) and by M. Budraitis (unpublished) for Lithuania (Fig. 8a right).Both models have been compiled using results of available DSS projects (e.g.EU-ROBRIDGE, CELEBRATION, POLONAISE, BABEL, Sovietsk -Kohtla-Järve, etc.) carried out around Poland and Lithuania.We calculated the crustal TT corrections using the following equation: where TT model is TT through the crustal velocity models by Majdański (2012) or by M. Budraitis (unpublished), TT iasp is TT through the IASP91 velocity model and TT diff is TT difference.Although the crustal TT corrections for individual seismic stations do not take into account the bending of the seismic rays in the crust, the result is reliable as the rays hit the surface almost vertically and the crust is thin compared to the entire velocity structure.

Model parameterization
In the teleseismic tomography inversion as an input model, we used the 1-D IASP91 velocity model (Kennett and Engdahl, 1991) and transformed it into the 3-D velocity model with 16 layers of different thicknesses down to 700 km.As the resolution of the inversion is governed by spacing between seismic stations, frequency content of the seismic signals and seismic ray geometry, we used spacing of 50 km between the nodes of the model grid in horizontal directions (Fig. 4).We performed a number of inversions with different values of smoothing and damping in order to assess the optimal parameters of the inversion.After careful analysis, we found that the same value as spacing between the grid nodes in horizontal directions (i.e.50) is applicable for the diagonal elements of the smoothing matrix, while the optimal damping value 80 was determined investigating the trade-off curve between the data variance and model variance (Fig. 9).The inversions with both the synthetic and the real data sets were performed using the defined optimal parameters of smoothing and damping for 10 layers between 60 and 350 km depths.

Resolution
The resolution assessment includes calculation of spatial resolution and standard deviations of the model parameters and helps to evaluate the precision of inversion results.In our study, we used the hit matrix method to assess the resolution, and the synthetic checkerboard test with the real station configuration in order to indicate the parts of the study area which could be reasonably resolved.The hit matrix is based on a calculation of the number of rays which transverse a particular cell (Fig. 10).The compiled synthetic checkerboard velocity model contains 200 km wide blocks in the horizontal directions and four layers thick with ±4 % velocity difference compared to the IASP91 velocity model (Fig. 11a).The inversion results show that the synthetic velocity structure is fairly well resolved in horizontal directions in the areas with good station coverage, while the vertical smearing is quite significant (Fig. 11b).The W-E smearing dipping to the east is most likely due to the seismic rays coming mainly from the NE-E-SE, because we use more EQs from the region to the east (i.e.Japan, Kamchatka, Sumatra and Aleutian regions) due to higher seismic activity compared to the region to the west of the study area.

Synthetic "geological" model
We compiled a synthetic "geological" 3-D velocity model using the velocity model by Wilde-Piórko et al. (2010) as a base, but we modified both thicknesses of different layers (because of different model grid) and some values of the seismic velocities regarding some other studies (e.g.Griffin et al., 2003).In our synthetic model (Fig. 12a), we introduced the seismic P -wave velocities 2-6 % higher compared to the IASP91 velocity model at different depths beneath the craton.In the TESZ area, we introduced the shape of a ramp-type of the lithosphere-asthenosphere boundary (LAB) dipping to the NE with seismic velocity values close to those of the cratonic part but up to 2 % smaller in the upper layers down to about 180 km.At the depths between 270 and 350 km, we introduced velocities 2 to 4 % lower compared to the IASP91 velocity model, and the higher velocity area (about 3 % higher P -wave velocities compared to the IASP91 velocity model) in the NE part, which implies that we expect the deeper cratonic roots in this part of the study area.
With the synthetic data set, we performed inversions without (Fig. 12b) and with (Fig. 12c) the crustal TT corrections, as used with the field data.The crustal corrections were applied in order to obtain similar raypaths in the upper layers, and to estimate the effects of the crustal corrections to the signal amplitudes and the depth to which this effect is significant.The inversion results with the crustal TT corrections (Fig. 12c) show in total about 2.5 % higher signal amplitudes (both positive and negative) compared to the results obtained without crustal corrections (Fig. 12b).This high value of signal amplitudes is caused because the synthetic data set was compiled using theoretical TT only, while the crustal corrections used with the field data reduces the signal amplitudes (Fig. 7).Moreover, we indicate that the effect due to the crustal TT corrections is significant (up to 0.5 %) down to about 120 km, while, going deeper, the effect is negligible (Fig. 12b, c).Both results obtained using the synthetic data set with and without crustal TT corrections (Fig. 12b, c) show reasonably resolved ramp-type shape of the LAB and the deep cratonic roots going down to 350 km in the NE part of the study area.

Results and discussion
In teleseismic tomography, many factors -which are related to either the model parameterization or the field data setinfluence the observed signal amplitudes: (1) in teleseismic tomography only a part of the ray path through the velocity model is inverted.However, the rays experience distortions along their full paths from source to receiver (i.e.outside the velocity model) which are mapped in the final results of the inversion, and add up to positive or negative signals; (2) the TELINV code used in this study implements the "flat-earth" transformation which has an effect on the apparent velocities when dealing with large study areas.Our study area is 800 km in the longest direction and the model is set to the depth of 700 km.Thus, the model at the bottom is horizontally stretched by about 11 %.Regarding the incidence angles of our data set, the effect on velocity perturbations due to the flat-earth transformation is about 1.5 % of the observed velocity contrast.However, we did not apply any corrections for the spherical model; (3) the damping value has an influence on the velocity contrast -the larger the value, the smaller the velocity contrast.On the other hand, too high a damping value would result in reduction of lateral variations.The optimal damping value (which is 80) used in the inversion was set after careful analysis (Fig. 9); (4) the precise crustal TT corrections are essential in teleseismic tomography in order to eliminate (or reduce) the crustal effects.We applied the crustal TT corrections assuming vertical ray  propagation in the crust.However, its effect to signal amplitudes is negligible.The larger effect (up to 1 %) can be caused by the used crustal models because they have their own limit of precision; (5) the increase in signal contrast in the results due to temperature variations could be up to about 1 % because, in the study, the heat-flow variations are significant; (6) some anisotropy studies of share-wave splitting (e.g.Wüstefeld et al., 2010;Vecsey et al., 2013;Sroda et al., 2014) show relatively small anisotropy for the EEC compared to the territories to the west of the TESZ.Thus, its effects to the observed velocity contrast should be quite small (up to about 0.5 %).Taking into account all the above-listed causes, we should consider the velocity contrast not up to ±6 %, which we observe in our results (Figs. 12d,13), but close to ±3 %.

Lithosphere structure
In our study area, we resolved structure of the upper mantle from 60 km down to 350 km (Figs. 12d,13).Beneath the EEC (Fig. 12c) we obtain up to 3 % higher seismic velocities compared to the IASP91 velocity model.The higher velocities in the upper mantle can be traced going down to the depth of about 180 km beneath NE Poland which coincides with the results of Wilde-Piórko et al. (2010) and Majorowicz et al. (2003).Going further to the NE the lithosphere thickness increases and beneath Lithuania it is at least 300 km or more, which coincides well with observations by Koulakov et al. (2009).Thick lithosphere was previously reported for other cratonic areas, i.e. the Fennoscandian Shield (Sandoval et al., 2004), but no evidence of the seismic LAB was found anywhere within the depth of 300 km beneath the Fennoscandian Shield (Bruneton et al., 2004).The shear-wave studies of Legendre et al. (2012) showed no deep cratonic roots below about 330 km in the EEC.There is a good correlation between our results obtained with the real data set and the synthetic data set (Fig. 12), which implies that the lithosphere thickness may increase going from the TESZ towards the NE and could be larger than 300 km in the EEC.Moreover, the results with the synthetic data set show that the P -wave velocities beneath the craton down to 180 km could be about 3 % higher compared to the IASP91 velocity model.
In the EEC beneath western Lithuania (i.e. the WLG) down to 90 km, we observe the lower seismic velocities (Fig. 13) which could be related to an upper-mantle dome.Motuza et al. (2000) proposed that the mantle dome could be related with delamination processes because, beneath the WLG, the heat flow, which is significantly higher compared to the adjacent areas, was observed (Kepezinskas et al., 1996;Rasteniene et al., 1998) and the high-density reflectors in the upper mantle have been found (Giese, 1998;Motuza et al. 2000).These high-density bodies can potentially represent delaminated slices of the crust which sank into the mantle (e.g.Defant and Kepezhinskas, 2002).In our results (Figs. 12c, 13), we do not find any well-defined high velocity reflector in the upper mantle.On the other hand, below the discussed low velocity area (i.e. the proposed upper-mantle dome), we observe area of velocities which are significantly higher than those of the surroundings.As the delamination processes occur locally, the lower and the higher velocity areas observed in our results beneath the WLG could possibly be related to the local upper-mantle dome and the delaminated high-density rocks.
Beneath the TESZ, we find about 2 to 3 % smaller seismic velocities compared to the IASP91 velocity model, except for the northern TESZ (northern Poland), where we observe the values of seismic velocities close to those of the craton down to about 150 km.In their work, Knapmeyer-Endrun et al. ( 2013) observe an increase in TT of Ps conversions across the mantle transition zone which they think could be caused either by a temperature reduction or an increase in water content.
In the northern part of the TESZ, we find the seismic LAB at a depth of about 180 km (Fig. 12d).We also indicate the seismic LAB of a ramp-type dipping towards the NE, which coincides with the inversion results obtained with the syn-thetic data set (Fig. 12).Hansen and Balling (2004) also reported on a number of mantle reflectors beneath the Baltic Sea along the TESZ dipping towards the N-NE.
The velocity model by Wilde-Piórko et al. (2010) proposed the higher P -wave velocity values compared to the IASP91 velocity model for the TESZ and the cratonic area for depths more than 250 km.Our results with the real and the synthetic data sets (Fig. 12) indicate that the seismic velocities at these depths are 1-2 % smaller or similar compared to the IASP91 velocity model.

Traces of the crustal units
The crustal units are not well resolved in our results.There are no clear indications of the structures (Fig. 2) in the upper mantle (the uppermost inverted layers of the velocity model) which could be related with the crustal units in the study area (Figs. 12c,13).However, in the uppermost inverted layers we find correlation between the Moho depth and velocity variations: the positive signal amplitudes are usually observed in the areas with thicker continental crust beneath Poland and Lithuania, while the negative ones are in the areas with thinner crust beneath the TESZ.This could be related either to the imperfect crustal TT corrections used or to different geological conditions.We may infer only one possibly resolved boundary between the EL and the WLG beneath Lithuania, which could be related to the local lower velocity areas at the depths from 120 km to at least 150 km.This area was interpreted by Motuza (2005) and Motuza and Staškus (2009) as a palaeosubduction zone.The other possible explanation for the lower velocity area beneath the southernmost region of Lithuania -NE Poland at a depth of 100-120 km -is an effect due to an anorogenic granitoid massif, the Mazury Complex (Fig. 7), which is 40 km wide and 6.5 km thick, extending 200 km from the Baltic Sea through the Kaliningrad District of Russia into NE Poland.A number of studies (e.g.Bruneton et al., 2004;Beller et al., 2013) showed that the upper mantle beneath anorogenic granitoid massifs inside cratonic crust is different from that of the surrounding cratonic mantle.There is another anorogenic granitoid massif, the Riga Pluton (Fig. 7), in western Latvia which, in our results, could be related to the lower velocity area down to about 150 km in the NE region of our study area.As the granitoid massif lies on the edge of the study area where resolution is quite poor, we cannot assert its effects on our results.Both the Mazury Complex and the Riga Pluton are of rapakivi type and have formed 1.6-1.5 Ga (Rämö et al., 1996;Dörr et al., 2002).

Conclusions
-Beneath the EEC, we obtain up to 3 % higher seismic velocities compared to the IASP91 velocity model.The lithosphere thickness increases towards the NE from about 180 km beneath NE Poland to at least 300 km or more beneath Lithuania.
-Beneath the TESZ, we find the seismic velocities 2-3 % smaller compared to the IASP91 velocity model, and only in the northern TESZ do we observe higher seismic velocities down to about 150 km, which show that the northern part of the TESZ is more craton-like.
-The seismic LAB beneath the northern part of the TESZ is at a depth of about 180 km, and it is most likely of the shape of a ramp.We did not find the seismic LAB beneath the EEC.
-The seismic velocities in our study area at the depths more than 250 km could be 1-2 % smaller or similar compared to the IASP91 velocity model.
-The observed local lower and higher velocities beneath western Lithuania might be related to an upper-mantle dome.
-In our results, we did not find a strong correlation between the separate crustal units and the upper mantle.However, we observe that the positive signal amplitudes coincide with the areas with thicker continental crust, while the negative ones coincide with areas with thinner crust.We also recognize the trace of the palaeosubduction boundary between the EL and the WLG beneath Lithuania.
-In our results, we possibly identified the Riga and the Mazury anorogenic granitoid plutons.

Figure 3 .
Figure 3. Models of the crust and the uppermost mantle along the EUROBRIDGE transect (EB'94 and EB'96), the POLONAISE'97 profiles P4 (northern part), P5 and P3, and CELEBRATION 2000 profile CEL05 (after Grad et al., 2006).Values of the P -wave velocities are given in kilometres per second.Arrows indicate positions of the shot points; the crossing points with other profiles are marked in blue.For other explanations see Fig.2.

Figure 4 .
Figure 4. Map of the seismic stations (triangles) used in the study, and locations of nodes of the model grid (dots).The area in between the dashed lines indicates the TESZ.

Figure 5 .
Figure 5. Map of the epicentres of 101 EQs used in teleseismic tomography inversion.Grey rectangle indicates the study area.

Figure 6 .
Figure 6.Example of manual picking of the P -wave arrivals.Filtered seismograms of an EQ on 02.08.2007 at 03:21 UTC.Out of all the seismograms we picked the best trace (the reference station) with relatively high signal-to-noise ratio (SNR) and picked the absolute P -wave arrival (P_abs) (the onset of the P wave) and the relative P -wave arrival (P_ref) of some well-expressed minima or maxima of the seismic signal on the same trace (i.e.station PP81).Then we compared the waveform of the reference seismogram with the waveforms of other seismograms, and picked the relative P -wave arrivals there.For some EQs, we observed more than one type of the waveform.Thus, we grouped the events with similar waveforms and picked absolute and relative P -wave arrivals for each group separately.Every pick was assigned with a quality factor according to the picking error from 1 (best quality) to 3 (poor quality).The purple picks (stations PP81, PA81, PB50, PD83 and PJ42) indicate P -wave arrivals of quality factor 1, while the red ones (station PF47) indicate P -wave arrivals of lower quality (either 2 or 3).In the data of some stations, we indicated inverted polarities (station PD83).

Figure 8 .
Figure 8.(a) Moho maps compiled by Majdański (2012) (left) and M. Budraitis (unpublished) (right) used to estimate the crustal TT corrections.The Moho depths in the depicted area vary from 27 to 57 km.(b) Estimated crustal TT corrections for individual seismic stations.Values are expressed in seconds relative to the IASP91 velocity model.

Figure 9 .
Figure 9. Data variance versus model variance obtained during inversions with damping values from 10 to 360.The optimal damping value was set at 80.

Figure 10 .
Figure 10.Resolution at depth of 90 km for the field data set.Low and high values show, respectively, poorly and well-resolved areas.White triangles mark the seismic stations.

Figure 11 .
Figure 11.Results of the synthetic checkerboard test.We added small random perturbations to the compiled synthetic data set.Horizontal slice at a depth of 90 km and two vertical slices parallel to the main PASSEQ transect of the target area.(a) Initial velocity model with synthetic blocks of 200 km wide in the horizontal directions and ±4 % P -wave velocity difference compared to the IASP91 velocity model.(b) Inversion results with the synthetic data set.Dashed lines indicate the TESZ.Triangles indicate the seismic stations, and on the vertical slices they indicate seismic stations ±50 km around the depicted transects.

Figure 12 .
Figure 12.The initial synthetic "geological" velocity model (a), the inversion results with the synthetic data set with (c) and without (b) the crustal TT corrections, and the inversion results with the real data set with applied crustal TT corrections (d).The P -wave velocity perturbations on the horizontal slices at different depths and the vertical slice along the main PASSEQ transect.The bluish and reddish areas show, respectively, the higher and the lower P -wave velocities compared to the IASP91 velocity model.The thin dashed lines on the horizontal slices indicate the TESZ.Triangles indicate the seismic stations, and on the vertical slices they indicate seismic stations ±50 km around the main PASSEQ transect.The solid thin lines on the horizontal slices (right side) indicate boundaries of different tectonic units (for detailed explanation see Fig.7).Interpreted velocity anomalies on horizontal slices: 1 -upper-mantle dome; 2 -effect from the Riga batholith; 3 -palaeosubduction boundary between the WLG and the EL; 4 -higher velocity anomaly beneath the northern part of the TESZ; and 5 -effect from the Mazury Complex.Solid and dashed lines on vertical slice mark the interpreted LAB.

Figure 13 .
Figure 13.Vertical slices perpendicular to the main PASSEQ transect close to the eastern edge of the TESZ (left), and about 350 km to the NE from the TESZ (right).The thick lines indicate possibly resolved boundary between the EL and the WLG, and the mantle dome beneath the WLG.

Table 1 .
Data set compiled during the manual picking procedure of the P -wave arrivals.