One- and Three-dimensional P- and S-wave Velocity Models of Central and Southern Sweden Based on SNSN Data

Examensarbete vid Institutionen för geovetenskaper ISSN 1650-6553 Nr 294

One- and Three-dimensional P- and S-wave Velocity Models of Central and Southern Sweden Based on SNSN Data

One- and Three-dimensional P- and S-wave Velocity Models of Central and Southern Sweden Based on SNSN Data

Ne Xun Chan

Ne Xun Chan

Uppsala universitet, Institutionen för geovetenskaper Examensarbete E1, 45 hp i Geofysik

ISSN 1650-6553 Nr 294

The velocity structures of southern and central part of Sweden have

been derived with the local tomography (LET) method. The region has

been divided into two study areas and the datasets come from the P-

and S-wave traveltimes recorded by the Swedish National Seismic

Network (SNSN). Man-made explosions and earthquakes occurring

over the period of 5 years and 10 years, respectively, within the study

areas have been used. One-dimensional starting models were derived

based on an a priori model obtained from the SNSN, that were later

used for starting models in the inversion for the 3-D crustal structures

of the study areas. Attempts were also made to invert for Moho

topography in the areas. The study areas are found to have an upper-

crustal thickness of approximately 20 to 25 km and the Moho

boundaries vary from 42 to 46 km in depth. The V p /V s ratios varies

from about 1.68 to 1.78. The LET method appear to resolve the

different between the Sveconorgwegian and Svecofennian orogen

regions, but the stations and sources are too sparsely distributed for

higher resolution models. The seismicity in the study areas are

distributed in two distinctive depth ranges. The focal depth of the

SNSN catalogued earthquakes concentrated in approximately 5 km

and 15 - 20 km depth. Relocations of the earthquakes using a global

search method reduce this tendency. The results also show that using

3-D models produces less biased results than using 1-D models with

the same relocation method.

Supervisors: Ari Tryggvason Examensarbete vid Institutionen för geovetenskaper

ISSN 1650-6553 Nr 294

One- and Three-dimensional P- and S-wave Velocity Models of Central and Southern Sweden Based on SNSN Data

Ne Xun Chan

Copyright © /F
9VO
$IBO and the Department of Earth Sciences Uppsala University

Abstract

The velocity structures of southern and central part of Sweden have been derived with the local tomography (LET) method. The region has been divided into two study areas and the datasets come from the P- and S-wave traveltimes recorded by the Swedish Na- tional Seismic Network (SNSN). Man-made explosions and earthquakes occurring over the period of 5 years and 10 years, respectively, within the study areas have been used.

One-dimensional starting models were derived based on an a priori model obtained from

the SNSN, that were later used for starting models in the inversion for the 3-D crustal

structures of the study areas. Attempts were also made to invert for Moho topography in

the areas. The study areas are found to have an upper-crustal thickness of approximately

20 to 25 km and the Moho boundaries vary from 42 to 46 km in depth. The V p /V s ratios

varies from about 1.68 to 1.78. The LET method appear to resolve the different between

the Sveconorgwegian and Svecofennian orogen regions, but the stations and sources are

too sparsely distributed for higher resolution models. The seismicity in the study areas

are distributed in two distinctive depth ranges. The focal depth of the SNSN catalogued

earthquakes concentrated in approximately 5 km and 15 – 20 km depth. Relocations of

the earthquakes using a global search method reduce this tendency. The results also show

that using 3-D models produces less biased results than using 1-D models with the same

relocation method.

Contents

Abstract i

List of figures iii

1 Introduction 1

1.1 Study objectives . . . . 1

1.2 Local earthquake tomography (LET) . . . . 2

2 Theory 4 2.1 Forward calculations . . . . 4

2.2 The inverse problem . . . . 5

2.3 LET method . . . . 6

2.3.1 Joint inversion for P- and S-wave velocities . . . . 7

2.3.2 Hypocentre-structure coupling . . . . 9

2.3.3 Regularisation . . . . 11

3 SNSN Data 15 3.1 SNSN . . . . 15

3.2 Data selection . . . . 16

4 Modelling 19 4.1 One-dimensional modelling . . . . 19

4.2 Three-dimensional modelling . . . . 26

4.3 Checkerboard tests . . . . 28

4.4 Moho depth modelling . . . . 31

5 Results 33 5.1 Three-dimensional crustal models . . . . 33

5.2 Moho topography . . . . 34

5.3 Seismicity . . . . 34

6 Discussions 43 6.1 Upper-crustal velocity structures . . . . 43

6.2 Moho depth variations . . . . 45

6.3 Seismicity . . . . 45

7 Conclusions 47

Acknowledgements 48

References 49

List of Figures

1 Recorded local events by the SNSN in 2011. . . . 14

2 Locations of current SNSN stations . . . . 17

3 Map view of the central and southern study areas. . . . 18

4 Wadati plots with best-fitting lines for S traveltimes vs P traveltimes. . . . 21

5 Comparison of traveltime curves of the initial explosion data and theoretical traveltime curves calculated with the 1-D starting models. . . . 22

6 Plots comparing the traveltime residuals against offset after relocation of explo- sion events using the a priori 1-D model and the modified starting models. . . . 23

7 Comparison of theoretical epicentral traveltime curves of earthquake events and traveltime curves of relocated earthquake events . . . . 24

8 Comparison of SNSN a priori model and the modified starting models for the two study areas. . . . 25

9 Trends of weighting factors of regularisations and change in residual RMS from the first to the last iteration. . . . 27

10 Horizontal layer across the checkerboard models and the synethic inversion results. . . . 30

11 Map views of the resultant relative velocity structures and V _p /V _s ratio distribu- tions of 4 – 7 km depth layer of the two study areas. . . . 36

12 Cross-sections of the crustal relative velocity structures and V _p /V _s ratio distri- bution of the two study areas. . . . 38

13 Cross-sections of the absolute P-velocity structures of the two study areas and the results from the FENNOLORA study. . . . 39

14 Veritcal cross-sections of absolute P-velocity structures of the Moho regions in the two study areas. . . . 40

15 Three-dimensional Moho topographies showing the Moho variations of the two study areas. . . . 41

16 Histograms comparing the distributions of the source depth after relocations using different velocity models in the two study areas. . . . 42

17 A general geotectonic map view of central and southern Sweden. . . . 44

1 Introduction

Seismic tomography is an imaging method for modelling the velocity structures of the Earth’s interior using seismic signals, which can be generated by natural earthquakes or artificial explosions.

