Modelling short-period crustal phases at regional distances for the seismic source parameter inversion

Physics of the Earth and Planetary Interiors, 47 (1987) 159—178 Elsevier Science Publishers B.V., Amsterdam — Printed in The Netherlands

159

Modelling short-period crustal phases at regional distances for the seismic source parameter inversion Won-Young Kim Seismological Department, Uppsala University, Box 120, 19,750 12 Uppsala (Sweden) (Received January 27, 1986; revision accepted September 12, 1986)

Kim, W.-Y., 1987. Modelling short-period crustal phases at regional distances for the seismic source parameter inversion. Phys. Earth Planet. Inter., 47: 159—178. Waveforms of short-period crustal phases observed at regional distances are modelled in an attempt to retrieve the source parameters of small earthquakes. For a simple shield type crustal structure, the synthetic seismograms calculated for elementary source mechanisms show that the waveforms of crustal phases are quite distinct for different source mechanisms, and suggest that these features can be used to constrain the source mechanisms. Waveforms of the crustal phases are also strongly affected even by small variations in the crustal thickness as well as in the focal depth. A generalized, linearized waveform inversion technique is used to invert the regional crustal phases for a deviatoric moment tensor source. The inversion technique is applied to regional short-period records from a small earthquake in central Sweden for determining its source parameters. The waveform inversion, which mainly used Lg waves, produced a satisfactory fit between the observed and synthetic records, suggesting that the source mechanism and focal depth estimates were reasonable. The focal depth of 40 km obtained for the event studied suggests the occurrence of a deep crustal earthquake and demonstrates the usefulness of crustal phases in studying small intraplate earthquakes.

1. Introduction Determining the source parameters of small regional earthquakes (local magnitude generally less than 5) is an important problem. In many intraplate regions with low seismicity worldwide, these earthquakes provide the only clue to the tectonics of a region. Further, a widespread geographic occurrence of these events provides valuable data to constrain regional seismic models. However, the events in qucstion are usually recorded only by a few short-period seismograph stations at regional distances. Thus, it is difficult to determine the source parameters of small earthquakes, due to the paucity of high-quality data as well as to the difficulties in interpreting the regional short-period records. Seismograms from nearby earthquakes (i~< 100 kin) often exhibit relatively simple waveforms and P and S waves are easy to discern (see Fig. la). As the epicentral 0031-9201/87/$03.50

© 1987 Elsevier Science Publishers By.

distance increases, the seismograms become more complex and the characteristic crustal phases, such as Pn, Pg, Sn, Sg and Lg, are usually observed at regional ranges (Fig. ic). However, these phases no longer arrive as discrete impulses but consist of successive arrivals distributed over intervals of a few seconds to a few tens of seconds (see Figs. 1 and 2). In recent years, our understanding of the excitation and propagation of the short-period crustal phases at regional ranges (~ = 100 1500 km) has significantly increased, primarily due to the advent of theoretical methods in calculating complete synthetic seismograms (see, e.g., Kind, 1979; Kennett, 1980; Bouchon, 1981), as well as to the advances in digital three-component seismographs. Therefore, it now seems possible to model the entire seismic wave train on regional short-period records for retrieving the source parameters, at least for a horizontally stratified —

160

(~)

UME, 78km

KJF,

-~*~ 1 ~°~t~-;~ / / /~ 0

(b)

496 km

06

~___

0.5

P

S

/ ~ I

ER

\

~2~/KJF~

/‘~

Sn

Sg

(~NUR~1\ )

upp~”.~~

3~’

~

500 km

UPP,

~(d)

(C)

492km

05

-

1 2

06

To

s,.,

~.

4.. ~~

3

—, —

~

~IiM4~Pn~

60 Sec

Fig. 1. Typical observed three-component seismograms from the three largest earthquakes in Sweden since 1962, and the locations of epicentres and recording stations. In the map, the stars indicated by nos. 1, 2 and 3 represent epicentres of the events on 28 September 1962 (ML = 4.0), 29 September 1983 (ML = 4.1) and 15 June 1985 (ML = 4.6). Propagation paths corresponding to seismogram components shown are indicated by the dashed lines. In each three-component seismogram pair, the peak to peak trace amplitude of each component normalized to T-component is shown at the end of each trace. (a) short-period seismograms at WWSSN station Umeâ (UME) from event 1, (b) short-period seismograms at WWSSN station Kajani (KJF) from event 2, (c) broad band digital recording (Wielandt—Streckeisen type) at Uppsala (UPP) from event 3, (d) short-period seismograms at DWWSSN station Bergen (BER) from event 2 for the Lg wave portion (earlier arrivals were not triggered to be recorded).

seismic velocity model (see, e.g., Vered and BenMenahem, 1974; Bouchon, 1982). In this study, the effects of the source mechanism, focal depth and crustal thickness on the waveforms of short-period crustal phases are demonstrated using the synthetic seismograms. A generalized linear inversion technique is used to retrieve the source Darameters of small earthquakes by inverting the waveforms of the crustal phases. An attempt has been made to invert regional

short-period records from a small earthquake in central Sweden for the source parameters. For the sake of brevity, ‘regional record’ is henceforth used to indicate a conventional WWSSN short-period seismogram recorded at regional ranges and ‘crustal phase’ is used to indicate a short-period crustal phase on the regional records. Z-, R- and T-component is used to indicate vertical-, radial- and transverse-component of the regional records, respectively.

161

Solberg, Sept. 29, 1983 (Z)

Lg

442km

KIR

466km

KJF

496km

AP~1WMA~Vt~M

NUR A/.~vnJ ~

0

~

520km

3Osec

Fig. 2. Observed Z-component seismograms at stations in Fennoscandia from event 2 shown in Fig. 1. The epicentre and station locations are depicted in Fig. 1. Seismogram traces are plotted with a reduction velocity of 8.12 km s~ and each trace is normalized to the maximum trace amplitude. Trace peak to peak amplitude (in millimetres) for the WWSSN short-period seismograph (magnification of 25 K) is shown at the end of each trace. The dashed line indicates a group velocity of 3.6 km ~ corresponding to Lg wave arrivals. Note that on all traces, Pn, Pg and Sn phases are hardly visible, while Lg waves dominate. The waveforms of Lg waves are relatively simple and consist of a set of distinct arrivals separated by about 8 to 10 s.

2. Modelling crustal phases In the past few years, several methods have been proposed to synthesize the entire seismic wave train from a source embedded in a layered half-space observed at regional ranges. These indude extended reflectivity (Kind, 1979), discrete wavenumber summation (Bouchon, 1981) and reflection matrix methods (Kennett, 1980), among others. Applications of synthetic seismogram techniques showed remarkable results. For example,

