POLENET / LAPNET teleseismic P-wave traveltime tomography model of the upper mantle beneath northern Fennoscandia

Introduction Conclusions References Tables Figures


Introduction
Recently, dense two-dimensional (2-D) networks of broadband seismic instruments have proved to be a most effective means to study the three-dimensional (3-D) structure of the lithosphere (Trampert and Van der Hilst, 2005).One such network was the POLENET/LAPNET broadband seismic network (http://www.oulu.fi/sgo-oty/lapnet)deployed in northern Fennoscandia (Finland, Sweden, Norway, and Russia) during the third International Polar Year 2007-2009.The project was a part of the POLENET (Polar Earth Observing Network, http://www.polenet.org)consortium.The network consisted of 37 temporary and 21 permanent seismic stations (Fig. 1).All the stations, except two temporary stations, were broadband with flat frequency response at least between 0.01 and 25 Hz.The network registered waveforms from teleseismic, regional, and local events during May 2007-September 2009.The average distance between the stations was 70 km.
One of the main targets of POLENET/LAPNET was to obtain a 3-D seismic model of the upper mantle in the northern Fennoscandian Shield, in particular, beneath its Archaean domain (Fig. 1), as the area has not been studied previously by dense broadband seismic networks.Since 1980Since -1990, after the discovery of two large diamond deposits in its northeastern margin (close to the city Arkhangelsk in Russia), the area has been considered to be worth further investigating for diamondiferous kimberlitic rocks.This supposition is based on three empirically established factors necessary for diamond preservation: Archaean bedrock, a low geothermal gradient, and a thick lithosphere (Clifford, 1966).
As shown by Snyder et al. (2004), recordings of teleseismic events can be used to explore the mantle lithosphere to depths of several hundred kilometres, and through geological interpretations, to assess its potential as a diamond reservoir.Generally, such modelling of the upper mantle needs to include P and S wave velocity models estimated by teleseis-  Lahtinen et al. (2015) in northern Fennoscandia.The Fennoscandian crustal block is outlined by a black dotted line and our study area by a rectangle.(b) The simplified geological map is based on the 1 : 2 000 000 geological map of Fennoscandia (Koistinen et al., 2001).The seismic stations are shown with black stars, triangles, and dots, and main geological provinces are named.mic body wave tomography, the position of major boundaries in the crust and upper mantle estimated by controlled source seismics (CSS), and/or receiver function techniques and strength and orientation of seismic anisotropy that can be assessed by shear-wave splitting, receiver function analysis, or ambient noise tomography.Such a combined seismic model makes it possible to evaluate the thickness and composition of the mantle lithosphere beneath the study area.It may further be used for 3-D mapping of lithospheric architecture and structures responsible for the formation, accumulation, and preservation of economically significant mineral deposits (Laznicka, 2014;Mole et al., 2014).

Published by Copernicus
In our study, we explore the 3-D P wave architecture of the upper mantle beneath the northern Fennoscandian Shield using high-resolution teleseismic tomography.This technique was used by Sandoval et al. (2004), who obtained a 3-D P wave velocity model of the upper mantle beneath the south-eastern Fennoscandian Shield, documenting a deep lithospheric "keel" structure beneath that region.A similar technique was also used by Shomali et al. (2006) and Janutyte et al. (2014), in order to study the upper mantle structure in the transition zone between the Precambrian Fennoscandian Shield and the Palaeozoic lithosphere of Western Europe.Eken et al. (2007) used teleseismic tomography with data of the Swedish National Seismic Network in order to estimate the upper mantle structure beneath Sweden.The data of the POLENET/LAPNET experiment were used by Plomerová et al. (2011) and Vinnik et al. (2014), who estimated seismic anisotropy in the upper mantle beneath the POLENET/LAPNET study area.Vinnik et al. (2016) report on variations of S and P wave velocities and Vp/Vs ratio beneath the POLENET/LAPNET network using joint inversion of P and S receiver functions.Our study thus complements the previous studies of the upper mantle in the Fennoscandia by the body-wave tomography technique.It is also a valuable contribution to the combined seismic model of the upper mantle beneath the northern Fennoscandian Shield.

Tectonic setting
Our study area is located in the Fennoscandian Shield in the northern part of the East European Craton.The area consists of the Archaean Karelian craton in the eastern part of the study area, subdivided into Karelian and Kola provinces and the Belomorian Mobile Belt in between, the Svecofennian Norrbotten craton in the western part, and the Caledonides in the north-western corner (Fig. 1b).
The Lapland-Kola orogeny was preceded by subduction of the new oceanic crust and by island arc accretion at 1.95-1.91Ga.The orogeny was a transpressional continentcontinent collision between Kola province and Karelian province and produced only a minimal amount of juvenile crust (Lahtinen et al., 2008).The juvenile material seems to be dominated by Archaean crustal origin though mantle contribution cannot be ruled out (Heilimo et al., 2014).Although in general, the Svecofennian orogeny formed a large unit of new Palaeoproterozoic crust, similarly to the Lapland-Kola orogeny, its northern part is mainly comprised of reworked Archaean material.The orogeny began from the north, where the Karelian craton and Norrbotten Craton collided at ca. 1.92 Ga (Lahtinen et al., 2008(Lahtinen et al., , 2015)).
After these two orogenies, there have been no major tectonic events in our study area, but there is still some smaller volume magmatism with ages from Palaeoproterozoic to Devonian (Downes et al., 2005).These rocks are often from relatively deep sources and, though quite diverse in rock types, are some of the most important carriers of mantle xenoliths (Woodard, 2010).

Data set
The main data set used in this study is the data of the POLENET/LAPNET network that was recorded during May 2007-September 2009 (Kozlovskaya et al., 2007).For this study, we selected 96 teleseismic events that were located at epicentral distances between 30 • and 90 • from our station network and were clearly recorded by most of the stations.In addition to the new data set from the POLENET/LAPNET network, we included previous data from the northern part of the SVEKALAPKO network (Hjelt et al., 2006) located 30°6 0°9 0°F igure 2. Location of the 111 teleseismic events used as sources for the tomography study.The centre of our array is shown with the blue star, while events recorded by the POLENET/LAPNET array are shown with red stars (96 events) and those recorded by the SVEKALAPKO project with yellow stars (15 events).
south of the POLENET/LAPNET network but still within our study area in our study.These data consist of 15 events.Most of the selected events have magnitudes larger than 6.0, but events with magnitudes larger than 5.5 were also considered in order to improve the azimuthal coverage.In this way, we obtained a good coverage over all back azimuths (see Fig. 2); the largest azimuthal gap in events recorded by POLENET/LAPNET network is smaller than 15 • .The first arrivals of P waves were picked manually using Seismic Handler software (Stammler, 1993).We used the Worldwide Standard Seismological Network (WWSSN) short period simulation filter with maximum displacement magnification at approximately 1 s for easier comparison between waveforms recorded by different types of sensors.The travel time residuals were calculated using the iasp91 reference model (Kennett and Engdahl, 1991).While picking the arrivals, the uncertainty of each arrival time was also estimated.The 3167 picked arrival times were then divided into three quality classes based on individual uncertainty estimates.Table 1 shows the error estimates for each class and the number of arrival times attributed to different classes.Example waveforms with arrival time picks of different quality are shown in Fig. 3.
The absolute residuals calculated as differences between observed and theoretical arrival times, using the iasp91 velocity model and ISC catalogue hypocentre parameters, exhibit large amplitudes while showing only small azimuthal variations between stations (see Fig. 4a).This is the combined effect of deep mantle velocity variations located outside our study volume and of near-source structure or  teleseismic hypocentre parameter uncertainties.To separate these effects from our observations, the average travel time residual over all stations was calculated for each event and subtracted from all corresponding residuals.The resulting relative residuals contain the effects of the velocity variations within our study volume with reference to the iasp91 model.Figure 4b shows the azimuthal distributions of the residuals for three selected stations representing different regions within our study area.
The additional data recorded by SVEKALAPKO network included 15 events recorded at 31 stations yielding 360 P wave residual times (Plomerová et al., 2006).To ascertain the compatibility of the two data sets, we compared the residuals recorded at four seismic stations operational during both SVEKALAPKO and POLENET/LAPNET projects and found no significant difference in residual amplitudes or distribution.The comparison figures are available in the Supplement.The quality class distribution of the SVEKALAPKO data is shown in Table 1.All seismic stations included in this study are shown in Fig. 1b and the event distribution is shown in Fig. 2.

The inversion method
In seismic tomography, the basic system of equations relating velocity perturbations inside the study volume to travel time residuals is where d is a data vector composed of travel time residuals, m is a model parameter vector composed of velocity perturbations in each cell of the velocity model, and G is a matrix of partial derivatives defining the coupling between data and model parameters (e.g.Kissling, 1988).The method searches the velocity perturbations inside a defined  3-D model in order to explain the observed travel time residuals.Velocities outside and initial velocities inside the study volumes are approximated using a reference model.The reference model used in this study is iasp91 (Kennett and Engdahl, 1991).
For the inversion of the P wave residual data, we used TELINV, a non-linear tomographic inversion program, originally developed by Evans and Achauer (1993), who applied the ACH tomographic method by Aki et al. (1977) to lithosphere investigations at regional scale, and later modified by several authors (e.g.Weiland et al., 1995;Arlitt et al., 1999;Lippitsch et al., 2003;Sandoval et al., 2004;Karousová, 2013) to include 3-D ray tracing and iterative non-linear inversions.
The initial 3-D model is defined as an orthogonal grid of nodes approximating the study volume, where the velocity is defined at the node points.The travel time calculation is based on the 3-D simplex ray-tracing technique (Steck and Prothero, 1991).The inversion problem is formulated as a weighted damped least square problem and solved using a singular value decomposition technique.The overall non-linear inversion scheme iteratively inverts the travel time residuals for velocity changes relative to the 3-D velocity model of the previous iteration, beginning with the initial reference 3-D model.Each iteration involves a complete onestep inversion, including both ray tracing (forward problem) and inversion for an update of the velocity distribution in the model.
The basic inversion equation for the TELINV code can be written as where m est are estimated model parameters, W D is the weighting matrix of the data, ε 2 is a damping factor, and W M is the smoothing matrix of the model (Menke, 1984).
The capabilities of the ray geometry and model parameter grid to resolve the velocity perturbations can be formally estimated by the resolution matrix (Menke, 1984): (3) The best results for resolution estimates, however, are obtained by synthetic data testing in combination with hit count, derivative weighted sums, and resolution matrix (Kissling, 1988).

Model parameterization and regularization
The POLENET/LAPNET study area is approximately 550 km in the east-west direction and 600 km in the north-south direction.For this study area, we selected an 80 km × 80 km × 60 km inversion grid.The horizontal size of the cells is slightly larger than the average distance between stations to guarantee at least one station for each surface cell within the network space.Thus, the minimum dimension of a resolvable anomalous body beneath the POLENET/LAPNET network is approximately 100 km.The uppermost 60 km is mostly made of crust and is not included in the high-resolution teleseismic inversion (e.g.Arlitt et al., 1999).
In the case of real data, a priori crustal corrections are applied to travel time residuals (see Sect. 4.2).For the inversion of synthetic data (see Sect. 4.3) no crustal correction was applied as the crustal layers were assumed to have velocities according to the iasp91 reference model.The first inverted layer, below the 60 km crustal layer, is 40 km thick, and all deeper layers are 60 km thick.For better performance of the ray tracer, the study volume was surrounded on all sides and at the bottom with large cells of respective iasp91 velocities that were kept fixed during inversion.The P wave velocities of the initial model correspond to the iasp91 reference model.
Because of the uneven ray geometry the kernel matrix (G in Eq. 1) can be singular and needs to be regularized.The regularization can be done by selecting an appropriate damping value for the data set (ε 2 in Eq. 2).We analysed the effect of the damping values by performing multiple inversions with The number of iterations was also optimized (dots with the number of iterations marked to their left side) for the selected damping value of 100.The final number of iterations, 4, is marked with a red circle.The grey line denotes the overall data uncertainty, 0.017 s 2 , calculated as the average uncertainty of all observations.different damping values both with synthetic data and real data to find the best damping value for our data set.The final damping value of 100 was selected by investigating the trade-off curve between model and data variance (see Fig. 5 for the trade-off curve of the real data and Fig. 6 for examples of the crustal corrected residuals).At the same time, the smallest singular value to be used (minSV) was estimated using a similar method that was used with damping.Finally, a value of 70 was selected.

Crustal correction model
Due to steep incidence angles and subsequent near-lack of cross firing, teleseismic rays are largely incapable of resolving upper lithosphere structure such as Moho topography and 3-D crustal velocity variations.Lateral variation of crustal structure, however, may significantly affect travel times of teleseismic rays (Arlitt et al., 1999).Hence, when illuminat- ing the structures at upper mantle depths with high-resolution teleseismic tomography, it is necessary to correct the data a priori for the effect of crustal structures as documented for southern Finland by Sandoval et al. (2003).A new map of the crustal thickness was established by Silvennoinen et al. (2014) based on previous and new controlled source seismic and receiver function results in our study area.This map in combination with an averaged 1-D velocity-depth function is used in this study to construct a 3-D crustal model for the purpose of correcting travel times for crustal effects.
The 3-D crustal model used in this study was established by using Bloxer software (https://wiki.oulu.fi/display/~mpi/Block+model+maintenance).Bloxer is a piece of software designed to build 3-D block models with different parameters (density, seismic velocity, magnetization etc.).The Moho depth map by Silvennoinen et al. (2014) was imported to Bloxer.All blocks below the Moho to a maximum model depth of 60 km were given a P wave velocity value 8.05 km s −1 that is the P wave velocity of the uppermost mantle of the iasp91 reference model.
As our study area does not have enough P wave velocity information available to estimate a true 3-D distribution of seismic velocities in the crust between CSS profiles, we used an average P wave velocity for crustal depths defined from major refraction seismic profiles, namely FEN-NOLORA (Guggisberg, 1986;Guggisberg et al., 1991;Luosto and Korhonen, 1986) in Sweden, POLAR (Luosto et al., 1989;Janik et al., 2009) in Finland, and PECHENGA-KOSTOMUKSHA (Azbel et al., 1989(Azbel et al., , 1991;;Azbel and Ionkis, 1991) in Russia.A map of crustal correction times for vertical incident rays arriving at each station is shown in Fig. 6a.The crustal effect established by Sandoval et al. (2003) for southern and central Finland has an almost perfect fit with the crustal effect established in this study for the overlapping southern part of our study area (south of 65.5 • N).
To evaluate the travel time through this pseudo-3-D model with averaged velocity-depth function and locally variable Moho depths in comparison to the travel time through a local 2-D velocity distribution, derived from a wide-angle reflection and refraction profile, we calculated the travel times trough our model to the depth of 70 km along the POLAR profile (Janik et al., 2009).This profile is located near the centre of our study area.From the comparison of travel times through our crustal correction model and the 2-D model by Janik et al. (2009), we found the effect of local crustal velocity-depth variations to be minor compared to the effect of the variation in Moho depth (see Fig. 6b).
Figure 7 shows an example of crustal corrected travel time residuals for the same stations that were used in Fig. 4. The crustal corrected residuals for station LP62 generally exhibit slightly late arrivals for almost all azimuths, suggesting generally lower velocities beneath the station.This contrasts with the azimuthally pronounced variation of residuals for stations KIF and MSF that again corresponds with the station location above the boundaries between high-velocity and low-velocity regions in several layers.

Resolution and the synthetic tests
To evaluate the reliability of tomography inversion results, it is necessary to identify which parts of the resulting model are resolved well and which parts may contain artefacts caused by the method.In our study we analysed the resolution by scaling the diagonal elements of the resolution matrix (RDE) using layer-wise cross firing of rays and a series of synthetic tests, to find the optimal RDE limit to correspond to fairly well-resolved part of our model.The resolution assessment scheme is described in detail by Kissling (1988) for local earthquake tomography.Instead of the RDE, Kissling (1988) used the whole resolution matrix but this is not necessary for teleseismic tomography since the bias due to predominantly near-vertical rays and lack of horizontal rays is well known.
The resolution matrix can be used to evaluate the capabilities of the ray geometry and model parameter grid to resolve the velocity perturbations by inversion of the travel time data set.The diagonal elements of the resolution matrix correspond to a posteriori variances of the model parameters (Menke, 1984).Therefore, the closer the diagonal element is to a value of 1.0 and the smaller all off-diagonal elements are, the better the correspondent model parameter is resolved.Based on this criterion, the resolution is fair to reasonably good below the station network at all depths.Figure 8 shows the diagonal elements of the resolution matrix mapped at two horizontal sections through our study volume.The best resolution is below the centre of the network in the depth range of less than 360 km (see Fig. 8).The area with the best res- olution moves eastward at larger depths, corresponding with the region of most cross firing from events in the NE, E, and SE (see Fig. 8b).This is understandable as the stations in the western part of the network belong to the permanent Swedish National Seismic Network (SNSN) in Sweden, which was being updated during the POLENET/LAPNET data acquisition period.Additionally, the most common recorded back azimuth direction pointed towards east.As a consequence, the number of rays crossing the cells in the western part of the study area, especially from south to north, was smaller and the ray coverage sparser.The sensitivity of our data set to velocity heterogeneities in the upper mantle was tested with the "chequerboard" test.We constructed a model consisting of alternating positive and negative anomalies placed regularly through the study area in both vertical and horizontal directions.The magnitudes of the P wave velocity anomalies were ±2 % compared to the iasp91 reference model and the anomaly size was 160 km × 160 km × 120 km or 2 × 2 × 2 cells in the inversion grid.Between the anomalies, we left layers with no velocity perturbations as suggested by Sandoval et al. (2004).
Examples of the model at two depths (120 and 360 km) are shown in Fig. 9.
A synthetic data set was calculated using this model and the ray parameters of the real data set, and the resulting synthetic data set was inverted back to a velocity model.In the horizontal direction, the anomalies were generally well recovered below the station network at all depths.Also, in the vertical direction the recovery was good in the central and eastern parts of the study area.In the western part, however, there was some smearing especially in the north-south direction.Figure 9 compares the model used for computing the synthetic data set and the results after inversion with our selected damping value at selected depths of 120 and 360 km as well as through two vertical sections.
Additionally, the resolution was evaluated using synthetic tests with different structures simulating large-scale anomalies in the upper mantle.We speculated that there could be such structure in the upper mantle roughly below the Belomorian Mobile Belt (see Fig. 1).Hence, all three tests have an anomalous body there and an additional body crossing it diagonally to help us visualize how anomalies can affect each other.The synthetic test models and results are shown in Fig. 8. From the results, we can see that while the anomalies are recovered quite well, there is some leakage both upand downward.Similarly to chequerboard tests, we see that the resolution in the western part of the study area is not as good as in the eastern part.In the areas with fairly good resolution, leakage may generally extend into the next layer up or down from the anomalous body (i.e. up to 60 km), while in areas with poor resolution, the leakage can extend to two layers (up to 120 km).From the test (Fig. 8) we can also see that we obtain slight positive anomalies around the negative anomalous body and vice versa, marking an "overswinging" effect (Kissling et al., 2001).
The fairly well-resolved regions for each layer of our study area are denoted in Fig. 8.The regions are mainly based on an RDE value of 0.3, which was selected as the optimal RDE limit of our model based on ray coverage and synthetic test results.The regions may have slight deviations from the RDE based on contracting information from ray coverage and synthetic tests.Similar analysis for the rest of the depth layers can be found in the Supplement.
Synthetic tests were used to derive the damping parameter and singular value cut-off appropriate for data error and model parameterization.Results of those tests yielded a damping parameter of 100 and a minSV of 70 as the best choice for our data set.

Results
The main results of inversion with the real data are shown in Figs. 10 and 11.Our study revealed a highly heterogeneous lithospheric mantle beneath the northern Fennoscandian Shield, without any large high P wave velocity area that might indicate the presence of thick depleted lithospheric keel revealed beneath the southern part of the shield in Sweden (Shomali et al., 2006) and beneath the SVEKALAPKO study area (Sandoval et al., 2004), as well as in some other shield areas (e.g.Pasyanos and Nyblade, 2007;Priestley et al., 2008;Villemaire et al., 2012).The anomalies in the well-resolved part of our study area are concentrated in the upper part of the model volume, especially in the layer at 120 km depth (see Fig. 10).
While teleseismic tomography is not effective in resolving vertical variations of seismic velocities, the recent result of joint analysis of P and S wave receiver functions by Vinnik et al. (2016) revealed velocities close to those in the iasp91 model at upper mantle depths.Additionally, their results show that absolute average values of seismic velocities in the lithospheric mantle beneath southern Finland are generally higher than those beneath the northern Finland.
In our model, we can recognize several anomalies with higher velocities in the upper part of the lithospheric mantle (down to depths of about 120-200 km) that spatially correlate with the Karelian, Kola, and Norrbotten cratons (anomalies I, II, and III, respectively, in Fig. 10b the study area, the high velocities of anomalies II and III are stable results of inversions obtained with different inversion parameters.We interpret these higher velocity anomalies as non-reworked fragments of cratonic lithosphere preserved since the Archaean aeon.The non-reworked part of the Karelian craton (anomaly I in Fig. 11b) can be recognized as a higher velocity anomaly (+3 %) from the Moho down to a depth of about 160 km.Below this depth we observed velocities that are low compared to the overall velocity at the depth.At the margin of the Kola craton, we also see a positive (+3.5 %) anomaly roughly down to 200 km.Below the anomaly, velocity values are close to the average (Fig. 10d).The Norrbotten craton is also seen as a high-velocity anomaly (+1.5 %), starting from the Moho down to a depth of about 160 km (Fig. 10).
The most significant feature seen in the velocity model is a large negative velocity anomaly (up to −3.5 %, anomaly IV in Fig. 10b) in the central part of our study area that can be followed down to a depth of 160-200 km.In the upper part of the model this low-velocity area separates the three highvelocity regions corresponding to the cratons and it extends to the greater depth below the Karelian craton (Fig. 10b, c,  d).
The previous teleseismic tomography results close to the study area are based on SVEKALAPKO network data in southern and central Finland (Sandoval et al., 2004) and Swedish National Seismic Network (SNSN) data in Sweden (Shomali et al., 2006;Eken et al., 2012) south and southwest of our study area, respectively.While the comparison of the three models is not entirely straightforward, taking into account the relative nature of the teleseismic tomography method, particularly the anomalies obtained by Sandoval et al. (2004) fit our results very well in the overlapping part of the study areas.There is some similarity between anomalies obtained in our study and those revealed by SNSN data (Eken et al., 2012), but their comparison is less reliable due to the smaller overlapping area of both studies and lower resolution in that part of our study area.

Discussion
As is shown in several previous studies (e.g.Pedersen et al., 2013;Plomerová and Babuška, 2010;Lebedev et al., 2009), the lithospheric thickness in the POLENET/LAPNET study area in northern Fennoscandia must be at least roughly 150 km, which makes it unlikely that the upper boundary of the low-velocity anomaly IV (Fig. 10b) seen in our model at depths shallower than 100 km is the lithosphereasthenosphere boundary.In Fig. 12 we show a north-southdirected cross section through a model, based on surface waves by Pedersen et al. (2013), together with the location of the low-velocity zone found in our study.
The relative sizes of the anomalies recovered in the upper mantle of the POLENET/LAPNET study area are relatively large for a shield area; at 120 km depth, the negative anomaly IV is −3.5 %, while positive anomalies I and II are +3.5 %, resulting in a total of 7 % maximum velocity perturbation over the layer.Generally, the velocity variation in shield areas fits within 4 % variation (e.g.Kennett, 1997;Sandoval et al., 2004;Frederiksen et al., 2007;Youssof et al., 2015).On the other hand, large variations in P wave velocity in cratonic areas have been found recently close to our study area in the upper mantle beneath Sweden by Eken et al. (2012).
These upper mantle heterogeneities can generally be explained by three major factors, namely temperature variations, compositional variations, or seismic anisotropy.In addition, a water content in the mantle can enhance seismic attenuation and decrease seismic velocities (Karato, 2003).Seismic anisotropy beneath the POLENET/LAPNET area  2011) is collocated with the western part of the upper mantle high-velocity anomaly I of our model, with both northern and western boundaries fitting well together.North of the Archaean mantle block they found an area with fast S wave direction roughly towards west.This area is roughly collocated with our low-velocity anomaly IV though the northern boundaries of the domain by Plomerová et al. (2011) and our anomaly IV does not exhibit as good a fit as the southern ones do.Vinnik et al. (2014) analysed the southern half of the POLENET/LAPNET study area using joint inversion of P wave receiver functions and SKS recordings.In contrast with the work by Plomerová et al. (2011) they obtained a laterally averaged depth distribution of the magnitude and azimuthal direction of the anisotropy in the area but not lateral variations.They found an anisotropic layer with approximate S wave anisotropy of 2.5 % from below Moho to approximately 110 km depth with fast direction of 40-60 • , as well as another, slightly less anisotropic layer with the same azimuthal direction starting from approximately 220 km depth.Above this lower layer, they modelled another anisotropic layer with a contrasting azimuthal direction of 110 • .
While we do have reasonably good azimuthal coverage in our data set with seismic energy arriving from all directions, the events from a roughly east-west direction have on average slightly larger distances, and hence smaller incidence angles, than events from a south-north direction.This could cause the anisotropy to affect the results of the tomography that could explain some portion of the lateral velocity variations in our model.Eken et al. (2012) shows that these effects are likely to concentrate on less well-resolved areas at depths larger than 200 km but in our case they might have some effect in the uppermost mantle depths, too, as Plomerová et al. (2011) estimated the anisotropy beneath our uppermost mantle anomaly IV to have a contrasting fast direction of 55 • from north when compared with the roughly 0 • of the units both north and south of it.In summary, we may conclude that at least in some regions our -for the mantle lithosphere -significantly strong lateral velocity anomalies could in part be attributed to seismic anisotropy.
As shown by Lee (2003), James et al. (2004), andSchutt andLesher (2006), compositional variations between man-low P wave velocity anomaly mapped in this study  tle peridotites that are depleted in Fe and more fertile can explain up to 1-2 % velocity anomalies of P wave velocities.In order to explain larger anomalies, one would need the combined effect of major element chemistry and temperature as shown by Hieronymus and Goes (2010).That is why the lowered seismic velocities in the lithospheric mantle of the central part of our study area (anomaly IV compared to anomalies I, II and III) are probably most due to the combined effect of anisotropy, a more fertile composition, and a generally higher temperature.
The low-velocity zone spatially overlaps in the east with the Kola alkaline province (see Fig. 13), in which alkaline magmas have intruded the crust (Downes et al., 2005) during several metasomatic events.Kempton et al. (1995) proposed that at least one of them was ancient while the latest occurred in the Devonian.The same events would also result in extensive reworking and refertilization of the originally depleted Archaean mantle keel while the latest Devonian magmatism would also explain higher mantle temperatures.
On the other hand, the area of low seismic velocities in the central part of the network separates Norrbotten, Kola and Karelian cratons from each other and is spatially correlating with the northern part of the 1.9-1.8Ga N-S trending Baltic-Bothnia Megashear (BBMS) stretching below the Baltic Sea and the Gulf of Bothnia and continuing to the Caledonides in Norway (Berthelsen and Marker, 1986, see Fig. 13).The northern part of the BBMS was later named Pajala shear zone by Kärki et al. (1993).This shear zone is about 40 km wide, represented by a complex set of N-S striking shear and thrust zones.Lahtinen et al. (2015) proposed that the Pajala shear zone originated as a divergent plate boundary due to the collision of two Archaean continental units (Norrbotten and Karelian), and it was multiply reactivated after the continental collision with both lateral and vertical movements before 1.83 Ga.As the horizontal sections in the upper mantle obtained in our modelling suggest a depth distribution of the velocity perturbations beneath the megashear similar to that revealed by Bruneton et al. (2004) and Sandoval et al. (2004) (see Fig. 13), we suggest that the Pajala shear zone may continue to the south beneath the Gulf of Bothnia as was originally proposed by Berthelsen and Marker (1986).
However, the metasomatic processes and refertilization of the upper mantle in the Palaeoproterozoic alone would not produce such low seismic velocities as we observed in our model and combined effect of temperature and composition would be necessary (Hieronymus and Goes, 2010).Thus,

Figure 1 .
Figure 1.The map of the study area.(a) Simplified map of orogenies based on Lahtinen et al. (2015) in northern Fennoscandia.The Fennoscandian crustal block is outlined by a black dotted line and our study area by a rectangle.(b)The simplified geological map is based on the 1 : 2 000 000 geological map of Fennoscandia(Koistinen et al., 2001).The seismic stations are shown with black stars, triangles, and dots, and main geological provinces are named.

Figure 3 .
Figure 3. Example travel time picks of different quality from the same event observed at three different stations.The recorded waveforms have been filtered with the WWSSN short period simulation filter.The travel time observation from the station OUL was attributed a quality class 1, the one from station LP62 a quality class 2, and from station SGF a quality class 3.For OUL, the absolute travel time (P abs ) was picked as well as the relative travel time (P ).For each pick, the timing uncertainty estimate of the quality class is shown with a grey rectangle.

Figure 4 .
Figure 4. Distribution of travel time residuals relative to standard Earth model iasp91(Kennett and Engdahl, 1991) at three example stations selected from different parts of the POLENET/LAPNET network.The observed absolute residuals are on the left-hand side and the relative residuals after the average over all stations for each event has been removed are on the right-hand side.

Figure 5 .
Figure5.A comparison between data and model variance with different damping values.The selected damping value was first found using one inversion round only and the results are shown with diamonds with the damping value used shown on the right side of the symbols.The number of iterations was also optimized (dots with the number of iterations marked to their left side) for the selected damping value of 100.The final number of iterations, 4, is marked with a red circle.The grey line denotes the overall data uncertainty, 0.017 s 2 , calculated as the average uncertainty of all observations.

Figure 6 .
Figure 6.Crustal correction map.Panel (a) shows crustal correction values for all of the stations.The values are based on a travel time of a vertical seismic ray travelling through the crustal correction model described in Sect.4.2.Panel (b) shows the comparison of the travel times through a vertical section of the pseudo-3-D crustal correction model used in this study (black line) and the southern part of the POLAR wide-angle reflection and refraction profile (red line, Janik et al., 2009).The location of the comparison is shown by a black line in (a).
Distribution of crustal corrected travel time residuals at a few example stations from different parts of the POLENET/LAPNET array.

Figure 8 .
Figure 8. Two examples of data used to evaluate the lateral resolution of the inversion results with our data set at depths of (a) 120 km and (b) 360 km.The evaluation was based on diagonal elements of the resolution matrix (RDE) in the upper left corner of both subplots but also on ray coverage, shown in the upper right corner of both subplots, and the three synthetic test ST1, ST2, and ST3.Red dots and triangles in RDE plots show the locations of the POLENET/LAPNET and SVEKALAPKO stations, respectively, and the yellow line shows the boundary of the fairly well-resolved area.The boundary mostly follows a 0.3 RDE line but has occasionally been adjusted based on ray coverage and synthetic test results.In the right side of the figure, west-east vertical sections through the three synthetic models are shown, clarifying the location of the anomalous bodies in the vertical direction.The resolution analysis figures at all depths can be found in the Supplement.

Figure 9 .
Figure 9. Chequerboard test results.Two horizontal sections are shown at the depths of 120 (a) and 360 km (b) as well as two vertical sections, one in SN direction (c) and one in EW direction (d).For each section, both the model used to compute the synthetic data set and the result after inversion are shown.The model plot of subplot (a) shows the grid used and the locations of the vertical and horizontal sections with thicker lines, respectively.The locations of the anomalies are bordered with rectangles.The areas that were not inverted are marked with grey blocks.

Figure 10 .
Figure 10.Horizontal sections through our final inversion result.The depth of each section is marked above each plot.The grey corners show areas that were crossed by no rays and consecutively were excluded from inversion.The Roman numbers in (b) are referred to in the text when discussing the anomalies.

Figure 11 .
Figure 11.Vertical sections through our final inversion result.The locations of the sections are shown as black lines on the map also showing the horizontal grid.The depth grid is shown in plot of A-A .The areas where no inversion was done are shown in grey.

Figure 12 .
Figure 12.Comparison to previous lithosphere thickness results.The figure has been modified from Pedersen et al. (2013) and it shows the S wave velocity structure obtained in that study in comparison to the global S wave velocity model by Debayle and Ricard (2012).The area outlined with the yellow line shows the wellresolved part of our model at the depth of the low-velocity anomaly found in the upper mantle.

Figure 13 .
Figure13.Our results at the depth of 120 km together with cratonic units in our study area and locations of the Kola alkaline province, an upper mantle low-velocity anomaly found in SVEKALAPKO data byBruneton et al. (2004) andSandoval et al. (2004), and Baltic-Bothnia megashear.

Table 1 .
The travel time data quality classes.The table shows the data quality classes, the corresponding error estimates, and the number of travel time residual in the quality class for POLENET/LAPNET data, SVEKALAPKO data, and the total.The bottom row shows the total number of travel time residuals over all quality classes for both projects, and finally the total number of travel time residuals in our database.