Geophysicists have been studying the seismic information from earthquakes since the early 20th century. The famous Mohorovi˘ci´c discontinuity (Moho) was first observed by Croatian seismologist Andrija Mohorovi˘ci´c in 1910, based on the seismic data collected from an earth- quake that occurred in October 1909, near Zegreb, Croatia (Mohoroviˇci´c, 1910). His discovery of the discontinuity of seismic wave velocity at the depth of about 50 km helped to establish the boundary between the crust and the mantle, which is now called Moho (Herak, 2012).

Seismology has since helped us to better understand the interior structure of the planet and explore the subsurface.

Pioneering seismic tomographic works in the 1970s focused on body-wave traveltime to- mography (e.g. Aki and Lee, 1976; Aki et al., 1977; Dziewonski et al., 1977), but developments had since followed with the use of different types of seismic data, such as waveforms, ampli- tudes, attenuation, or surface-waves. This study is using only traveltime data of body-waves to perform the tomography.

1.1 Study objectives

This thesis studies the crustal velocity structures and thicknesses in the selected southern and central part of Sweden with data obtained from a local seismic network, the Swedish National Seismic Network (SNSN).

Although Sweden is not lying on any major fault lines, earthquakes around the country

are observed in a daily basis, albeit small and mostly unnoticed, with moment magnitude

generally less than M _W = 2.5 (Böðvarsson et al., 2006). In addition, explosion sources gener-

ated from mines, quarries and construction sites provide most of the seismic data recorded by

the network (Fig. 1). Therefore, the SNSN has built a substantial traveltime database since its establishment in the late 1990s.

This study exploits the traveltime data from both explosion (on the surface) and earthquake sources (at depth) in the SNSN catalogue to perform a local earthquake tomography (LET) of two study areas respectively. An one-dimensional starting model is first derived for each study area using explosion events. Those 1-D models are also refined with the earthquake data. A three-dimensional crustal velocity model of each study area is then constructed with the newly-derived starting model. With the resultant 3-D models, I attempted to invert for the Moho depths within the areas and to investigate the earthquake depth distributions in south and central Sweden.

1.2 Local earthquake tomography (LET)

LET solves for the velocity distribution of a region within a receiver network, which usually covers a horizontal area of dimensions of several hundred kilometres or less.

The first arrival times of the body-waves are commonly employed as the data, for their

ease of extracting from the seismogram and their straightforward relationship with the wave

velocity. While the earthquakes processed in a teleseismic study are generally thousands of

kilometres away from the receivers, the earthquake sources involved in a LET normally come

from the targeted volume underneath the receiver array. The determination of the uncertain

source location is also included in the inversion. Aki and Lee (1976) demonstrated a simul-

taneous inversion of velocity model and hypocentres of earhtquakes, with first P arrival time

data from local earthquakes in California. This distinctive feature of determining hypocen-

tres and the velocity structure together is important. Thurber (1992) refers to this feature as

hypocentre-velocity structure coupling, and shows that neglecting the coupling results in po-

tentially significant bias of the models. This thesis is using the singular-value decomposition

method proposed by Pavlis and Booker (1980) to separate the two sets of parameters, source

and velocity model parameters, and allow the process of smaller matrices without decoupling the problem.

Previously, the P-wave (compressional wave) first arrival times were often the only choice for the LE tomographer given the poor quality of the S-wave (shear wave) arrival record- ings in some older networks with mostly vertical-component seismographs. However, the present three-component seismic networks can generally overcome this difficulty and record reasonably good S-wave data (Tryggvason et al., 2002). The inclusion of S-wave traveltime is beneficial as better hypocentre locations can be obtained (Gomberg et al., 1990), which give us an improved velocity model due to the hypocentral-velocity coupling. With comparable S- wave data quality and quantity as the P-wave data, a joint inversion of P- and S-wave velocity model is generally more desirable than an inversion using only P-wave data.

The LET computation in this study is done by a computer algorithm, PStomo_eq (Tryg-

gvason, 2012), which is based originally on fdtomo_eq (Benz et al., 1996). The finite difference

(Podvin and Lecomte, 1991) and the LSQR algorithm (Paige and Saunders, 1982) are included

in the algorithm. For full description of the algorithm, I refer to the user manual of the pro-

gramme (Tryggvason, 2012). The data selection and processing are handled by tools and filters

provided in the programme package. All the figures and plots showing the data and mod-

els in this thesis have been created by the Generic Mapping Tools (GMT) (Wessel and Smith,

2014).

2 Theory

The LET method involves solving forward and inverse problems. The joint inversion of P and S waves are described here, as well as the coupling of relocations and inversions for velocity structures.

2.1 Forward calculations

The forward problem in the seismic tomography calculates the traveltime of a seismic wave from sources to receivers by ray tracing methods. The traveltime of a ray can be determined given a velocity model of the medium, and the origin time and location of the source. Each ray is the normal to a wavefront. The Fermat’s principle states that a ray follows the path that takes the minimum time to travel (Lay and Wallace, 1995). The traveltimes are, therefore, the shortest time that a ray travels from the source to the receiver. The LET in this thesis concerns only with local scale, and therefore a flat Earth is assumed.

Different ray tracing techniques have been developed over the years. For example, ray bending (Thurber, 1983) and shooting (Benz and Smith, 1984) techniques can be used to find a ray path between a source and a receiver by refining an initial path using iterative methods.

They are efficient in finding the first arrival path in gradually varying structures, but have difficulties in highly contrasted velocity structures (Benz et al., 1996).

A finite-difference method (Vidale, 1988) is used here. The study region is parametrised into cells, and the wave propagation is described using a finite difference approximation of the eikonal equation

(5 t ( x )) ² = s ( x ) ² (1)

where x is a positional vector, s ( x ) is the slowness (reciprocal of the wave velocity) of the

medium and the t ( x ) is the arrival time of the wavefront at point x. The eikonal equation

relates the wavefronts with the seismic velocity distribution. The traveltime can be obtained

by solving the partial differential equation. One can addressed the problem with Huygens’

principle (Podvin and Lecomte, 1991), which states that every point on a wavefront can be seen as the secondary source with a wavelet that travels outwards (Lay and Wallace, 1995).

The tangent of the expending waves can be connected to give the next wavefront. If the targeted volume is parametrised into cells of the same size, each point on the boundary can be trended as the secondary source that emits an impulse according to Huygens’ principle.

The first arrival from the set of secondary sources is then picked. Therefore, the traveltimes are computed at some new points using values from the adjacent points and the gradients in Eq. (1) are approximated (see Podvin and Lecomte, 1991, for detailed formulations).

Solving the forward problem requires a given velocity model, and initial locations and times of the sources. Therefore, an initially guessed model is provided and is improved using an iterative process in the inverse problem. The forward calculation provides data predictions for each iterative step.