using the complete synthetic seismograms computed for distances of up to 500 km, Bouchon (1982) and Campillo et al. (1984) have shown that the Lg waves are guided waves made up of S waves supercritically incident on the Moho and multiply reflected within the crust. In this section, we will consider the effects of the source mechanism, focal depth and crustal thickness in generating the crustal phases by using synthetic seismogram calculations.

162

the synthetic record section from an event with a vertical strike-slip (SS) mechanism, while Figure 4 shows the record section from a vertical dip-slip (DS) mechanism event. Both sections are calculated for a receiver at an azimuth of 450 and for a focal depth of 30 km in a two layered crustal model with moderate level of attenuation (see Table Ib; Leong, 1975).

2.1. Source mechanism

To understand the role of the source mechanism in generating the crustal phases, we calculate the Z-component synthetic seismograms for epicentral distances of up to 900 km and for two elementary source mechanisms using the extended reflectivity method (Kind, 1979). Figure 3 shows T I

I

I

0

I

I

-

RANGE/7.84

I

I~

I

~

I

I

I

I

I

I

22

10.0

t*fff__ 300

I

60

(SEC) I

lBS

~0

364

S

2.26

STRIKE

SLIP

0.66

755

026

0.17

PnPgSnLg

10 13

Fig. 3. Z-component synthetic record section from an event with a vertical strike-slip mechanism at a focal depth of 30 km. The seismograms are calculated for the ground displacements convolved with the WWSSN short-period seismograph response and for a cutoff frequency of about 2 Hz. Traces are plotted with the maximum amplitude indicated at the end of each trace and the peak to peak amplitude of each trace is normalized to that at a range of 100 km, which has been assigned a value of 10. Vertical radiation patterns of a double-couple point source along the direction of receiver (450) for both the P and SV waves are depicted at the upper right corner where Z 0 indicates the source level (horizontal). Pn, Pg, Sn and Lg arrivals are indicated by the dashed lines.

163 T I I

I

0

ISO

3ØØ

-

RANG E/7.84

(sEC) I

1111111!

66

.~~~————

I

I

I

I

26

10 0

4~._J~f 4(~._._.

~

86

~

D I P

1 89

S L P

111

w

~

0.92

___________

0. 24

Fig. 4. Z-component synthetic record section from a vertical dip-slip mechanism event. Note that the trace amplitudes given at the end of each trace should be multiplied by a factor of about 1.6 to be directly comparable to those in Fig. 3. Notation is the same as in Fig. 3.

It is obvious from the initial look at the record sections of Figs. 3 and 4, that the seismograms from the SS event are more complex and most of the P and S wave trains are well developed (Fig. 3). Seismograms from the DS event are relatively simple and Lg waves are the only prominent phases (Fig. 4). If we compare the seismograms of both sections in detail, we notice that the waveforms of the crustal phases from the SS and DS events are quite different. When scanning across the record sections, we note the following differences,

Pg waves have significant amplitudes for the SS event (Fig. 3), while they are weak and hardly visible for the DS event (Fig. 4). The Pg consists of the first few crustal reflections of the P and SV-P converted waves which are leaving the source near the horizontal (see, e.g., Langston, 1982; 01sen et al., 1983; Haddon, 1984). Because of the P wave vertical radiation pattern lobes (see Figs. 3 and 4), strong reflections are produced by the SS event, and weak reflections by the DS event. It is noted that the waveforms and amplitudes of the P

164 TABLE I Crust and upper mantle velocity models used in calculations Thickness V V Density 3) (km) (km 5_i) (km s~) (g cm (a) Single layer crust 35 6.30 3.65 2.9 8.20 4.70 3.3 (b) Two layer crustal model for Sweden 19 6.22 3.58 2.7 19 6.64 3.69 2.9 7.84 4.55 3.1

—

Q i’

Qs

2000 4000

1~ 2000

500 1000 1500 (c) Crust and upper mantle model for the Baltic shield 12 6.07 3.59 2.7 500 18 6.51 3.78 2.9 900 14 6.64 3.83 3.1 1300 26 8.03 4.64 3.3 1800 8.30 4.79 3.5 2500

—

300 400 900

—

300 600 900 1200 1500

___________________________________________ wave train (between Pg onset and Sn arrival) show insignificant change with increasing distance when compared with later S wave arrivals (see Fig. 3). This feature indicates the inefficiency of Pg propagation in the crustal wave guide because of the relatively strong P— SV conversion of P energy upon incidence at the free surface and continuous leakage of P energy into the upper mantle in the form of S waves at the Moho (see Haskell, 1966; Olsen et al., 1983). The S waves between the Sn onset and the Lg arrival have significant amplitude for the SS event (Fig. 3), while they are small in amplitude for the DS event (Fig. 4). These S waves are mainly formed by the initial P waves leaving the source near the horizontal and being subsequently converted to SV waves as the waves propagate through multiple reflections in the crust. Again, because of the P radiation lobes, the SS event yields larger amplitude S waves than the DS event, For the SS event, the amplitude maximum of Lg appears to shift toward the later part of the trace as the distance increases (Fig. 3). Since, as the distance increases, successive higher multiples become prominent due to their take-off angles at the source approach to the SV radiation lobe. In case of the DS event, due to its SV radiation lobe at the horizontal (Fig. 4), most of the individual reflections become stronger as the distance increases. However, still the first few multiples

dominate within the ranges considered, due to the relatively deep focal depth used in this example. Thus, the first few multiples arriving at the earlier part of the Lg wave train have the maximum amplitude throughout the record section (Fig. 4). The companson between two elementary source mechanisms indicates that the relative amplitudes and waveforms of the crustal phases on regional records are quite different for the SS and DS events and that these features can be used as .

constraints in searching for the best mechamsms.

fit

source

2.2. Focal depth

To understand how the focal depth modifies the waveforms of the crustal phases, synthetic seismograms are calculated for various focal depths at a range of 440 km. Z-component seismograms calculated for a single layer crustal model (Table Ia) are shown in Fig. 5. In the figure, P 2 and ~2 denote the P and S waves reflected once at the Moho. Subsequent subscripts refer to multiple reflections at the Moho and at the free surface and the index refers to the number of legs of the corresponding geometrical ray. Note that with a single layer crustal model used, the successive arrivals on the synthetic traces are discrete and regularly separated in time. In Fig. 5, the synthetic trace for the focal depth of 2 km shows several distinct groups of arrivals. The first group (A) following Pn by about 10 s consists of Pg, P2, P3 and sP2 phases. The second group (B) may be identified as being made up of pPS, PS and other P—SV or SV—P converted phases (see Haddon, 1984). As the focal depth increases, phases starting upward at the source tend to arrive later (e.g., sP2), while the downward propagating phases arrive earlier (e.g., P2). Thus, the upward and downward propagating phases in the same group start to separate from each other, and at intermediate focal depths they appear as distinct arrivals (see the trace for 20 km depth). At greater depths, the upward propagating phases eventually merge into the downward propagating phases of the following group and form a new group of arrivals (C in the trace for 30 km). The Lg wave portions on synthetic traces show a similar pat-

165

~2,3

545

S

67

S8

2 ~

10 0

5 ~

W

~

10

~

15

______________k4_~e

20

__

44

5 .8

6.6

~_~4

mn~1

6.1

62

H 25

~

66

30

____________4r.__)~I,t~~