2.2 The inverse problem

LET is a highly non-linear inverse problem because the path taken by a seismic wave from the source to the receiver is also a function of the velocity model. It is typical to linearise the problem by calculating the misfits between the observed data and the predicated data. The relationship between the dataset and the model can be illustrated by

d = g ( m ) (2)

where d is a vector containing the seismic traveltime data, and g ( m ) is a function of the velocity model, m. The predicated data can be calculated using an a priori velocity model. The difference between the observed data and the predicated data is

∆d = _d ^obs − _g ( _{m 0} ) ₍₃₎

where m ₀ is the reference model, and g ( m ₀ ) is the predicted data, which are computed using the reference model by a forward problem. Using this data misfit, one can compute a model perturbation, ∆m, by

∆d = g ( _∆m ) (4)

The tomography inverse problem is typical neither purely over-determined nor under- determined, as the number of data is usually more than the number of parameters, yet some of the model cells may not be covered by any ray. That results in a mixed-determined problem.

The objective function, Φ ( m ) , with a regularisation term is in the form of (Menke, 1984)

Φ = ( _g ( _m ) − _d ^obs ) ^T _C ⁻ _d ¹ ( _g ( _m ) − _d ^obs )+

ε ( _m − _{m 0} ) ^T _C ⁻ _m ¹ ( _m − _{m 0} ) + η ( _m ^T _D ^T _Dm )

(5)

where m ₀ is the starting model, D is a smoothing operator, C ⁻ _d ¹ and C ⁻ _m ¹ are the covariance for the a priori data and model respectively, and g ( _m ) and d ^obs are the predicted and the observed data, respectively. The two parameters, η and ε, are controlling the degree of damping and smoothing of the model respectively. These damping and smoothing factors determine the trade-off between the model resolution and smoothness. The model perturbation, ∆m is updated for every iteration, and it is followed by a data update from the forward calculation.

2.3 LET method

The inversion formulations presented here follow mostly from Tryggvason (1998). The scheme

has considered the joint inversion of P- and S-velocities and the coupling of hypocentres and

velocity structures.

2.3.1 Joint inversion for P- and S-wave velocities

The traveltime of a body-wave can be related to the velocity structure by a line integral of the slowness (reciprocal of velocity) of the medium between the source and the receiver. The non-linear traveltime equation for an individual ray path from earthquake source i with an origin time, τ, to receiver j, can be expressed as

T _ij = τ _ij +

Z

raypath u ( r ) dl (6)

where r is a three-dimensional position vector and dl is the differential length of the ray path.

The equation is non-linear because the raypath that is taken by the seismic wave energy also depends on the velocity model. The raypath as a function of the slowness can be represented by l [ u ( r )] . The linearisation of the traveltime equation is carried out using the difference between the observed data and the theoretical data, which is calculated based on the reference velocity model and the starting earthquake locations. This difference is called the traveltime residual, and for earthquake i to receiver j, the residual is

r _ij = τ _i + t ^obs _ij − t _ij ^pred (7)

where t ^obs _ij and t ^pred _ij are the observed and theoretical traveltime respectively. The model is parametrised into uniformly-sized blocks with constant slowness. By the first order Taylor series approximation, a linearised relationship is established between the traveltime residual and the velocity perturbation about the reference model:

r _ij =

∑ 3 k = 1

∂T _ij

∂x _ik ∆x ik + _{∆τ i} +

∑ L l = 1

∂T _ij

∂u _l ∆u l (8)

where T _ij is the traveltime from earthquake i to receiver j, and ∆x ik , ∆u l and ∆τ i are the perturbations to the earthquake location, slowness and origin time respectively. ^∂T _∂x ^ij

ik and ^∂T _∂u ^ij

l

represent the partial derivatives of traveltime with respect to the coordinates of the earth- quake location and to the slowness model respectively. L is the number of model blocks that particular ray has propagated through.

Under this parametrisation scheme, a ray should travel with constant velocity and in a straight path within an individual model block. The traveltime of earthquake i to receiver j within the n-th cell is expressed as

T _ijn = u n l n (9)

where u _n is the slowness of the n-th cell and l _n is the length of the ray inside the cell. If one assumes the raypaths in the cell are not affected by the other paths of the rays outside the block, the partial derivative of traveltime with respect to the slowness model in Eq. (8) can be approximated for every block, as

∂T _ij

∂u n

≈ ^{∂T ijn}

∂u n

= ^∂ ( u _n l _n )

∂u n

= l _n (10)

A similar assumption of Eq. (10) can apply to the cell k with source e inside, the length of the ray i from the source to the exit-point of the cell will be

l _e = v u u t

∑ 3 k = l

( ε _k − x _ik ) ² (11)

and the partial derivative of traveltime with respect to the location can be approximated by

∂T _ik

∂x _ik ≈ ^{∂T ije}

∂x _ik = ^{∂u e l e}

∂x _ik = − ^{x ik u e}

l _e (12)

For all the observations of an individual earthquake i, Eq. (8) can be grouped into a system of equations in the matrix form

γ _i = A _i ∆h i + B _i ∆u (13)

where γ _i is the vector containing all the traveltime residuals for earthquake i, A _i is the matrix containing the traveltime partial derivatives with respect to the hypocentre locations, ∆h i is the vector of hypocentre perturbations (∆x ik and ∆τ i ), B _i is the matrix of partial derivatives of

traveltime with respect to the slowness perturbations, and ∆u is the vector with the slowness perturbations.

If the total number of sources is P, Eq. (13) can be arranged into one matrix equation as

























A ₁ 0 · · · 0 B ₁ 0 A 2 .. . B 2

.. . . .. 0 .. .

0 · · · 0 A _P B _P

























































∆h 1

∆h 2

.. .

∆h P

∆u

































=























 γ ₁

γ ₂ .. .

γ _P

























(14)

To simultaneously inverse the S-wave velocity, Eq. (13) should also contain the partial derivative matrix of the the S-wave traveltimes, and expressed as

γ ^P,S _i = A ^P,S _i ∆h i + B ^P _i ∆u ^P + B ^S _i ∆u ^S (15)

where the superscripts P and S indicate P- and S-wave model respectively. A clearer indication of the coupling can be shown in a diagrammatic form of Eq. (15)







 γ _i ^P γ ^S _i









=









A _i ^P B _i ^P 0 A ^S _i 0 B ^S _i

























∆h i

∆u ^P

∆u ^S

















(16)

2.3.2 Hypocentre-structure coupling

The difficulty with hypocentre-structure coupling is that the size of the kernel matrix increases

rapidly with observations and earthquakes. From Eq. (16), one can see that the dimension of