6 5

~ Pn

S2 S34

Pg ‘2

8 0 ~5,6

S79

20 SeC I

I

Fig. 5. Z-component synthetic seismograms calculated for a vertical strike-slip event at various focal depths. The solid line segments indicate the downward propagating waves at the source and the dashed line segments indicate upward propagating ones. Group of arrivals in the Lg wave train, e.g., S34 denotes S3 and S4 phases, etc. The traces are plotted with the maximum amplitude shown at the end of each trace and the peak to peak amplitude of each trace is normalized to that for a focal depth of 2 km, which has been assigned a value of 10.

tern. It is seen in Fig. 5, that the waveforms of the entire seismic wave train are strongly affected by the focal depth, and that the shallow as well as the deep sources tend to produce relatively fewer groups of arrivals, and consequently the corresponding seismograms appear to be simpler.

2.3. Crustal thickness

In Fig. 6, we present an observed T-component record and corresponding synthetics for various crustal thicknesses. The observed record is from the station Kajani (KJF) in Finland from a small

166

2:

43

m

~

~

km/sec

4

5

67

8

9

-

2.53

~

3

V,

~6O

1.69

loseC

Fig. 6. Observed T-component seismogram (top trace) at the distance of 496 km and corresponding synthetic seismograms (bottom traces) for various crustal thicknesses. The synthetic traces are for a focal depth of 40 km in the crust and upper mantle model shown in Fig. 7. The solid line segments indicate the phases with constant arrival time, while the dashed line segments denote the phases with time delays for increasing Moho depth. The traces are normalized with the maximum amplitude shown at the end of each trace. Note that the trace peak to peak amplitude decreases with increasing Moho depth.

80

100

VS

VP

-

Fig. 7. Crust and upper mantle velocity model for the Baltic shield (Table Ic) used in the synthetic seismogram calculations shown in Figs. 6—12. V~,and V~denote the P and S wave velocities, respectively.

constant arrival times for various Moho depths and are indicated by the straight line segments in the figure. Phases refracted at Moho arrive with a

event in central Sweden (see Fig. 1). The synthetic seismograms are calculated with the discrete wave-

little time delay for deeper Moho (e.g., Sn and sSn). As the Moho depth increases, the time delays of the arrivals with a larger number of Moho reflections increase due to successively increasing length of the propagation paths. For example, the phases S5 and S7 which are reflected twice and

number summation method (Bouchon, 1981) for a crust and upper mantle model for the Baltic shield (Noponen et al., 1967; Der and Landisman, 1972). The model used is given in Table Ic and is also depicted in Fig. 7. The various crustal thicknesses are represented by placing the Moho at different depths while keeping other model parameters unchanged. Note that in Fig. 6, there are phases with no time delays for increasing Moho depth. These phases start as the upward propagating phases at the source and arrive at the receiver without being reflected at Moho (e.g., Sg). Thus, the phases show

three times, respectively, at the Moho show relatively large time delays for increasing Moho depth (Fig. 6). The amplitude of Lg waves decreases by about 70% for deeper Moho (46 km) when compared with the shallow one (43 km), due to the increased distance from the Moho interface. It is seen in Fig. 6 that rather simple observed waveforms, especially the separation between the arrivals S5 and S7, can be best fit by a synthetic trace calculated for a thick crust of about 46 km in thickness. This example demonstrates the important role of the total thickness of the crust in developing the waveforms of the crustal phases.

167

2.4. Lateral heterogeneity

A suite of synthetic seismograms are calculated to assess the effects produced by the lateral heterogeneity in the crustal structure upon the propagation of the crustal phases. The effect of signalgenerated noise associated with the wave propagation in a nonstratified medium is simulated by introducing a random element in the transform domain as in the work of Kennett (1985). The complex amplitude of each frequency-wavenumber component has been multiplied by 1 + c( r1 + ir2) where r1, r2 are random numbers with a uniform distribution on the interval (0, 1) and c is

a variable constant. The simulated synthetic traces are shown in Fig. 8 for both the low- and high-level of signal-generated noise, which are indicated by the perturbation constant c 0.1 and c 0.5, respectively. The synthetic traces correspond to those in Fig. 9 and 11, except that the thickness of the second layer has been reduced from 18 to 13 km (Fig. 7). In case of the high-level of noise, the =

Ps

r

Z .

~

=

amplitude may vary up to 50% and the phase by a quarter period, for each point in transform domain. Thus, the resulting synthetic traces can indude significant deviations from the unperturbed case (see Fig. 8). By following the traces shown in the figure, it may be noted that the small amplitude phases such as, Pn and Sn, are completely masked with noise, when the signal-generated noise is increased, while, the large amplitude Lg waves can still be recognised. Although, it is a simplified procedure, the resulting perturbed seismograms show a close resemblance to the appearances of actual records (e.g., see Fig. 2). If we compare the records in Fig. 2 with those in Fig. 8, we notice that the observed records from stations UPP and KJF are similar in nature to those with the low-level of noise, while the record from NUR is comparable to the one with a high-level of noise. It may be an indication of a relatively homogeneous crustal structure in the Baltic shield. At the same time, it may also suggest that there is a substantial degree of lateral heterogeneity existing along different propagation paths within such a relatively homogeneous region.

Sg

392

3. Moment tensor inversion of regional data c~ci _n____l~_

____~__++4#14n~~~_422

The main purpose of modelling the crustal

~

as source mechanism, seismic moment and focal phases is to determine the source parameters, such depth. In recent years, moment tensor representation of the seismic wave excitation by Gilbert T form inversion methods to determine the source (1970) has been extensively used in various wave~ parameters (e.g., Gilbert and Dziewonski, 1975; C~C 1 Mendiguren, 1977; Stump and Johnson, 1977; Langston, 1981; Dziewonski et al., 1981; DoomC~0.5 0.69 bos, 1982). Even though Wallace et al. (1981) and loses Lg Wallace and Helmberger (1982) successfully inverted the long-period P wave portion of regional Fig. 8. Simulation of the effect of signal-generated noise on the records for the source mechanism, the short-period crustal phases. The traces are calculated for the moment tensor crustal phases have not yet been used in such an component Myz and for a source at a depth of 40 km in the inversion scheme. model depicted in Fig. 7. In both the Z- and T-component In this section, the generalized inverse techC’0.5

shown, the unperturbed trace is shown at the top and corresponding perturbed traces are shown below. The level of noise is indicated by the parameter c.

nique along with the moment tensor formalism for a point shear-dislocation is briefly discussed to

168

utilize the technique for determining the source parameters from regional records. Our general strategy is to use simplest possible parameterization of the seismic source wavefield. 3.1. Moment tensor formalism

Using the moment tensor representation, a component of ground displacement in cylindrical coordinates excited by a point source may be written as (e.g., Dziewonski et al., 1981) 6

u(r, z,

~,

t)

=

~

G~(r, z,

~,

t)* m

1(t)

(1)

f—I

where m3 represent the six independent components of moment rate tensor, Gj are the Green’s functions corresponding to the individual components of the moment rate tensor and asterisk denotes convolution. The moment rate tensor in ________________________

_________________________

______

-~

496km

eq. 1 consists of an isotropic and a deviatoric part and is defined by the six independent time histories. To limit the number of free parameters in eq. 1, we first consider that the source is purely deviatoric. This approach is justified since we are studying tectonic earthquakes with no net volume changes, which have been modelled by a shear-dislocation with great success (Engdahl and Kanamori, 1980). For such a purely deviatoric source, a linear constraint that the associated moment tensor have zero trace is easily imposed. Further, we assume that all the components of the moment rate tensor have a common source time function and that the form of the time function is known. This assumption may be justified in seismic source studies, since the Earth structure is not perfectly known and the data always have errors. Thus tradeoffs will occur between source and structure through the errors in the data and in the model. This does not warrant us to use a more complicated source model.

Obs

R

Lg

496km

PS

-

4.68

0 b s. Lg

r IA

lAss

~

236

~

246

0.41

F

Mss ~

Msz

~

~

-~_.~

0.33

Ms~

3.81

Msz

Fig. 9. Z-component seismogram observed at station KJF (top trace) and corresponding Green’s functions calculated for the components of deviatotic moment tensor (bottom traces). The station is at the range of 496 km with the azimuth of 85° and the source is at a depth of 40 km in the crust and upper mantle velocity model (Fig. 7) with the Moho depth of 45 km. Each synthetic trace consists of 1024 points and has a cutoff frequency of about 3 Hz. Peak amplitudes are given at the end of each trace,

M~

28 023

~

~

0.’B

~

~

200, .

~.

0SeC

Fig. 10. Observed R-component seismogram (top trace) and corresponding Green’s functions (bottom traces). Note the relatively large P wave amplitudes on the synthetic traces compared with those on the observed trace. Notation is the same as in Fig. 8.

169

By factoring out the common source time function from the moment rate tensor in eq. 1, a component of ground displacement for a point shear-dislocation (i.e., deviatoric moment tensor) source can be written as (cf. Stump and Johnson, 1977) 5

u(r, z,

~,

t)

~ [G’(r,

=

~,

t) * s(t)] ~

~,

j=I

(2) where s is the source time function, G,1’ are the new Green’s functions corresponding to the individual components of the deviatoric moment tensor (see Langston, 1981), and M~represent the five components of the deviatoric moment tensor (M~ Mxx, M2 Myy, M3 Mxy, M4 Mxz, M5 Myz), which are explicitly given in Aki and Richards (1980, p. 117). Our moment tensor inversion is based on eq. 2. The Green’s functions may be computed using any available method. In this study discrete wave=

=

=

=

=

~

number summation (Bouchon, 1981) and extended reflectivity (Kind, 1978, 1979) methods were used. Example calculations of the Green’s functions for the five components of the deviatoric moment tensor (eq. 2) are shown in Figs. 9, 10 and 11 for the Z-, R- and T-component, respectively. The Green’s functions are calculated for a focal depth of 40 km and are shown for the ground displacements convolved with the WWSSN short-period seismograph response. The source time function used is a smooth ramp type time dependence with a rise time of 0.2 s. 3.2. Generalized inverse method The objective of the inversion is to determine the components of the deviatoric moment tensor given in eq. 2. The observed data are inverted by using the generalized, linearized inverse technique (Wiggins, 1972; Barker and Langston, 1982), in which a set of parameters that minimize the differences between the observed and synthetic seismograms is found in a least-squares sense. The problem, in linearized form, can be written as (Wiggins, 1972) wher

I I

I v

M~r

~

~-—

_____

~

—

1.08

1.08

~

6Ixs

12.30 I

~

_____-.

4.88

II

--

Msz .. ~.

--. ~.

synthetic, C1, seismograms, A is the matrix of partial derivatives (A. aC./8P.) of the synthetic J seismogram with respect to each parameter, P~, and 1~p is the vector matrix of the unknown parameter corrections to be solved. The matrix A can be factored into the product . =

I I

. -

sthe vector matrix containing the dif-

.

.

I

A UAV’ where A is the diagonal matrix of the non-zero singular-values of A, U is the matrix of eigenvectors associated with the columns (data) of A, V is the matrix of eigenvectors associated with the rows (parameters) of A, and denotes the transpose. The parameter corrections are obtained from the generalized inverse =

‘

.~.

~

~.

___________

Fig. 11. Observed T-component seismogram (top trace) and corresponding Green’s functions (bottom traces). Note that the Green’s functions for the moment tensor components, Mxx and Myy have identical waveforms and amplitudes but differ by their polarities. Notation is the same as in Fig. 9.

tsp

=

At tsc

=

V A’U’ tsc

(4) .

.

where denotes the generalized inverse which simultaneously minimizes IA tsp — tsc 2 in the data space and I tsp 2 in the parameter space

170

(Wiggins, 1972). The parameter corrections obtamed in (4) are added iteratively to the starting model parameters and the inversion is repeated until some termination criterions are met. The convergence of the inversion process can be checked by the least-squares error I A tsp tsc I 2 A numerical measure of the fit is given by the rms error 2 rms error (1/N) (0 C )2 1/ —

=

~

I —

where N is the number of the data points used. The umqueness of the solution is measured through a parameter resolution matnx formed from V, the matnx of ej~envectorsassociated with the parameters .

R

=

V V’

.

(5)

If R is an identity matrix, the solution is unique and each parameter is perfectly resolved. The deviation of R from the identity matrix indicates non-uniqueness of the solution and the locations of non-zero off-diagonal elements indicate which parameters are trading off in the inversion (Wiggins, 1972; Barker and Langston, 1982). To assess the resolution in data space and to estimate information contents among the data, another resolution matrix called ‘information density’ matrix (Jackson, 1972), also known as ‘hat’ matrix in the statistical literature (Cook and Weisberg, 1982), is formed by using the eigenvector matrix u ‘6’ H UU ~ The rank (H) m, and if m k, where k is the number of independent variables, then this condition guarantees the existence of a solution (see Jackson, 1972). The diagonal element h., of the matrix H is a number such that 0 h 1 1 and has a mean value of rn/N. Also the sum of the diagonal elements h,, m. The diagonal element h11 can be used as a quantitative measure of the relative contribution that the i-th datum makes to the solution. Minster et al. (1974) referred to it as the ‘importance of the i-th datum’. This property of h1, lies in the fact that the information density matrix H does not depend on the actual values of the data. =