the equation depends on the number of earthquakes and receivers. If the total number of observations is M and the number of model parameters is N, the dimension of the combined matrix is M × ( 4P + N ) . Any addition of data quickly increase the size of the problem as every earthquake takes up four columns and each observation takes one row. A large number of sources is generally required in an LET for the raypaths to sufficiently cover the study region.

As a result, the problem is almost intractable.

As mentioned previously, decoupling the problem is not desirable. One method is to sep- arate out the hypocentre parameters from the velocity parameters. The problem is, however, that it is not apparent what proportion of the data misfit goes to earthquake location and what goes to the velocity model. Although there are usually more than enough data to locate earth- quakes, as each of which needs only four parameters to describe, the chances are that there are parts of the model that are not covered, or poorly constrained and under-determined. The inverse problem is, therefore, mixed-determined. Pavlis and Booker (1980) shows a method that uses singular-value decomposition procedure to separate the earthquake data that cov- ers the hypocentres and the remaining data for the velocity inversion. This allows locating the hypocentres and then estimating a linearised perturbation to the velocity model without the need of decoupling. They make use of the orthogonal properties of the hypocentre partial derivative matrix A _i in Eq. 13 to reduce the number of unknowns. In doing so, the dimensions of the matrices can be reduced and the efficiency of the inversion process can be improved.

Matrix A _i can be decomposed by singular-value decomposition in the principal axis system by

A _i = _UΛV ^T (17)

where U is an M × M square matrix containing the orthogonal eigenvectors of AA ^T , Λ is

an M × N diagonal matrix with the singular values of A, and V is a N × N square matrix

containing the orthogonal eigenvectors of A ^T A (Menke, 1984). Then, U can be partitioned

and be further defined as

U = [ U p U _l U 0 ] (18)

where U p contains the eigenvectors of the first four discrete parameters concerning the hypo- central location, U _l and U ₀ contain the eigenvectors of the remaining parameters without information about the hypocentral location. U _l contains the parameters that are poorly con- strained and cannot be determined for the hypocentral location and U 0 contains the eigenvec- tors associated with the zero singular values.

Using the orthogonal properties of U,

U ^T _{0 i} A i = _{0 (19)}

and operate it on to Eq. (13) gives

U ₀ ^T _i γ = _U ^T ₀

i A _i ∆h i + _U ^T ₀

i B _i ∆u (20)

The dependence on ∆h i is thus eliminated:

γ ⁰ _i = _B ⁰ _{i ∆u (21)}

which shows the transformed data and traveltime derivative with respect to the slowness perturbation. The result leads to the reduction of the dimensions of the matrix in Eq. (14), from M × ( 4P + N ) to ( M − 4P ) × N. The prime indicates the transformed kernel.

2.3.3 Regularisation

Regularisations are used to control between the roughness and resolution of the solution as shown in Eq. (5). There are three major constraints in this LET inversion.

The smoothing constraint is formulated by the Laplacian of the slowness perturbation field

of the adjacent cells of cell at i, j, k.

6∆u i,j,k − ( _{∆u i − 1,j,k} + _{∆u i + 1,j,k} + _{∆u i,j − 1,k} + _{∆u i,j + 1,k} + _{∆u i,j,k − 1} + _{∆u i,j,k + 1} ) = _{0 (22)}

where ∆u i,j,k is the slowness perturbation of the cell at i, j, k and the rest are the slowness perturbations of the adjacent cells. This smooths out the contrast between adjacent cells.

The ratio of P- and S-wave velocity (V _p /V _s ) can be used as a regularisation parameter.

The V _p /V _s ratio is damped to prevent strong variations, and is controlled with the damping equation of the form of

σ∆u ^P _i − _∆u ^S _i = 0 (23)

where u ^P _i and u ^S _i are the P- and S-wave slowness of block i respectively, and σ is an a priori V _p /V _s ratio. The weighting parameter to the equation determines the degree of deviation of the ratio is allowed.

The V _p /V _s ratio, however, can bias the model due to the influence from the a priori ratio (Tryggvason and Linde, 2006). Tryggvason and Linde (2006) propose a structural constraint in the form of a cross-gradients function, which is a cross product of the gradients of P- and S-wave slowness perturbations. The cross-gradients function t is defined as (Gallardo and Meju, 2003; Tryggvason and Linde, 2006):

t ( x, y, z ) = ∇ m ^P ( x, y, z ) × ∇ m ^S ( x, y, z ) (24)

where ∇ m ^P and ∇ m ^S are the P- and S-wave slowness model perturbations repectively.

The cross-gradients are zero when the two gradient vectors are parallel, i.e. spatial changes

occurring when both the P- and S-wave slowness are pointing in the same or opposite direc-

tion with each other. On the contrary, the cross products are largest when the two gradients

are perpendicular to each other. It is assumed a boundary is sensed by both velocity struc-

tures in a common orientation regardless of the amplitude of the physical property changes.

By minimising the cross-gradients, one forces the P- and S-slowness to change in similar directions. Additionally, t is still minimised, even when ∇ m ^P and ∇ m ^S are zero.

After adding the regularisation terms, the system of equations (Eq. 13) is then given by























 γ ⁰

0 0 0

























=























 B ⁰ kL lS µC

























∆u (25)

where B ⁰ is matrix with the transformed partial derivative with respect to the slowness. L is the Laplacian operator controlling the roughness of the slowness, and S controls the V _p /V _s damping, and k and l are weighing parameters for the smoothness constraint and the V _p /V _s damping respectively. C is the cross-gradients function and µ is a constant weighting factor for the cross-gradient.

A conjugate gradient solver LSQR (Paige and Saunders, 1982) is then used to solve Eq. (25).

The LSQR solves for sparse linear equations and is stable for under-determined or mixed-

determined systems (Tryggvason, 1998).

Figure 1: Recorded local events by the SNSN in 2011. The picture is taken from the SNSN

report 2011 (Böðvarsson, 2012).

3 SNSN Data

All traveltime data used in this study were recorded by the SNSN. They are first arrivals of P- and S-waves. There are two study areas, which are determined according to the catalogued event and station distributions. A set of arrival times is then selected for each of the study areas from events within that area.

3.1 SNSN

SNSN is a permanent seismic network with more than 60 seismic stations installed across Sweden (Fig. 2). The SNSN has been in operation since 1998 by the Uppsala University. All data are digitally recorded, and send to a central database. The data acquisition system is similar to the one developed for the SIL seismic network in Iceland (Böðvarsson et al., 1999;

Böðvarsson, 1999). Local earthquakes are automatically located as any triggering phase from any station is reported to the data centre, and data from all the stations are combined to create a catalogue of located events. The output event lists are improved and updated by manual inspection of waveforms, in order to optimise the location accuracy.

Most of the local events located by the network are explosive-induced earthquakes, and

only a small number are located as true earthquakes, along with undetermined events outside

the network. According to Böðvarsson (2012), there were nearly 6000 explosions out of 9000

detected events and only approximately 550 events were identified as true earthquakes in the

year 2011 by the network (Fig. 1). Most of the explosions are concentrated in Southern region

due to mining and construction activities and the earthquakes are mostly located along the

Baltic Sea coast, which can be accounted to the post-glacial rebound of the region (Böðvarsson

et al., 2006).

3.2 Data selection

Data are selected from both the explosion and earthquake events. The central region catalogue contains events in Sweden of latitudes between 59.8N and 62N, while the southern events are located from the latitude of about 55.5N to 59N. The earthquake data are gathered from the SNSN database between the year 2002 and 2012 without any further processing before the data selection procedure. However, a list of automatically picked explosion events from the catalogue, which were located within the southern and central region of Sweden between the year 2005 and 2010, were collected and repicked manually with a locating software, Lokimp, to improve the traveltime data. All events catalogued as explosions are assumed to come from quarries and constructions, and are referred to as untimed explosions (Uxp) in the follow- ing. Their initial time and hypocentres are uncertain, but it is reasonable to assume they are generated at the surface, or very close to the surface. I, therefore, take the coordinates of the epicentres of the Uxp and match them with the topography of the area from a gridded data file. The depths of Uxp are then moved to the elevation of the epicentres. The depths of the Uxp will remain the same throughout the inversion process.

By showing the events in map views, two study areas are then selected to cover most of the events and stations in that area (Fig. 3). The central region has an area of 560 × 400 km ² , within which there are 24 stations. There are 740 explosions and 1497 earthquakes located in the area. The local moment magnitudes of the earthquakes range from -1 to 3.7 and their depths range from the surface to about 40 km depth. The southern region has an area of 480 × _{480 km} ² , and contains 27 stations. The number of explosions and earthquakes located in the area are 1703 and 684 respectively. The magnitudes of the earthquakes range from -1 to 4.4. The depths of the sources range from the surface to 41 km.

Each event requires 7 or more recorded phases (P and S) to be selected for the tomography.

Because of the sparse distribution of events and stations, not all events are recorded with

sufficient phases. Eventually, there are 707 Uxp (4.5% reduction) and 785 earthquakes (47.6%

12.0˚E 16.0˚E 20.0˚E 24.0˚E 56.0˚N

60.0˚N 64.0˚N

68.0˚N

Figure 2: Locations of current SNSN stations (in red triangles)

reduction) selected for the central region, resulting in a total of 10576 and 5006 P and S phases,

respectively. The number of Uxp and earthquakes in the southern region were reduced to 1599

(6.1% reduction) and 477 (30.3% reduction), respectively. There are a total of 15498 P- and 9010

S-wave phases in this region. The explosion data have all been manually revised, which may

account for the much smaller reduction rates of the Uxp sources.

(a)

10.0˚E 10.0˚E

12.0˚E 12.0˚E

14.0˚E 14.0˚E

16.0˚E 16.0˚E

18.0˚E 18.0˚E

20.0˚E 20.0˚E

59.0˚N

59.0˚N 60.0˚N

60.0˚N 61.0˚N

61.0˚N 62.0˚N

62.0˚N

(b)

10.0˚E 10.0˚E

12.0˚E 12.0˚E

14.0˚E 14.0˚E

16.0˚E 16.0˚E

18.0˚E 18.0˚E 20.0˚E

55.0˚N

55.0˚N 56.0˚N

56.0˚N 57.0˚N

57.0˚N 58.0˚N

58.0˚N 59.0˚N

59.0˚N 60.0˚N

60.0˚N

Figure 3: Map showing the (a) central and (b) southern Sweden respectively. The study areas are indicated with the rectangular boxes, where each cross indicates the origin of that area. Red and blue circles are the earthquakes and explosions in the catalogue respectively.

Triangles are the SNSN stations.

4 Modelling

The modelling of the study areas is done in two steps. First, an 1-D starting model for each study area is derived by forward traveltime computations using Uxp sources, and is further refined with earthquake sources. Then, the 3-D models of the upper crusts are inverted iteratively with the 1-D models as the starting models. In addition, the Moho variations are modelled starting with the resultant 3-D models.

4.1 One-dimensional modelling

The SNSN uses an 1-D velocity model for the routine event location in the whole network.

The velocities in the 1-D starting models increase gradually with depth. The P-velocity model is used here as an a priori starting velocity model for both of the study areas. However, it is discretised into 3 km thick layers in my modelling instead of being a continuous model. The S- velocity models are derived from the P-wave velocity model by dividing a V _p /V _s ratio, which is determined by the Wadati plot (Fig. 4a and 4b). The slopes of the best-fitting lines show that the regional V p /V s ratios are 1.72 and 1.73 for the central and southern region respectively.

They are also used as the σ in the V p/Vs damping equation (Eq. 23). It is important that the starting models are appropriate in order to produce a well-resolved 3-D structure (Kissling et al., 1994). As the a priori starting model has a different representation from the original SNSN model, the model may not be appropriate. To see if the model can produce similar traveltime results as the initial SNSN data, it is used to calculate the traveltimes of P-arrivals of the Uxp events in the two study areas respectively. The theoretical traveltime curves for both study areas are plottedtogether with the respective observed traveltime curves (Fig. 5a and 5b). In both study areas, the theoretical traveltime curves do not match with the observed traveltime curves. After relocation of the Uxp events using the a priori models, the traveltime residuals are computed and plotted against the offset of the sources from the stations in Fig.

6 (black points). A reasonably good starting model is expected to have the mean residuals

close to zero, and no significant trend with offset (Tryggvason, 2012). The resultant plots from the SNSN model, however, show that the residuals deviate from zero as the offset increase for both P and S traveltimes in both study areas. The a priori SNSN model, therefore, does not seem to be an appropriate starting model in this representation. Thus, it is necessary to improve the 1-D model. The a priori model is then modified respectively for each study region in a trial and error procedure to match the relocated travel time curves with the initial travel time curves. After some trials, one model is derived for each study area. The comparing traveltime curves in respective study areas fit much better with each other than that with the a priori model (Fig. 5c and 5d). The traveltime residuals are calculated using the modified models, and are plotted against offset in Fig. 6 (red). The new residuals deviate less from zero then previously (black) for both study areas. The linear trends are also less obvious then before. For the central study area, the root-mean-square (RMS) values of all residuals are 0.675 and 1.168 for the P- and S-wave model respectively, using the a priori model. My model gives a residual RMS of 0.370 and 0.391 (45.0% and 66.5% improvement) for P and S respectively.