=

=

~

=

Moreover, h.1 is an additive quantity, thus ‘the importance of a subset of data is simply the sum of the importances of the individual data’ (Mmster et al., 1974). The reliability of the solution is measured by its covariance matrix, provided that an inverse of A exists (see Jackson, 1972; Aki and Richards, 1980) 2VA2V’ (7) Cov(tsP)=o where o2 errors is the(Wiggins, variance of randomly distributed, problem 1972). Assuming that all the components of the data share the same vari2 ance, a in (7) is substituted by the sample variance of the data. The diagonal elements of the covanance matrix are the vanances of parameters resolved and the off-diagonal elements are the covanances between these parameters. Thus, it allows a quantitative estimate of errors and the trade-off between these parameters. The deviatoric moment tensor obtained from the inversion is decomposed into a double-couple and a compensated linear vector dipole (CLVD) as the remaining non-double-couple component (Geller, 1976). The percentage contribution of the CLVD component is used as a measure of error in the moment tensor inversion due to unmodelled effects, such as source complexity and imperfect crustal model. 3.3. Experiment with synthetic data A set of experiments has been performed to test the reliability of the inversion method and to assess the importance of different components of ground motion as well as the effect of signal-generated noise upon the moment tensor inversioll, The experiments were carried out by using the Green’s functions calculated for a known source at a depth of 40 km in the crustal model depicted in Fig. 7. The corresponding synthetic seismogram data are generated for a set of four stations listed in Table II and a low-level of noise (c 0.1) is included in the calculation (see Fig. 8). The source time function is assumed to be known and not included in the inversion, so that the generalized inverse becomes a linear least squares process. Note that the rms error given in this study has units of millimetres and its value corresponds to =

171

noise, the solution deteriorates further compared with the previous case and the rms error is increased to about 1.325 and the CLVD contribu-

TABLE II Station information Station code

Distance (km)

Azimuth (°)

Moho depth (km)

Component used

UPP KJF NUR

442 496 520

180 85 135

46 43

Z, R, T Z, R, T

the trace amplitude on a record written by the WWSSN short-period seismograph with peak magnification of 25 K (i.e., 1 mm 4 x 10—2 micron of ground motion at 1 s period). The first experiment examines the extent to which we can recover the true solution from records which are corrupted by signal-generated noise. Two trials are made, one using the synthetic data with a low-level ol signal-generated noise and the other by using the data with a high-level of noise (c 0.5). In both trials, only Z-component data from four stations are used to invert for the —

=

moment tensor. When the data with the low-level of noise is inverted, there is a very slight deterioration in the solution, which is indicated by the rms error of 0.265 and CLVD component of 3.6%. Table III lists the results of the inversions. In the table, the moment tensor elements are given together with their standard deviations and M0 serves as an amplitude factor. In case of the inversion using the data with the high-level of

tion reaches to about 15%. Nevertheless, the moment tensor components recovered are reasonably good and the double-couple mechanism obtained is virtually not changed. This experiment indicates that the inversion for the moment tensor is stable even when there is a relatively high degree of signal-generated noise in the records. The effect of the lateral heterogeneity upon the inversion may be viewed as to develop deficiency in the model (Green’s functions) and consequently to increase errors in the solution. Another trial is carried out to test the robustness of the inversion by introducing an incorrect crustal model into the inversion. For this trial, Green’s functions are calculated for a model with a crustal thickness of 43 km, while the corresponding synthetic seismograms are generated for a model with a thickness of 45 km. In this case, we were able to align only the phase with largest amplitude on the synthetic seismogram (i.e., S3, see Fig. 6) but the other phases on the synthetic seismogram were not properly aligned against those on the Green’s functions. Such a misalignment might simulate part of the effects produced by the unknown crustal model in the inversion. The inversion results are given in Table III. They may be considered as reasonable, considering the misalignments of phases. The wrong model used is reflected through a large rms error and a CLVD component.

TABLE III Inversion results for data with signal-generated noise Correct solution Mxx Myy Mxy Mxz Myz M0 (dyne-cm) rms error CLVD (%) Strike Dip Rake

0 0 0 — 0.643 0.766 1.0

— 0.0 220° 90° 270°

Perturbed (c 0.1) 0.007±0.060 —0.004±0.050 —0.034±0.096 — 0.645 ± 0.038 0.761 ±0.05 3 1.05 0.265 3.6 220.3° 89.6° 270.0°

Perturbed (c 0.5) 0.029±0.079 0.014±0.066 0.145±0.126 — 0.655 ±0.049 0.748 ±0.069 1.24 1.325 15.4 221.2° 88.2° 270.1°

Moho depth (43 km) 0.181±0.119 —0.204±0.098 0.365±0.200 — 0.514 ± 0.053 0.828 ±0.077 2.33 2.893 40.2 212.3° 83.5° 270.9°

172

The second experiment was designed to provide information about the relative importances of different phases on the record and different components of ground motion in the moment tensor inversion. In the first trial, complete three-component records from four stations are used and the P and S waves are windowed separately to assess their relative importance. Information distribution among the data considered is assessed by using the diagonal elements of the information density matrix (eq. 6) and information carried by a subset of data is obtained by summing the individual elements. For this particular station geometry used, information carried by the R-component is only about 12%, while T-component carries about 57% of the total information distribution. Information carried by all subsets of data is given in Table IV. The total contribution of P wave (between Pn onset and Sn arrival) is only about 8% (see Table IV). This indicates that most of the information is carried by the S waves, especially the large amplitude Lg waves, Three trials are carried out to understand the role of a particular component of ground mOtion in the inversion. If only Z- or R-component from all four stations are used, the correct solution is recovered with small rms errors and CLVD cornponents. However, analysis of the covariance matrix indicates that there are relatively high correlations among the moment tensor elements, in particular Mxx and Myy, and that Mxy is poorly resolved. In this case, all four stations have almost equal importance and for the R-component, the P waves carry about 50% of the total information, When only the T-component data are used, the