For the southern study area, the values are reduced from 0.993 to 0.759 (23.6% improvement) for the P-wave, and from 1.745 to 1.440 (17.5% improvement) for the S-wave.

The modified models are used to relocate the earthquake data, in order to test their per- formance at depth. The traveltime curves of the relocated data are compared with that of the theoretical traveltime curves with epicentral distance (Fig. 7). The coloured points represent the traveltime curves of the relocated data, and the coloured lines are the theoretical traveltime curves at the corresponding depth ranges. Each colour represents the sources within a 5 km depth range, from 0 to 30 km depth, respectively. Without further alternation, the relocated travel time curves show good agreements with the travel time curves of the relocated data (Fig. 7) in both study areas.

Fig. 8 shows the comparison between the SNSN models and the final modified models in

the two areas. The final models are similar to the a priori model, but generally have higher

velocity for each layer. The two modified models have shallower Moho discontinuity than the model. The Moho in the central region appears around 42 km, and the Moho discontinuity for the southern region appears to be at a shallower depth of 40 km. The respective model for the central and southern regions are otherwise similar. The starting S-velocity models are derived by dividing the P-models with the designated a priori V p /V s ratio of the study areas respectively.

(a)

0 5 10 15 20 25 30 35 40 45 50 55 60 65 70

S−wave traveltime (s)

0 5 10 15 20 25 30 35 40

P−wave traveltime (s) Central

(b)

0 5 10 15 20 25 30 35 40 45 50 55 60 65 70

S−wave traveltime (s)

0 5 10 15 20 25 30 35 40

P−wave traveltime (s) South

Figure 4: Wadati plots with best-fitting lines for S traveltimes vs P traveltimes. The slope for

(a) the central study area is 1.72, and that for (b) the southern study area is 1.73.

(a)

−4.0

−3.5

−3.0

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5

Traveltime (s) /6

0 50 100 150 200 250

Distance (km) Central (SNSN model)

(b)

−4.0

−3.5

−3.0

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5

Traveltime (s) /6

0 50 100 150 200 250

Distance (km) South (SNSN model)

(c)

−4.0

−3.5

−3.0

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5

Traveltime (s) /6

0 50 100 150 200 250

Distance (km) Central (modified model)

(d)

−4.0

−3.5

−3.0

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5

Traveltime (s) /6

0 50 100 150 200 250

Distance (km) South (modified model)

Figure 5: The reduced traveltime curves of the initial explosion data (black points) are com- pared with the reduced traveltime curves calculated with the 1-D starting models (red points).

(a) and (b) indicate the traveltimes curve plotted from the a priori SNSN model for the cen-

tral and southern study area respectively. (c) and (d) are showing that computed from the

modified 1-D models for the central and southern study area respectively.

(a)

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5 2.0 2.5

P−traveltime residual (s)

0 20 40 60 80 100 120 140 160 180 200 220

Distance (km) Central (P)

(b)

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5 2.0 2.5

S−traveltime residual (s)

0 20 40 60 80 100 120 140 160 180 200 220

Distance (km) Central (S)

(c)

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5 2.0 2.5

P−traveltime residual (s)

0 20 40 60 80 100 120 140 160 180 200 220

Distance (km) South (P)

(d)

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5 2.0 2.5

S−traveltime residual (s)

0 20 40 60 80 100 120 140 160 180 200 220

Distance (km) South (S)

Figure 6: The traveltime residuals are plotted against offsets. They are calculated from reloca-

tion of initial data using the parametrised a priori SNSN model (black points) and that using

the modified models (red points) respectively. (a) and (b) are P and S residuals in central

study area respectively. (c) and (d) are P and S residuals in the southern study area respec-

tively. The blue lines indicate the ideal median (zero) for the residual spreading, while the

green lines represents the residual thresholds in the 3-D modelling.

(a)

−4.0

−3.5

−3.0

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

Traveltime (s) /6

0 50 100 150 200 250

Distance (km) Central

Depth range 5−10 km 10−15km 15−20km 20−25km 25−30km

(b)

−4.0

−3.5

−3.0

−2.5

−2.0

−1.5

−1.0

−0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

Traveltime (s) /6

0 50 100 150 200 250

Distance (km) South

Depth range 5−10 km 10−15km 15−20km 20−25km 25−30km

Figure 7: The reduced epicentral traveltime curves are calculated for the relocated earth- quakes using the modified starting models, within particular depth ranges (coloured points).

Coloured lines represent the expected reduced traveltime curves calculated using the same

models, respectively, for the central and southern study area, and that of the relocated earth-

quakes in different depth ranges, for (a) the central and (b) southern study area respectively.

0 5 10 15 20 25 30 35 40 45 50 55 60 65 70

Depth (km)

3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5

P−wave Velocity (km/s)

Starting models SNSN

Central modified Southern modified

Figure 8: Comparison of the modified P-starting models for the two study areas and the SNSN

1-D model in discrete representation.

4.2 Three-dimensional modelling

A 3-D velocity model is computed separately for each study area. Each model contains cells of size 10 × 10 × 3 km ³ (horizontal × horizontal × vertical) from -2 to 7 km depth (negative sign being above sealevel). Under 7 km depth, the cell size is increased to 40 × 40 × 3 km ³ , in order to increase the coverage across the deeper crust and obtain a relatively smooth model. The upper crustal modellings use only arrivals from the crust (P g and S g phases). The maximum offset for any ray used in the inversions is limited to the crossover distance of 140 km. The modelling for each region was done over five iterations.

Model regularisations are controlled by respective weighting factors of the smoothing equation (Eq. 22), the V _p /V _s damper (Eq. 23), and the cross-gradients (Eq. 24). The smoothing weighting takes a large value at the first iterations, and subsequently reduced in later itera- tions (from 200 to 25 in five iterations). The V p /V s ratio is constrained by the V p /V s damper.

The dampers are set to 50 in the beginning and are reduced in each of the successive iterations (10 in the last iteration). The cross-gradient constraints are used in the second to the fourth iterations. Tests on different weighting parameters concluded these weight factors. Addi- tionally, the velocity models are forced to give increasing velocities with increasing depth to maintain smooth models without low-velocity zones.

Residual thresholds are introduced to control the deviations of the misfits. The P residual thresholds are set to 1.0 s for the first two iterations and tightened to 0.5 s for the rest of the iterations. The S residual thresholds are set to 1.8 s for the first two iterations and 0.8 s for the rest. There are more allowances for the S residual since S-waves travel slower (hence having larger traveltimes) than P-waves.

Eventually, the LSQR solver takes up to 99 steps to solve the inversion problem. The last

models are then used to calculate the final traveltime residuals. Fig. 9 show the number

of traveltimes and residual RMS after every iteration. For the central study area, the final

number of P- and S-wave traveltimes are both reduced by 12.0% (7122 to 6268 and 3767 to

3313 respectively). The number of traveltimes are reduced in successive iteration because phases excessing the residual limits are excluded. The reduction is most significant in the third iteration as residual thresholds are lowered from previous iterations. The traveltime RMS decreased from the first iteration of 0.238 to 0.113 s (52.5% reduction) for the P-waves and 0.609 to 0.131 s (78.5% reduction) for the S-waves. For the southern study area, the P- and S-traveltimes are reduced by 10.9% (9654 to 8602) and 8.1% (6420 to 5887) for the same reason as in the central study area. The P-wave residuals decreased from 0.228 to 0.138 s (39.5%

reduction) and the S-wave residuals decreased from 0.490 to 0.197 s (59.8% reduction).

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

RMS data fit (s)

1 2 3 4 5 6

Iteration

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

Max residual (s)

6000 6500 7000 7500 8000

Number of P−traveltimes

3000 3250 3500 3750 4000

Number of S−traveltimes

0 50 100 150 200

Smoother

0 50 100 150 200

Vp/Vs damper

Central

(a)

(b)

(c)

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

RMS data fit (s)

1 2 3 4 5 6

Iteration

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

Max residual (s)

8000 8500 9000 9500 10000

Number of P−traveltimes

5500 5750 6000 6250 6500

Number of S−traveltimes

0 50 100 150 200

Smoother

0 50 100 150 200

Vp/Vs damper

South

(d)

(e)

(f)

Figure 9: (a) and (d) are plots showing the trends of the weighting factors of the smoother

(solid lines) and the V _p /V _s damper (dotted lines) in the central and southern study area

respectively. (b) and (e) are plots showing number of respective P- and S-phases from the

first to last iteration in the central and southern study area respectively. (c) and (f) show the

changes in residual RMS (solid lines) and residual thresholds (dashed lines) for central and

southern study area respectively. Blue lines indicate properties for P-waves and black lines

are for S-waves.

4.3 Checkerboard tests

A checkerboard test is used to test the model resolution of the particular raypath geometry of a tomography. The test starts with creating a 3-D synthetic model by adding velocity perturbations to the starting model in a checkerboard cell pattern. A set of synthetic traveltime data are generated for the checkerboard model, with the same source-and-receiver geometry as the real dataset. Noise is added to the traveltime dataset to make the test more realistic.

Then, the synthetic data are inverted to resolve the checkerboard velocity structure with the same regularisation parameters and number of iterations as those in the actual inversion. The resultant models indicate the recovery of the initial synthetic model and is to some extent a test for which areas the models are reliable.

The 3-D synthetic checkerboards models are generated by superposing 5 × 5 × 5 number of model cells (50 × 50 × 15 km ³ checkerboard cells) ± 10% velocity perturbation on the start- ing 1-D models. The models are tested with different noise levels added to the traveltimes.

The noise levels vary from traveltime standard deviations of 50 to 500 ms. The respective inversion results are similar until 500 ms, which is a much larger uncertainty than the final RMS residuals in the 3-D inversions in both areas. The test results of both study areas with 200 ms noise level are shown in Fig. 10. The reconstructed models are smooth but in good agreement with the synthetic model. The resolution is similar for both P- and S-wave models.

The resolution of the models is best in the middle of the study area. and, particularly in the

central region, where stations density is highest.

0 50 100 150 200 250 300 350 400

Y (km)

0 50 100 150 200 250 300 350 400 450 500 550 X (km)

10.0˚E 12.0˚E 14.0˚E 16.0˚E 18.0˚E 20.0˚E

59.0˚N 60.0˚N 61.0˚N

(a)

Depth: 4−7km

0 50 100 150 200 250 300 350 400

Y (km)

0 50 100 150 200 250 300 350 400 450 500 550 X (km)

10.0˚E 12.0˚E 14.0˚E 16.0˚E 18.0˚E 20.0˚E

59.0˚N 60.0˚N 61.0˚N

1.64 1.66 1.68 1.70 1.72 1.74 1.76 1.78 1.80

Vp/Vs ratio

(b)

Depth: 4−7km

0 50 100 150 200 250 300 350 400

Y (km)

0 50 100 150 200 250 300 350 400 450 500 550 X (km)

10.0˚E 12.0˚E 14.0˚E 16.0˚E 18.0˚E 20.0˚E

59.0˚N 60.0˚N 61.0˚N

(c)

Depth: 4−7km

0 50 100 150 200 250 300 350 400

Y (km)

0 50 100 150 200 250 300 350 400 450 500 550 X (km)

10.0˚E 12.0˚E 14.0˚E 16.0˚E 18.0˚E 20.0˚E

59.0˚N 60.0˚N 61.0˚N

0.88 0.90 0.92 0.94 0.96 0.98 1.00 1.02 1.04 1.06 1.08 1.10 1.12

Relative velocity

(d)

Depth: 4−7km

0 50 100 150 200 250 300 350 400

Y (km)

0 50 100 150 200 250 300 350 400 450 500 550 X (km)

10.0˚E 12.0˚E 14.0˚E 16.0˚E 18.0˚E 20.0˚E

59.0˚N 60.0˚N 61.0˚N

(e)

Depth: 4−7km

0 50 100 150 200 250 300 350 400

Y (km)

0 50 100 150 200 250 300 350 400 450 500 550 X (km)

10.0˚E 12.0˚E 14.0˚E 16.0˚E 18.0˚E 20.0˚E

59.0˚N 60.0˚N 61.0˚N

0.88 0.90 0.92 0.94 0.96 0.98 1.00 1.02 1.04 1.06 1.08 1.10 1.12

Relative velocity

(f)

Depth: 4−7km

0 50 100 150 200 250 300 350 400 450

Y (km)

0 50 100 150 200 250 300 350 400 450 X (km)

12.0˚E 14.0˚E 16.0˚E 18.0˚E

56.0˚N 57.0˚N 57.0˚N 58.0˚N 59.0˚N

(g)

Depth: 4−7km

0 50 100 150 200 250 300 350 400 450

Y (km)

0 50 100 150 200 250 300 350 400 450 X (km)

12.0˚E 14.0˚E 16.0˚E 18.0˚E

56.0˚N 57.0˚N 57.0˚N 58.0˚N 59.0˚N