correct solution is not recovered and the solution found is not unique. It is due to the complete correlation between the moment tensor elements Mxx and Myy in the inversion (see Fig. 11; also Barker and Langston, 1982). Further trials were carried out to give an assessment on the redistribution of information among the data by using the partial combinations of components or stations. These trials were carried out by using only S waves. If any single Z- or R-component record is added individually to the T-component data from all four stations, the correct solution is recovered in all attempts. Any single Z- or R-component record added separately carries at least 20% of the total information, and also the importance of the corresponding Tcomponent increases relative to other stations. In case, any single T-component data are added individually to the Z-component from all four stations, similar results are obtained as in the previous case. The distribution of information among the different sets of data is given in Table V. However, if any single R-component is added individually to the Z-component from all four stations, then any R-component added carries less than 5% of the total information. Moreover, the importance of the corresponding Z-component decreases relative to the other stations (see Table V). The contribution of the R-component data to the moment tensor inversion is doubtful, when only the S wave portions are used which may be quite common in actual data set. The result of this set of experiment underlines the important role of the different components of ground motion in the moment tensor inversion.

TABLE IV Information distribution when all stations and components are used Station

Z

R (P

UPP KIR KJF NUR Sum

(%)

0.587 (0.038, 0.460 (0.030, 0.273 (0.018, 0.224 (0.011,

5)

T (P

5)

Station sum

(5)

0.548) 0.430) 0.254) 0.214)

0.231 (0.126, 0.105) 0.174 (0.092, 0.081) 0.109 (0.050, 0.059) 0.076 (0.025, 0.051)

0.846 0.735 0.551 0.735

1.544 (0.097, 1.446) 30.9 1.9 28.9

0.590 (0.293, 0.296) 11.8 5.9 5.9

2.867 57.3

=

1.664 (33%) 1.369 (27%) 0.933 (19%) 1.035 (21%)

=

5.00

= = =

173 TABLE V Redistribution of information when single component is added Trial

Z-component UPP

KIR

KJF

NUR

1 2 3 4

1.312 a 1.201 1.188 0.902

0.875 1.021 0.782 0.963

0.749 0.724 1.222 a 0.707

0.401 0.482 0.535 1.202 a

1.164 1.572 1.272 1.226

5 6 7 8

1.249 1.305 1.329 1.330

1.083 1.049 a 1.124 1.115

1.227 1.236 1.079 a 1.203

1.211 1.221 1.194 1.099 a

0.230 0.189 0.255 0.253

b

a

a

b b b b b b b b

Station

Component

added

added

UPP KIR KJF NUR

T T T T

UPP KIR KJF NUR

R R R R

Denotes importance of Z-component of the station added by a T- or R-component. Denotes the importance of a station added.

Further, it demonstrates the potential power of the information density matrix to assess the distribution of information among the data being considered. 3.4. Application to regional data The method described above is applied to regional records from the 29 September 1983 Solberg, central Sweden, event (ML 4.1) observed at four stations in Finland and Sweden. The data set consisted of three-component seismograms from stations, KJF and NUR, and the Z-component seismograms from stations UPP and KIR. The station locations are shown in Fig. 1. Epicentral distances of the stations used are in the range of 440—520 km, while azimuthal coverage of the observation is far from ideal and covered about 165°. The station information is listed in Table II. Each N—S and E—W pair of seismograms was rotated into its back azimuth to obtain the Tand R-component records. To calculate Green’s functions, the crust and upper mantle velocity model for the Baltic shield (Table Ic) is used for all stations and a moderate level of anelastic attenuation is included in the calculation using the attenuation factors for both the P and S waves. The source time function is assumed to be known and is not included in the inversion. The spectral analysis of Lg waves on the Z-component seismograms recorded at 10 seismograph stations in Fennoscandia suggested =

that the event had a source radius of about 0.7 km (Kim et al., 1985). Assuming a rupture velocity of about 2.5 km s~, the total rupture duration is 0.3 s. Thus, a ramp type source time function with rise time of about a half the total duration (about 0.2 s) was used to compute the synthetics. For the waveform inversion to be meaningful, the observed and synthetic data must be synchronized. However, it is impossible to synchronize all the crustal phases on both the observed and synthetic records, due to the uncertain source location as well as to the effects of unknown crustal structure. Therefore, we windowed the observed and synthetic records beginning at Lg wave arrivals to about 15 s, and lined up both data sets. This simple windowing procedure eliminates part of the contamination due to the imperfect modelling of the crustal structure and much of the noise present in the observed records. Although the relative amplitudes between the P and Lg waves are important to constrain the source mechanism, P waves could not be used for the inversion due to extremely small trace amplitudes below noise level (Fig. 12). Since the inversion is nonlinear with respect to the focal depth and structural parameters, the focal depth and Moho depth are constrained by a trial-and-error fit of the synthetic seismograms to the observations before the waveform inversion. By a trial-and-error fit of the T-component synthetic seismograms calculated for various focal depths, we found that the synthetics for a focal —

174

SOL BERG, (Z) Sept. 29, 1983 54. ~

~

1

442 k~

found for the station paths are in agreement with those reported by Bungum et al. (1980) for the Baltic shield. After the focal depth and Moho depths are

Obs --

~

36 6

0.5 KIR 466 km 46 en

-.

~

1

KJF

49665m - -~=

NUR