1.64 1.66 1.68 1.70 1.72 1.74 1.76 1.78 1.80

Vp/Vs ratio

(h)

Depth: 4−7km

0 50 100 150 200 250 300 350 400 450

Y (km)

0 50 100 150 200 250 300 350 400 450 X (km)

12.0˚E 14.0˚E 16.0˚E 18.0˚E

56.0˚N 57.0˚N 57.0˚N 58.0˚N 59.0˚N

(i)

Depth: 4−7km

0 50 100 150 200 250 300 350 400 450

Y (km)

0 50 100 150 200 250 300 350 400 450 X (km)

12.0˚E 14.0˚E 16.0˚E 18.0˚E

56.0˚N 57.0˚N 57.0˚N 58.0˚N 59.0˚N

0.88 0.90 0.92 0.94 0.96 0.98 1.00 1.02 1.04 1.06 1.08 1.10 1.12

Relative velocity

(j)

Depth: 4−7km

0 50 100 150 200 250 300 350 400 450

Y (km)

0 50 100 150 200 250 300 350 400 450 X (km)

12.0˚E 14.0˚E 16.0˚E 18.0˚E

56.0˚N 57.0˚N 57.0˚N 58.0˚N 59.0˚N

(k)

Depth: 4−7km

0 50 100 150 200 250 300 350 400 450

Y (km)

0 50 100 150 200 250 300 350 400 450 X (km)

12.0˚E 14.0˚E 16.0˚E 18.0˚E

56.0˚N 57.0˚N 57.0˚N 58.0˚N 59.0˚N

0.88 0.90 0.92 0.94 0.96 0.98 1.00 1.02 1.04 1.06 1.08 1.10 1.12

Relative velocity

(l)

Depth: 4−7km

Figure 10: Map views of checkerboard models and their corresponding inversion results

(right) at the depth of 4 – 7 km for the (a)-(f) central study area (on previous page), and the

(g)-(l) southern study area respectively. From top to bottom, the maps are showing V p/Vs

ratio, S-velocity model, and P-velocity model respectively.

4.4 Moho depth modelling

Any first arrival with offset more than the crossover distance in the previous crustal mod- ellings was discarded at the crossover distance, which prevented any P n or S n phase from the

"contaminating" the crustal model. P _n and S _n phases are first arrivals in the traveltime curves after the crossover distance. They are head waves refracted from the Moho boundaries. In this section, I try to resolve the Moho depth within the study areas with offset thresholds now set to 450 km to include P _n and S _n phases.

Initially, I tried to invert the velocity models in the same way as in the crustal modelling, but with the inclusion of P n and S n phases. The final models could reach as deep as 60 km.

However, with 3 km thick model cells, the resolution is too low to indicate any sharp Moho boundary. The model cells are now re-sampled to 40 × 40 × 1 km ³ . Instead of starting the inversion from the 1-D starting models, the Moho modellings start with the 3-D resultant crustal models from the previous inversions. The crustal structures are inverted for two extra iterations (with offset thresholds of 120 km) before the actual Moho modellings. The starting source locations are the relocated locations from the previous modellings, but they are fixed throughout the modellings. A possible Moho region between 40 and 50 km depth was decided for each of the study areas. Only the model cells within this depth range were allowed to change in the inversions, which means the upper-crustal models are fixed during the Moho modellings. The models are constrained similarly as the crustal modellings with different parameters. One way to increase the velocity gradients between crust-mantle boundaries is to get a priori information and "block" out a certain velocity range in the models. The blocking range of 7.5 to 7.8 km s ^{− 1} P-velocities is chosen by trial and error of different approximate ranges as previous studies (e.g. Guggisberg et al., 1991) has shown absent of this P-velocities range. S-wave velocities, scaled by the respective V _p /V _s ratio in the study area, are similarly blocked out.

Each study area takes nine iterations to get a stable model. For the central study area,

the numbers of P- and S-traveltimes are reduced from 1020 to 586 (42.5%), and from 405 to 264 (34.8%) respectively. The P- and S-residual RMS decreased from 0.563 to 0.558 (0.4%), and from 1.013 to 0.868 (14.3%) respectively. For the southern study area, the numbers of P- and S- traveltimes are reduced by 30.8% (1837 to 1272) and 37.5% (550 to 344)respectively, and the P- and S-residual RMS are reducd by 3.9% (0.643 to 0.618) and 4.9% (1.076 to 1.023) respectively.

The larger number of traveltimes in the Moho modelling than in the crustal modelling is due

to the fact that more phases are included from larger offsets. The larger residual RMS are the

results of larger residual thresholds for both P- and S- travel times (3.0 and 3.8 respectively),

for longer travelling distance of P _n and S _n phases.

5 Results

5.1 Three-dimensional crustal models

The results of the V p /V s ratio distributions and the P and S velocities (V p and V s ) relative to the starting models are presented on horizontal slices through the depth of 4 to 7 km in Fig. 11, and on vertical cross-sections along W-E and N-S oriented profiles in Fig. 12. The V p anomalies in the central area distribute on the south-eastern part of the area with low relative velocities. However, the V _s anomalies are showing relatively high velocities in the in the north-eastern part of the area. The V p /V s ratios, as a result, show a large low anomaly region in the east coast of the central area. In the southern study area, a region of low V p

occurs in the north-eastern part of the study area, but there is no V _s anomalies in that region.

This results in low V p /V s ratios (ca. 1.66 – 1.69) anomalies.

Fig. 12a shows vertical cross-sections of the central study area along a W-E and a N-S profile respectively. The W-E profile shows a high-anomaly of V _p /V _s ratio (1.73 to 1.76) from about 100 to 200 km along the profile from depth of about 20 to 40 km, and a low V p /V s

ratio region from 200 to 430 km in the upper-crustal structure with ratios from 1.62 to 1.68 with depth ranges from 5 to 15 km. The N-S profile shows the low anomalies of V _p /V _s ratios throughout the profile in the upper-crust. The cross-sections of the velocity structures in the southern study area show less variations in V _p /V _s ratios along the upper-crust than the central study area, and a generally higher ratio in lower-crust than the upper-crust.

Fig. 13a and b show the cross-sections of the absolute P-velocity structures along the B-

B’ and D-D’ profile of the central and southern study area respectively. The profiles follow

roughly the FENNOLORA survey line (Guggisberg et al., 1991), of which the results are shown

in Fig. 13c. The upper-crustal boundary is at about V p 6.5 km/s. The boundary depths of

20 to 25 km are similar to the FENNOLORA study. However, the southernmost part of the

southern profile is poorly covered in depth because of the lack of receivers in the south coast.