‘s49(~44*s844l

231

520 ks 9.2

Fig. 12. Comparison of observed Z-component records (lower traces) from the Solberg earthquake with synthetic seismograms (upper traces). Trace peak to peak amplitudes (in millimetres) for the WWSSN short-period seismograph (magnification of 25 K) is given at the end of the observed seismograms, while the numbers at the end of the synthetic traces indicate the ratio of the maximum trace amplitude of the synthetic traces to that of the observed ones. Sampling interval on both the observed and synthetic traces is about 0.17 s. Note that the observed record at station KIR is recorded by a Grenet short-period seismograph, while the synthetic trace is for the WWSSN short-period seismograph. The source parameters used are M 0 =1.OX1O2i dyne-cm, strike = 225°, dip = 90°,rake = 270° and the ramp type source time function with a rise time of 0.25 s.

depth of 40 km provided the best solution. Synthetic seismograms calculated for various Moho depths and focal depths of 40 km showed that the observed waveforms especially the separation between the prominent phases, S3 and S5 (Fig. 12), would be best fit by a particular Moho depth for each station path, and suggested that the crustal thickness might vary substantially to each station path. The Moho depth for each station path constrained by the best fit synthetics lies between 43 and 46 km (Table II). The crustal thicknesses

determined, the Lg wave portions of the Z- and T-component records are inverted for the moment tensor. The R-component data are not used, because their contribution is small and doubtful as discussed whenThe onlyresults the S of waves be used in theearlier inversion. threecan inversions using the Z-component from three stations (excluding KIR), all four Z-component and both the Z- and T-component, respectively, are listed in Table VI. In all cases, the uniqueness of the solution is ensured through the resolution matrix defined in eq. 5, while the reliability of the solution is assessed by analysing the covariance matrix (eq. 7). In case of the inversion using only Z components, the moment tensor elements Mxz and Myz are fairly well resolved, while Mxy is poorly resolved. Also the covariance matrix mdicates that there is trade-off between Mxx and Myy, resulting in relatively large standard deviations for these elements (see Table VI). The double-couple mechanism obtained suggests a predominantly pure thrust motion along a near vertical nodal plane for this event. The CLVD component ranges between 8 and 12%. When both the Z- and T-component are used, all the moment tensor elements seem to be fairly well resolved and the double-couple mechanism obtained is comparable to the previous inversions. However, the CLVD component is increased to over 40% and the rms error is also slightly increased (Table VI). Further trials using the T-component data suggest that the Green’s functions for the event studied are unable to satisfy both the Zand T-component data simultaneously. Although, there is apparent inconsistency between the Zand T-component data, the double-couple mechanisms obtained from three inversions are fairly consistent. A comparison of the Z-component observed records and synthetic seismograms computed for a source mechanism determined from the inversion is shown in Fig. 12. The overall agreement of relative amplitudes and waveforms of the major phases between the observed and synthetic traces

175 TABLE VI Inversion results of the Solberg event (Inversion cases a Z (3)

Z (4)

Z and T

Moment tensor (X 102i dyne-cm) Mxx 0.020±0.169 —0.035±0.185 0.005±0.071 Myy 0.092±0.136 0.181±0.148 0.166±0.063 Mxy 0.061 ±0.318 —0.258±0.343 —0.562±0.047 Mxz —0.882±0.097 —0.848±0.084 —0.703±0.073 Myz 0.468±0.112 0.440±0.116 0.539±0.111 Rms error 3.155 3.524 3.702 (nun) CLVD(%) 12.7 7.9 41.5 Strike 62.1° 62.9° 61.8° Dip 88.6° 84.7° 80.5° Rake 89.7° 103.7° 110.5° a Z(3) and Z(4) indicate the inversions using only 3 and 4 of the Z-component, respectively, while Z and T indicate inversion using both Z- and T-component data.

is fairly good, suggesting that the source mechanism and focal depth estimates are reasonable. Note that on most of the observed traces, the separation as well as relative amplitudes between the first large amplitude phases (S3) and the second ones (S~)in the Lg wave trains are reasonably well reproduced by the corresponding synthetics.

Solberg, (T)

Despite the overall satisfactory match, there are some discrepancies between the observations and synthetics. The Sg phases which also include primary Moho reflection (S2) on most of the observed records are weak and hardly discernible compared with those on the corresponding synthetics. This feature is probably due to the imperfect crustal model used and may indicate that the true S wave velocity contrast across the Moho is slightly different from the one used to generate the synthetic traces. At station NUR, the observed amplitude is smaller than the synthetic one by a factor of about three. The synthetic seismogram predicts a much larger amplitude than the observed one, since the station is at the direction of the maximum SV wave horizontal radiation pattern. Figure 13 shows a comparison between the observed and synthetic T-component traces. The synthetic traces were computed with the same parameters as those used for the Z-component synthetics. The overall agreement between the observed and synthetic traces is rather poor and only the arrival times of the prominent phases at KJF are somewhat matched. The relative amplitudes between the prominent phases on the observed traces are not very well reproduced by the synthet-

Sept. 29, 1983

Syn

KJ F

1.1 0bS. ~

37.4

Sn

NUR

Lg

6 ~6.2

0I I

I

Fig. 13. Comparison of observed T-component records (lower traces) with synthetic ones (upper traces). Notation is the same as in Fig. 12.

176

ics. Note that the observed amplitude at station NUR is larger than the synthetic one by a factor of about two, which is the opposite to the case of Z-component data discussed earlier. Although, inconsistency between the Z- and T-component data is also suggested in the inversion, such a large scatter of amplitude is difficult to explain, Anomalous crustal structure beneath this station (S. Pirhonen, personal communication, 1985) might be an explanation. Alternatively, the presence of lateral heterogeneity as well as a crustal pinch (e.g., graben structure) along the path may produce such an inconsistency, since it leads to a build-up of T-component amplitudes, due to a progressive net transfer of energy from Rayleigh to Love modes as the wave train propagates through such structure (see Kennett and Mykkeltveit, 1984).

4. Discussion and conclusions As shown in the example (Fig. 12), a satisfactory fit between the observed and synthetic regional records is obtained. We believe that the successful fit is achieved mainly due to ,the relatively homogeneous crustal structure in the Baltic shield studied which, in its turn, enables us to calculate reasonable Green’s functions. The criteria in searching for the source mechanism were based on the observations, that the waveforms and amplitudes of crustál phases excited by the events with different source mechanisms are quite distinct, at least for relatively simple crustal models (see Figs. 3 and 4). However, there are situations when such criteria may fail to constrain the source mechanisms. For example, the presence of lowvelocity sedimentary layers near to the surface can change the waveforms of the entire wave train. Especially, the Pg waves become more complex due to the increased wave guide efficiency provided by the low-velocity layer (Olsen et al., 1983). The presence of lateral heterogeneity along the path may also change the waveforms and amplitudes of the crustal phases due to the coupling between the Rayleigh and Love modes (see Kennett, 1984). Further complications of waveforms and amplitudes may also be introduced through

the localized site-response effect beneath the receiver (see e.g., Barker et al., 1981). Although a low-pass filtering may be applied to smooth part of the complications in the waveforms (see e.g., Barker and Langston, 1982), generally, it will be very difficult to use such data, since those effects can only be disentangled partially by using the few observations available. The inconsistency between Z- and T-component data (Figs. 12 and 13) suggests that the Green’s functions for the event studied are unable to satisfy both components data simultaneously. The overall good agreement between the observed and synthetic records obtained are encouraging for further studies of crustal phases. It has been shown that a satisfactory fit of synthetic seismograms to observed regional records is obtained, and that the source parameters of small earthquakes can be determined using the waveform inversion technique, at least for regions with weak lateral heterogeneities. Actually, many intraplate regions have such a relatively homogeneous crustal structure and release only small earthquakes. Thus, the crustal phases can be a useful tool for seismic source studies of small intraplate events. The focal depth of about 40 km deduced for the event studied here demonstrates the importance of these kind of studies, since it suggests the occurrence of a deep crustal earthquake, which has an important implication in the tectonics of the Baltic shield. We believe that the inversion method can be applied to small and shallow earthquakes which occur in tectonic regions where there may be considerable lateral heterogeneity in the crust. However, in those cases, it will be important to interpret the observed records correctly and to improve the method for calculating the correct Green’s functions. Numerical techniques such as a low-pass filtering may be used to remove unnecessary complications in the data and to stabilize the inversion (e.g., Barker and Langston, 1982). Further, the assessment of information contents among the data discussed in the earlier section can be extremely helpful to understand the nature of the data and to improve the stability of the inversion in such a complex case.

177

Acknowledgements I thank Ota Kulhánek, Giuliano Panza and Nitzan Rabinowitz for critically reading the text and suggesting vanous improvements to the manuscript. I also thank Brian Kennett for his helpful comments and for suggesting a technique to simulate the signal-generated noise. I am grateful to Michel Bouchon, Michel Campillo, Reiner Kind and Charles Langston for providing me with computer programs which were used at various stages of the present work. I appreciate the critical comments by anonymous reviewers. This research has been carried out at the Seismological Department, Uppsala, Sweden with a partial financial support of the Swedish Natural Science Research Council under contract G-GU 3164-132.

References Aki, K. and Richards, P.G., 1980. Quantitative Seismology. W.H. Freeman, San Francisco, Vol. 1, 557 pp. Barker, B.W., Der, Z.A. and Mrazek, C.P., 1981. The effect of crustal structure on the regional phases Pg and Lg at the Nevada test site. J. Geophys. Res., 86: 1686—1700. Barker, J.S. and Langston, C.A., 1982. Moment tensor inversion of complex earthquakes. Geophys. JR. Astron. Soc., 68: 777—803. Bouchon, M., 1981. A simple method to calculate Green’s functions for elastic layered media. Bull. Seismol. Soc. Am., 71: 959—971. Bouchon, M., 1982. The complete synthesis of seismic crustal phases at regional distances. J. Geophys. Res., 87: 1735—1741. Bungum, H., Pirhonen, S.E. and Husebye, E.S., 1980. Crustal thicknesses in Fennoscandia. Geophys. J.R. Astron. Soc., 63: 759—774. Campillo, M., Bouchon, M. and Massinon, B., 1984. Theoretical study of the excitation, spectral characteristics, and geometrical attenuation of regional seismic phases. Bull. Seismol. Soc. Am., 74: 79—90. Cook, RD. and Weisberg, S., 1982. Residuals and Influence in Regression. Chapman and Hall, New York-London, 230 pp. Der, Z.A. and Landisman, M., 1972. Theory for errors, resolution, and separation of unknown variables in inverse problems, with application to the mantle and crust in southern Afnca and Scandinavia. Geophys. J.R. Astron. Soc., 27: 137—178. Doornbos, Di., 1982. Seismic moment tensors and kinematic source parameters. Geophys. JR. Astron. Soc., 69: 235—251. Dziewonski, A.M., Chou, T.-A. and Woodhouse, J.H., 1981.

Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res., 86: 2825-2852. Engdahl, E.R. and Kanamon, H., 1980. Determination of earthquake parameters. EOS Trans. Am. Geophys. Union, 61: 62-64. Geller, R.J., 1976. Body force equivalents for stress-drop seismic sources. Bull. Seismol. Soc. Am., 66: 1801—1804. Gilbert, F., 1970. Excitation of the normal modes of the Earth by earthquake sources. Geophys. JR. Astron. Soc., 22: 223226 Gilbert, F. and Dziewonski, A.M., 1975. An application of normal mode theory to the retrieval of structural parameters and source mechanisms from seismic spectra. Philos. Trans. R. Soc. London, Ser. A: 278: 187—26?. Haddon, R.A.W., 1984. Computation of synthetic seismograms in layered Earth models using leaking modes. Bull. Seismol. Soc. Am., 74: 1225—1248. Haskell, N.A., 1966. The leakage attenuation of continental crustal P waves. J. Geophys. Res., 71: 3955—3967. Jackson, D.D., 1972. Interpretation of inaccurate, insufficient and inconsistent data. Geophys. JR. Astron. Soc., 28: 97—109. Kennett, B.N.L., 1980. Seismic waves in a stratified half space —II theoretical seismograms. Geophys. J. R. Astron. Soc., 61: 1—10. Kennett, B.N.L., 1984. Guided wave propagation in laterally varying media—I theoretical development. Geophys. J. R. Astron. Soc., 79: 235—255. Kennett, B.N.L., 1985. On regional S. Bull. Seismol. Soc. Am., 75: 1077—1086. Kennett, B.N.L. and Mykkeltveit, S., 1984. Guided wave propagation in laterally varying media—I! Lg-waves in northwestern Europe. Geophys. J.R. Astron. Soc., 79: 257—267. Kim, W.Y., Kulhânek, 0., van Eck, T. and Wahlstrom, R., 1985. The Solberg, Sweden, earthquake of September 29, 1983. Report 1-85, Seismol. Dep., Uppsala Univ., 49 pp. Kind, R., 1978. The reflectivity method for a buried source. J. Geophys., 44: 603—612. Kind, R., 1979. Extensions of the reflectivity method. J. Geophys., 45: 373—380. Langston, C.A., 1981. Source inversion of seismic waveforms: the Koyna, India, earthquakes of 13 September 1967. Bull. Seismol. Soc. Am., 71: 1—24. Langston, C.A., 1982. Aspects of Pn and Pg propagation at regional distance. Bull. Seismol. Soc. Am., 73: 457—471. Leong, L.S., 1975. Crustal structure of the Baltic shield beneath Umeã, Sweden, from the spectral behavior of long-period P waves. Bull. Seismol. Soc. Am., 65: 113—126. Mendiguren, J.A., 1977. Inversion of surface wave data in source mechanism studies. J. Geophys. Res., 82: 889—894. Minster, J.B., Jordan, T.H., Molnar, P. and Haines, E., 1974. Numerical modelling of instantaneous plate tectonics. Geophys. J. R. Astron. Soc., 36: 541—576. Noponen, I., Porkka, M.T., Pirhonen, S. and Luosto, U., 1967. The crust and mantle in Finland. J. Phys. Earth, 15: 19—24. Olsen, K.H., Braile, L.W. and Stewart, J.N., 1983. Modeling

178 short-period crustal phases (P, Lg) for long-range refraction profiles. Phys. Earth Planet. Inter., 31: 334—347. Stump, B.W. and Johnson, L.R., 1977. The determination of source properties by the linear inversion of seismograms. Bull. Seismol. Soc. Am., 67: 1489—1502. Vered, M. and Ben-Menahem, A., 1974. Application of synthetic seismograms to the study of low-magnitude earthquakes and crustal structure in the northern Red Sea region. Bull. Seismol. Soc. Am., 64: 1221—1237.

Wallace, T.C. and Helmberger, D.V., 1982. Determining source parameters of moderate-size earthquakes from regional waveforms. Phys. Earth Planet. Inter., 30: 185—196. Wallace, T.C., Helmberger, DV. and Mellman, G.R., 1981. A technique for the inversion of regional data in source parameter studies. J. Geophys. Res., 86: 1679—1685. Wiggins, R.A., 1972. The general linear inverse problem: implication of surface waves and free oscillations for earth structure. Rev. Geophys. Space Phys., 10: 251—285.