Sens Schoenfelder Dissertation

Exploring Scattered Waves in Seismology


Von der Fakultat fur Physik und Geowissenschaften

der Universitat Leipzig

genehmigte

DISSERTATION

zur Erlangung des akademischen Grades

doctor rerum naturalium

(Dr. rer. nat.)

vorgelegt

von Christoph Sens-Schonfelder

geboren am 22. November 1977 in Magdeburg

Gutachter: Prof. Dr. Michel Campillo

Prof. Dr. Heiner Igel

Prof. Dr. Michael Korn

Tag der Verleihung: 17. Dezember 2007

Dedicated to Susann

BIBLIOGRAPHICAL DATA

Sens-Schonfelder, Christoph


University of Leipzig, dissertation111 pages, 106 references, 39 figures, 4 tables

Abstract

This dissertation describes two different approaches to the use of scattered waves inseismology. In part one, envelopes of regional earthquakes are analyzed on the basisof radiative transfer theory to characterize the scattering and attenuation propertiesof the Earth’s crust and to obtain information about seismic sources. This includesa study of Lg-blockage in the western Pyrenees with elastic radiative transfer theoryin a model with deterministically described large scale structure and statisticallydescribed small scale heterogeneity. The concept of Green’s function retrieval fromambient seismic noise is used in part two. Here the method passive image interfer-ometry is introduced and verified by an application to data from Merapi volcano.It is also shown how the passively retrieved Green’s functions can be used to obtaininformation about the data quality of seismic networks.

Contents

Acknowledgments vii

Overview 1

I Radiative Transfer Theory 5

1 Acoustic radiative transfer theory 7

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.2 Envelope modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.3 Data and regional setting . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.4 Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

1.5.1 Attenuation parameters . . . . . . . . . . . . . . . . . . . . . 17

1.5.2 Source energy . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

1.5.3 Site response . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

1.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

1.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2 Monte-Carlo simulations of elastic wave energy 29

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.2 Elastic radiative transfer theory . . . . . . . . . . . . . . . . . . . . . 32

2.3 Monte-Carlo techniques to solve the radiative transfer equation . . . 36

2.4 Outline of the RadTrans algorithm . . . . . . . . . . . . . . . . . . . 37

2.4.1 The coordinate system . . . . . . . . . . . . . . . . . . . . . . 38

2.4.2 Description of a particle . . . . . . . . . . . . . . . . . . . . . 39

2.4.3 Particle emission from a point shear dislocation . . . . . . . . 39

2.4.4 Scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

2.4.5 Intrinsic attenuation . . . . . . . . . . . . . . . . . . . . . . . 48

2.4.6 Ray tracing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

2.4.7 Recording of particles . . . . . . . . . . . . . . . . . . . . . . . 51

2.4.8 Evaluation of probability distributions . . . . . . . . . . . . . 51

2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

v

vi CONTENTS

3 Elastic radiative transfer and Lg-blockage in the Pyrenees 533.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 553.2 Structure and evolution of the Pyrenees . . . . . . . . . . . . . . . . . 553.3 Data and observations . . . . . . . . . . . . . . . . . . . . . . . . . . 563.4 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.5 Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

3.5.1 The undisturbed case . . . . . . . . . . . . . . . . . . . . . . . 623.5.2 The disturbed case . . . . . . . . . . . . . . . . . . . . . . . . 62

3.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

II Phase Correlations 69

4 Passive image interferometry at Mt. Merapi 714.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.2 Green’s function retrieval . . . . . . . . . . . . . . . . . . . . . . . . . 744.3 Scattered Waves in the Green’s Function . . . . . . . . . . . . . . . . 754.4 Passive image interferometry . . . . . . . . . . . . . . . . . . . . . . . 754.5 Application to Merapi volcano . . . . . . . . . . . . . . . . . . . . . . 784.6 Hydrological modeling . . . . . . . . . . . . . . . . . . . . . . . . . . 794.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

5 Daily velocity variations and instrumental errors 855.1 Improving passive image interferometry . . . . . . . . . . . . . . . . . 87

5.1.1 Hourly Green’s functions . . . . . . . . . . . . . . . . . . . . . 875.1.2 Hourly velocity variations . . . . . . . . . . . . . . . . . . . . 895.1.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

5.2 Detecting instrumental errors with passive imaging . . . . . . . . . . 945.2.1 Cross-volcano Green’s functions . . . . . . . . . . . . . . . . . 945.2.2 Instrumental errors . . . . . . . . . . . . . . . . . . . . . . . . 955.2.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

Summary 101

Bibliography 103

Deutsche Zusammenfassung 113

vii

Acknowledgments

This thesis kept me busy during the last three years. This time would have beenmuch less fun without the support that I received from many people.

I am grateful to my supervisor Prof. Dr. Michael Korn. He provided guidancewhenever necessary and freedom when I desired.

I am also grateful to Prof. Dr. Michel Campillo for his advise and the oppor-tunity to work in his amazing group at the Laboratoire de Geophysique Interne etTectonophysique (LGIT) Grenoble. The LGIT provided computing facilities for theM-C simulations.

I like to thank Dr. Ulrich Wegler for introducing me to the topic of scatteringand for inspiring discussions.

Special thanks go to Dr. Ludovic Margerin for enlightening discussions, andthe idea to work on the Lg-blockage. I really appreciated his help in preparing mystay in Grenoble and his patience to help me with all the practical problems thatsuch a stay brings with it.

I am thankful for the numerous discussions with Dr. Jens Przybilla, SigwardFunke, Dr. Florent Brenguier, Pierre Gouedard, Laurent Stehly, and Dr. Eric Larosewhich were my colleagues in Leipzig and Grenoble. The experience of Jens Przybillawith Monte-Carlo simulations was very helpful. I thank Dr. Malte Westerhaus andBirger G. Luhr for discussion on Merapi and for meteorological data from Merapi.

I like to thank my wife and my kids for accompanying me during my work inGrenoble. My kids were always a good reason to keep me from work whenever itwas going to be too much. I thank my wife for a lot of good ideas that I had whileshe was carefully listening to my scientific problems.

I thank all the Institutes that provided data for this work. The German Seismo-logical Central Observatory (SZGRF) provided the seismograms analyzed in chap-ter 1. The French Atomic Energy Commission (CEA/DASE) granted access to theirgreat earthquake database. This data is used in chapter 3. The GeoForschungsZen-trum Potsdam (GFZ) acquired the seismic data from Merapi volcano that is usedin chapters 4 and 5. Special thanks for this interesting data go to Dr. Jo Wasser-mann. The US National Aeronautics and Space Administration (NASA) providedthe precipitation data used in chapter 4 and the relief data in figure 3.9.

I also appreciate financial support from the German Research Foundation(DFG, grant WE 2713/1-1), the German Academic Exchange Service (DAAD,Doktorandenstipendium), European Union through the NERIES program, and theLGIT Grenoble.

Overview

Beginning with the first recording of a teleseismic earthquake by Ernst von Rebeur-Paschwitz in 1889, seismology has lead to a variety of important discoveries. Suc-cessful methods have been developed to use seismic waves for various goals. Mostimportant among these is of course the seismic imaging, the extraction of infor-mation about the spatial distribution of structural properties. Another successfulbranch of classical seismology aims at the characterization of seismic sources.

These methods compare the observed wavefield with predictions made by adeterministic model. Since small scale structure is usually not contained in thesemodels the related occurrence of scattered waves is unwanted and sophisticatedtechniques are used to minimize their influence.

In this work a different approach is used. Instead of trying to eliminate scatteredwaves in seismic data they are utilized as source of information. The author presentsdifferent studies with different aims and results, but all are based on the multiplescattering nature of seismic waves. This work is divided into two parts according tothe different, and in a way complementary, theories that are applied.

Part I is based on Radiative Transfer Theory (RTT) and analyses the amplitudeof scattered waves in the seismic coda. This theory describes the propagation ofseismic energy in a scattering medium. It completely neglects phase informationof scattered waves. The first study presented in chapter 1 is published in Sens-Schonfelder and Wegler (2006b). It applies RTT in a simple approximation thatallows for an analytic solution of the Radiative Transfer Equation (RTE). Regionalseismogram envelopes of 11 earthquakes in Germany are compared with predictionsof RTT in a half space model with multiple isotropic scattering of acoustic waves.An inversion strategy is presented that uses amplitude information from ballisticshear waves and the seismic coda to simultaneously estimate the seismic moment ofthe source and site amplification factors of the seismic stations as well as attenuationand scattering properties of the crust. With a transport mean free path of about690 km and an intrinsic quality factor of 500 the distance dependent shape of thecoda envelopes can be well described. The seismic moments estimated with thesemedium parameters agree well with independent moment estimates using waveformmodeling techniques. A clear correspondence is found between the estimated siteamplification factors and geological conditions at the seismograph stations.

Despite the ability of the method presented in chapter 1 to estimate parametersof medium, source, and recording site, it relies on limiting approximations of theRTE. Seismogram envelopes can only be modeled in the S-wave coda. P-waves,mantle phases, and Moho reflections are not contained in this simple model. In

1

2

chapter 2 a tool is introduced that was developed to solve the RTE in a morerealistic model with less restrictive approximation. This algorithm uses a Monte-Carlo technique to numerically solve the elastic RTE with anisotropic scattering andconversions between compressional and shear waves. It allows to model seismogramenvelopes including ballistic P- and S-waves as well as the coda of both phases. To beapplicable for the modeling of energy propagation in the earth’s crust the algorithmcontains modules that treat the effects of the free surface, a velocity contrast atdepth and the resulting wave-guide effects. With the additional possibility to modelthe energy propagation in a model with vertical velocity gradient the, algorithm isable to simulate complete regional seismogram envelopes from the earliest arrival ofmantle phases until the late coda.

The ability of this algorithm to model the wave-guide effects of the crust isused in chapter 3 to study the phenomenon of Lg-blockage in the western Pyrenees.The term Lg-blockage refers to the strong attenuation of guided shear energy in alimited region of the crust. In contrast to other regions, the deterministic velocitystructure of the western Pyrenees cannot account for the Lg-blockage. Thus, thephenomenon remained unexplained. It was suggested that it is caused by scatteringat small scale heterogeneities. Chapter 3 starts with a detailed data analysis used toextract characteristic traces for propagation through the eastern and western partsof the mountain range. The parameters of the undisturbed crust are estimated bynonlinear inversion of the trace for propagation through the eastern Pyrenees. Forthe transport mean free path and the quality factor of intrinsic S-wave attenuationvalues of 750 km and 623 are found respectively. To simulate the Lg-blockage anadditional body with different scattering and attenuation properties is placed inthe crust below the western Pyrenees. The parameters of the medium inside thisbody are obtained by inversion of the trace for propagation through the westernPyrenees. The simulations support the hypothesis that the intrinsic properties ofthis body are responsible for the Lg-blockage. Increased scattering is necessary toattenuate the energy that propagates in the crust. Together with stronger intrinsicattenuation inside the Pyrenean body, the model correctly describes the observationsfor propagation through both the western and the eastern Pyrenees. It provides anexplanation for the Lg-blockage.

In part II a complementary theory to radiative transfer is applied that uses thephase information of scattered waves in the seismic coda. Contrary to part I the codasignals are not generated in the classical sense with an active source and a receiver.Instead the principle of recovering Green’s functions by correlation of random wavefields is applied here. This concept relays on the persistence of interference effects indiffuse wave fields. In chapter 4 a technique named Passive Image Interferometry isintroduced to use the continously present seismic noise to monitor velocity variationsof seismic waves. It is demonstrated that it is possible to retrieve scattered waves, i.e.seismic coda by correlation of seismic noise. With these scattered waves, spatiallyheterogenous velocity variations are observed with a temporal resolution of one day.A hydrological model is presented that explains the lapse time dependent velocityvariations on the basis of precipitation.

In chapter 5 two studies are presented that demonstrate the use of Green’s

3

function retrieval from seismic noise for different purposes. It is shown that thetemporal resolution of Passive Image Interferometry can be increased to observedaily changes. Evidence is presented for a connection between these changes andair temperature. In the other study a method is presented to identify and quantifyinstrumental errors in a seismic network. It is shown how the timing errors ofindividual instruments can be measured. This allows for the synchronization ofdata from stations with unreliable timing to the clocks of more precise stations.The accuracy of this synchronization can already exceed the precision of the internalinstrument clock after one day.

Part I

Radiative Transfer Theory

5

Chapter 1

Acoustic radiative transfer theoryfor the joint inversionof seismogram envelopes forsource- site- and path-effects

7

8 CHAPTER 1. ACOUSTIC RADIATIVE TRANSFER THEORY

Abstract

In the present study a new technique is proposed to obtain medium parameters, siteamplification factors, and source spectra of regional earthquakes from envelopes ofseismic coda. As compared to existing methods, this approach is based on radiativetransfer theory which physically describes the scattering process that produces theseismic coda. This allows the direct estimation of source parameters, without thenecessity to fix proportionality coefficients with reference events. I see an appreciableadvantage because the method is independent of the output from other techniques,such as reference events provided by moment inversions. The main component ofthis method is a joint inversion of the seismic records for source and site parameters,as well as for medium parameters assuming isotropic sources and isotropic, acousticscattering in a half-space.

The method is tested with recordings of 11 earthquakes (4 ≤ Ml ≤ 6) bythe German Regional Seismic Network at epicentral distances less than 1000 km.Traces are inverted in eight frequency bands between 0.2 and 24 Hz and it is shownthat estimates of the seismic moment are in good agreement with values obtainedin independent studies using waveform inversion techniques. In fact the estimatesof the seismic moment are better than approximations obtained from local magni-tudes using empirical relations specifically derived for the region under study. Theparameters that describe the scattering medium are mean free path that we foundto average around 690 km and the intrinsic quality factor for which IQ = 500 isobtained below 3 Hz.

1.1. INTRODUCTION 9

1.1 Introduction

Since it has been recognized for some time that amplitudes of regional coda areproportional to source excitation, attempts to extract information about the seis-mic source from the coda have a long history. Aki and Chouet (1975) were thefirst to study source spectra with coda waves. Due to the available instrumentsthey obtained source spectra above 3 Hz. A correction for attenuation is necessarythat involves an empirical quality factor describing both intrinsic and scatteringattenuation. The measurements are made relative to each other and have to beadjusted to reference events. Mayeda and Walter (1996) used two-dimensional mul-tiple scattering to approximate coda envelopes and to measure coda amplitudes.Above 0.2 Hz, additional empirical distance corrections were introduced. Mayedaand Walter (1996) proposed that the influence of body wave scattering, which isnot described by their 2D scattering model, required this empirical modification.They transformed the frequency resolved coda amplitudes to moment-rate spectrausing reference events. Mayeda et al. (2003) extended this empirical approach andended up with 12 free parameters describing the coda envelopes. Mayeda et al.(2003) tested the method and applied it to the Dead Sea Rift. Morasca et al. (2005)used it to analyze energy-moment scaling in the western Alps. The approach ofMayeda et al. (2003) and Morasca et al. (2005) is completely empirical and has noconnection to the physics of the scattering process. Dewberry and Crosson (1995)used the single-scattering model in a detailed analysis of seismic source spectra of 78north western U.S. earthquakes. They also provided coda Q and site amplificationestimates.

All these studies share the fact that they need reference events with knownsource spectra to fix proportionality coefficients to obtain absolute values for thesource spectra. In contrast, the approach presented here is independent of externalinformation because it uses radiative transfer theory which originates from a physicalmodel of the scattering and provides a direct relation between the amplitude ofthe coda and the source excitation, without any proportionality coefficient. To beable to solve the radiative transfer equation analyticly it is used in the acousticapproximation with isotropic scattering. In this approximation, the direct relationbetween source excitation and the coda level was previously used by Nakahara et al.(1998). Based on theoretical developments by Sato et al. (1997), they presented anapproach to study the source process of large earthquakes in great detail. Nakaharaet al. (1998) used a model of multiple isotropic acoustic scattering to invert for thespatial distribution of non-isotropic high frequency energy radiation on the faultplane. In this method, medium parameters are taken from studies that belong toanother branch of coda investigations.

The aforementioned studies focus on the properties of the propagation mediumrather than on the source and aim to separate the effects of intrinsic and scatteringattenuation to characterize the small scale heterogeneities and the dissipation of en-ergy. There are, for example, several studies (Fehler et al., 1992; Mayeda et al., 1992;Bianco et al., 2002) applying the Multiple Lapse Time Window Analysis (MLTWA)as developed by Hoshiba et al. (1991). MLTWA is based on multiple isotropic


acoustic scattering and uses ratios of energy integrated in three consecutive timewindows to separate intrinsic and scattering attenuation. An improvement to thegeometrical setting is due to Margerin et al. (1998), who model energy propagationin a layer above a half-space. Lacombe et al. (2003) use a similar model consistingof a scattering layer overlaying a transparent half-space to characterize the atten-uation properties of Lg-waves. All these studies apply coda-normalization (Aki,1980) to correct for unwanted effects of the source and site amplifications. Coda-normalization, however, fails for small events with a shorter coda, because the codacan be dominated by random seismic noise before the requirement of homogeneousdistribution of energy in space is fulfilled. In the present study, this disadvantageis overcome by including source excitation and site amplification directly into theinversion process.

I present an approach to estimate the source spectrum by jointly inverting theseismic record for medium parameters and source/site effects. The merit of thisintegrated approach is that

1. it is independent of external information such as the source spectrum of ref-erence events because a physical model of the scattering process is utilized forwhich the parameters are obtained in the inversion. This model provides adirect relation between coda amplitude and source excitation.

2. it is applicable for smaller events because the source excitation and the siteamplification factors are estimated within the inversion, thus no coda-norma-lization is needed.

This content of this chapter is published in Sens-Schonfelder and Wegler (2006b).It is organized as follows. I describe the modeling of seismogram envelopes in sec-tion 1.2 and introduce the data selection and the study area in section 1.3. Theinversion scheme is described in section 1.4 and results of the inversion are shownin section 1.5.

1.2 Envelope modeling

Envelopes are modeled using radiative transfer theory. Scattering of waves in 3Dspace is governed by the radiative transfer or Boltzmann equation (Margerin et al.,1998). Here I restrict myself to isotropic scattering of S-waves in a half-space withan isotropic source. In full-space the effective energy density Green’s function G(t, r)is given by the following integral equation

G(t, r) = F (t, r) + v0g0

∞∫

−∞

∫

V

F (t− t′, r− r′)G(t′, r′)dt′dr′ (1.1)

(Sato and Fehler, 1998, p. 175).

F (t, r) =1

4πv0r2H(t)δ(t− r/v0)e−v0got (1.2)

1.2. ENVELOPE MODELING 11

is the Green’s function for the coherent wave energy. Here v0 and g0 denote averagevelocity of S-waves and total scattering coefficient respectively. r = |r| is the sourcereceiver distance and H is the Heaviside step function. The exact solution of eq. 1.1involves a two dimensional Fourier transform (Zeng et al., 1991; Sato and Fehler,1998, p. 177).

To speed up the inversion an analytic approximation to the solution of 1.1 isused that was obtained by Paasschens (1997) by interpolating between the explicitsolutions of the Boltzmann equation in 2D and 4D. The solution of Paasschens(1997) reads

G(t, r) ' e−v0tg0

4πr2δ(r − v0t) +

(1− r2/(v20t

2))1/8

(4πv0t/(3g0))3/2

×e−v0tg0K

(v0tg0

[1− r2

v20t

2

]3/4)H(v0t− r) (1.3)

where K(x) ' ex√

1 + 2.026/x. Recently this interpolation was used by Abubakirov(2005) in a MLTWA. It is a good approximation to the solution of eq. 1.1 given byZeng et al. (1991). The deviation from the exact solution is below 5% where largerdeviations occur in the tail of the direct wave. The accuracy of the approximationfor the direct wave and its tail is reasonable in the present context because only anaverage value in a short time window following the direct wave will be used. Thistime window is largely dominated by the energy of the unscattered wave which iscorrectly described in equation 1.3.

The boundary condition for radiative transfer in a half-space is zero verticalenergy flux at the surface. Assuming total reflection, this condition can be accountedfor by introducing a mirror source above the surface for which the energy densityGreen’s function is G(t, r+), with r+ being the distance between receiver and mirrorsource. At the surface the upward flux of the real source equals the downwardflux of the mirror source, thereby satisfying the boundary condition. The energydensity Green’s function of the half-space is thus Gh(t,x) = G(t, r) + G(t, r+). Inthe following it is assumed that the receiver is at the surface where r = r+. In thiscase Gh(t,x) = 2G(t, r) (c.f. Wegler, 2004). With this treatment only the boundarycondition for energy transfer is considered. Specific effects arising from boundaryconditions of the wave equation like mode conversions, surface waves and angledependent reflection coefficients are not accounted for.

The energy density for an arbitrary source can be obtained by convolutionwith the source function. In this analysis I assume a source function of the formWδ(r)δ(t) where W is the spectral source energy in J/Hz. Intrinsic absorption canbe accounted for with an additional time dependent factor e−bt with the intrinsicabsorption parameter b. Finally I obtain

Emod(t,x) = WRiGh(t,x)e−bt (1.4)

= 2WRiG(t, r)e−bt (1.5)

for the energy density Emod(t, r) of the half-space model at time t and distance rfrom the source of energy W . Ri is the site response at station i.


13-4-1992

20-6-1995

21-10-1997

29-11-1997

23-6-200122-7-2002

22-2-2003 22-3-2003

22-5-2004

20-10-2004

5-12-2004

BFO

BRNL

BUG

CLZ

FUR

GRB1

GRC1

HAM

TNS

WET

GRFO/GRA1

BRG

CLL

MOX

STU

GSH

GEC2

IBBN

BSEG

RUE

HLG

RGN

NRDL

UBBA

4˚

4˚

8˚

8˚

12˚

12˚

16˚

16˚

48˚ 48˚

50˚ 50˚

52˚ 52˚

54˚ 54˚

Figure 1.1: Map of Germany and neighboring countries. Gray triangles: broadbandstations of the German Regional Seismic Network, the Grafenberg array and the GERESSstation GEC2 (station codes are added); Black stars: earthquakes used in this study (datesof occurrence are added).

1.3 Data and regional setting

Data for this study was recorded by the German Regional Seismic Network (GRSN),the Grafenberg array (GRF) and the reference station of the GERESS-array. Thestations are equipped with STS-2 or STS-1 (GRF stations) broadband seismometersand traces are sampled at 80 Hz. Refer to Korn (2002) for a detailed discussion fromthe GRSN. For this work data of 25 different stations is used (figure 1.1).

As sources I use local and regional events from Germany and adjacent areas inthe period from the installation of the GRSN in 1991 until December 2004. Eventswith local magnitudes estimated by the German Central Seismological Observatory(SZGRF) larger or equal to 4 are chosen. To reduce scatter in the estimates of theparameters events from the Alps are excluded because of differences in the geologicalconditions. One event that occurred on 18 October 1994 in the North Sea is excluded

1.4. INVERSION 13

because I observed suspicious differences in the arrival times, perhaps due to locationerrors or an unsuitable velocity model. In the end 11 earthquakes (figure 1.1) areused and more than 100 3-component broadband records. Locations, magnitudesand depths as provided by the SZGRF are listed in table 1.1. With the possibleexception of the Roermond 1992 earthquake, where depths even below 20 km havebeen estimated (Korn, 2002), the events occurred in the upper crust. The last threecolumns of table 1.1 contain seismic moments estimated by waveform modeling andthe method presented here. They will be discussed later.

Preparation of data involves filtering of the seismograms in eight frequencybands centered at 0.1875, 0.375, 0.75, 1.5, 3.0, 6.0, 12.0 and 24.0 Hz with a narrownormalized Gaussian filter. Choosing the filter such that

∫∞−∞ |B(f)|2df = 1 where

B(f) is the frequency response, ensures that the spectral energy density is preservedin the filter process (Wegler et al., 2006a). Envelope sections are selected accordingto the following criteria:

1. S/N - ratio greater than 4

2. No obvious disturbances e.g. aftershocks (checked manually)

3. Hypocentral distance greater than 60 km.

The hypocentral distance range is limited to distances greater than 60 km because ofthe simple half space velocity model, with a mean shear wave speed v0 = 3.5 km/s,which is inappropriate for shorter distances.

The observed energy density Eobs(t, r), i.e. the seismogram envelopes, is ob-tained from the bandpassed velocity seismogram u(t, r) using

Eobs(t, r) =ρ0

2

u2(t, r) +H2(u(t, r))

2. (1.6)

Here ρ0 is the mean density of the medium and H denotes the Hilbert transform.ρ0 = 2700 kg/m3 is used throughout this study.

1.4 Inversion

In the inversion scheme I estimate values for the parameters g0, b,W and Ri thatminimize the misfit function

ε =

N∑

i=1

endi∑

j=starti

[log

(Eobs(tj, ri)

Emod(tj, ri)

)]2

. (1.7)

Here N denotes the number of stations. starti and endi correspond to the indicesof the first and the last sample in the “coda” time window (figure 1.2) at the i− thstation. The “coda” time window starts after twice the S-wave travel time and endswhen the S/N drops below 4. All available traces of a certain event in one frequencyband are simultaneously inverted for g0, b,W and Ri. The inversion scheme is acombination of a one dimensional grid search for g0 and a linear inversion for theremaining parameters. Figure 1.2 illustrates the constraints of this process whichworks in the following five steps.


0 100 200 300 40010−15

10−10

10−5

directS−wave seismic coda

lapse time [s]

E [J

m−3

Hz−1

]

envelopemodelinversiontime windows

Figure 1.2: Example of an observed seismogram envelope (thin curve) on a logarithmicscale with the best fitting model (bold black curve). Time windows used in the inversionare indicated at the bottom of the figure. Note that in the time window of the directS-wave only the average values (black dot) of envelope and model are fitted.

1. A value of g0 is picked (a range of g0 is probed in a grid search) and theenergy density Green’s functions G(t, ri) is calculated for all stations i, usingthe approximation in eq. 1.3.

2. I equate Eobs (eq. 1.6) and Emod (eq. 1.5) and rearrange terms like

ln

(Eobs(tj, ri)

G(tj, ri)

)= ln (2WRi)− btj . (1.8)

Here ln denotes the base e logarithm. By fitting a linear curve to the lefthand side of eq. 1.8 as function of time the appropriate values of b and (WRi)are obtained. To reduce the effect of direct surface waves and anisotropicscattering the graph is only fitted in the ”coda” time window (figure 1.2)starting at twice the S-wave travel time. Additionally I apply the constraintthat average values of the modeled and observed envelopes within the ”directS-wave” time window (figure 1.2) are equal. At distance r from the source thistime window contains the samples with lapse time t in the range 3.5 km/s >r/t > 3.0 km/s. This constraint is applied to resolve the tradeoff between g0

and (WRi) that exists in the later coda.

3. With the values of b and (WRi) obtained in step two Emod is calculated ex-plicitly (eq. 1.5) to estimate the misfit ε (eq. 1.7).

4. Steps one to three are repeated for a range of g0 values. The final model forone event and one frequency band is that with the smallest value of ε.

5. For the best fitting model W and Ri are calculated explicitly by simply takingW to be the logarithmic average of all (WRi) values. This way it is implicitly

1.5. RESULTS 15

assumed that the average of the logarithms of the Ri is zero. This meansthat site amplifications are measured relative to the network mean which isassumed to be one.

The constraint that makes use of the ”direct S-wave” time window originatesfrom the following problem. In a weakly scattering medium the difference betweenthe coda of two models with different W and g0 can be very small provided theproduct Wg0 is equal for both models. The single scattering approximation (Satoand Fehler, 1998, p. 47) which is valid for small g0 predicts that the energy densityEmod is proportional to Wg0. In this range it is inherently impossible to separateW and g0. But also for moderate scattering that is not adequately modeled in thesingle scattering approximation it can be practically impossible to separate W andg0 due to noisy data and other simplifications in the model if only information fromthe coda is used.

In contrast the ballistic energy density is proportional to We−(v0g0+b)t. Herea tradeoff is observed between g0 and b which corresponds to the known fact thatintrinsic and scattering attenuation can not be separated using unscattered energyonly. However, the energy W and total attenuation can be separated because of thedifferent time dependence.

In the inversion the envelopes are fitted in the “coda” time window indicatedin figure 1.2. But I restrict the possible models to these which correctly predict theballistic energy in the “direct S-wave” time window (figure 1.2). I want to stressthat in a weakly scattering medium it is on the one hand not possible to reliablyestimate W from the coda alone without a priory knowledge of g0. On the otherhand it is impossible to estimate W solely from the ballistic wave because in orderto apply radiative transfer theory which provides the source energy W I need toseparate g0 and b. I note that information from the unscattered wave is implicitlyevaluated also by MLTWA to separate g0 and b since the first time window of theMLTWA generally contains the ballistic energy.

1.5 Results

As all available recordings are inverted simultaneously I implicitly assume that thereare no lateral variations in the structural parameters g0 and b. The results representintegral values averaged over possible regional variations. An example of the fit thatis achieved with the model is given in figure 1.3. The observed envelopes are plottedas thin lines whereas bold lines represents the model envelopes. Average energydensity values in the “direct S-wave” time window are indicated as gray and blackdots for observed and modeled envelopes respectively. The model generally fits theenvelopes in the “coda” time window. Deviations in the “direct S-wave” windoware probably due to directionality of the source.

In the following I will first present the results of the attenuation parametersas they are fundamental for the estimated source parameters that are presentedthereafter. I will also briefly mention the estimates of the site amplifications.


10−16

10−14

10−12

10−10

10−8

BFOr=123km

0 200 400 600

BRGr=594km

lapse time [s]0 200 400 600

BSEGr=671km

10−16

10−14

10−12

10−10

10−8

BUGr=346km

10−16

10−14

10−12

10−10

10−8

CLZr=469km

FURr=342km

0 200 400 600

GEC2r=520km

GRA1r=364km

10−16

10−14

10−12

10−10

10−8

GRB1r=382km

ener

gy d

ensi

ty E

[J m

−3 H

z−1]

GRC1r=363km

GRFOr=364km

IBBNr=446km

MOXr=439km

0 200 400 600

10−16

10−14

10−12

10−10

10−8

RUEr=681km

STUr=191km

TNSr=245km

WETr=464km

Figure 1.3: Comparison between the envelopes (thin curves) of the 22-2-2003 earthquakewith the model (bold black curve) in the 1.5 Hz band on a logarithmic scale. Time windowsused in the inversion are indicated with gray bars at the bottom of the plots. Note thatin the time window of the direct S-wave only the average values of envelope and modelare fitted. These levels are indicated as gray and black dots for envelope and modelrespectively.

1.5. RESULTS 17

1.5.1 Attenuation parameters

Scattering strength can be expressed in terms of the total scattering coefficient g0

which is the inverse of the mean free path or in terms of scattering attenuationparameter ScQ−1. The relation between both is

ScQ−1 =g0v0

2πf(1.9)

where f denotes frequency. Figure 1.4 displays g0 as function of frequency. Smallgray dots in figure 1.4 indicate measurements of individual earthquakes. The scat-ter of these values is due to uncertainties of the measurements but additionallyregional differences of the scattering strength will increase the variations in the g0

values. Black dots with error bars indicate the logarithmic averages of the individ-ual measurements together with their 95% confidence limits. Individual measure-ments range between 10−7 m−1 and 10−5 m−1. Logarithmic averages range between9·10−7 m−1 and 3·10−6 m−1. There is no significant frequency dependence of g0 andlogarithmically averaged over all measurements irrespective of frequency a value of1.45 · 10−6 m−1 is obtained corresponding to a mean free path of 690 km. I comparemy estimates with Abubakirov and Gusev (1990); Fehler et al. (1992); Mayeda et al.(1992) and Lacombe et al. (2003) who studied coda waves in Kamchatka, Japan,Hawaii/Long Valley/central California and France respectively. Table 1.2 summa-rizes the results of the different studies in frequency bands centered at 3.0 Hz. Referto Bianco et al. (2002) for graphical representation of most of these results. Thevalue obtained in the present study is low in this comparison but the difference canbe attributed to geological distinctions since the central European intraplate regionis likely less heterogeneous than the volcanic areas listed in table 1.2. Comparedto the results from France obtained by Lacombe et al. (2003) the difference mightpartially be due to a different model setup.

The intrinsic attenuation parameter IQ−1 is related to the absorption parameterb as

IQ−1 =b

2πf. (1.10)

Figure 1.5 displays the results of individual measurements with small gray dots andlogarithmic averages in the individual frequency bands together with 95% confidencelimits as black dots with error bars. Averaged over individual measurements inthe separate frequency bands values of IQ−1 between 2.6 · 10−4 and 2.2 · 10−3 areobtained. For frequencies below 3 Hz IQ−1 is approximately constant at 2 · 10−3.Above approximately 3 Hz a power law decrease like IQ−1 ∝ f−1 is observed. Thesevalues are within the range of results from different regions listed in table 1.2 andvery close to the results from France Lacombe et al. (2003). However, there maybe some bias because regional records are inverted with a half space model. Iprobably overestimate IQ−1 because leakage of energy through the Moho into aweakly scattering mantle is mapped into attenuation (Margerin et al., 1998; Wegler,2004).


10−1 100 10110−8

10−7

10−6

10−5

10−4

f [Hz]

g 0 [m−1

]

individual eventsmean with 95% confidence limits

Figure 1.4: Results for the total scattering coefficient of S-waves. Gray dots: scatter-ing coefficients measured by inversion of individual events; Black dots with error bars:logarithmic mean of individual measurements with 95% confidence limits.

1.5.2 Source energy

Since the present approach is on one hand based on a physical model of energypropagation and on the other hand does not require a coda normalization, theinversion parameter W can be interpreted directly. By measuring the source energyat various frequencies I obtain the source energy spectrum W (ω) of radiated S-waves. Assuming particle motion is caused by a double couple and is observed inthe far field one can obtain a relation between W (ω) and the seismic moment byintegrating the energy flux density over a sphere containing the source. Sato andFehler (1998, p. 152) state

W (ω) =ω4 |M(ω)|2

10πρ0v50

(1.11)

where v0, ρ0 and M(ω) denote mean S-wave velocity, mean density, and the Fouriertransform of the moment time function respectively. From eq. 1.11 one obtains thesource displacement spectrum ω|M(ω)|.

In figure 1.6 the source displacement spectra of the events listed in table 1.1are shown. Error bars correspond to the resolution of the inversion process andnot to the scatter of individual measurements as in figure 1.4 and 1.5. The misfitε in eq. 1.7 is the variance of the logarithmic difference between the samples ofthe observed and modeled energy densities. I use the F-test to decide whether twomodels are significantly different based on the ratio of their variances ε. The errorbars mark the range of values for the displacement spectra that can be obtained formodels which cannot be distinguished from the best model according to the F-test(Buttkus, 2000, p. 231) with a 5% significance level.

Most of the curves show the expected characteristics of a displacement spec-trum. A flat region can be observed towards the low frequency limit whereas the

1.5. RESULTS 19

10−1 100 10110−5

10−4

10−3

10−2

f [Hz]

I Q−1

individual eventsmean with 95% confidence limits

Figure 1.5: Results for the attenuation parameter of intrinsic absorption. Gray dots:attenuation parameter measured by inversion of individual events; Black dots with errorbars: logarithmic mean of individual measurements with 95% confidence limits.

displacement decays above a corner frequency. A different behavior is found for theOctober 2004 event that occurred in Northern Germany which is tectonically quiet.Neither low frequency plateau nor corner frequency can be observed. I regard thisas an evidence for a complex source process.

The low frequency plateaus can be compared with seismic moments Mwm0 in-

dependently obtained from waveform inversions. For most of the events the so-lutions of the Swiss regional moment tensor catalog that is available on-line athttp://www.seismo.ethz.ch/mt/ are used. Refer to Braunmiller et al. (2002) for ex-amples and a description of the method. Mwm

0 estimates for the three events before2000 can be found in Braunmiller et al. (1994) and Braunmiller (2002). For the eventin 1995 there is no seismic moment available. The waveform inversion estimates ofM0 are plotted as thick lines in figure 1.6. Obviously the low frequency values of thedisplacement spectra obtained in this study match the moments from the waveforminversion well. Assuming that the corner frequencies of the events used here exceed1.5 Hz I take the logarithmic average of the displacement spectra in the four lowestfrequency bands as a coda based estimate of the seismic moment denoted M coda

0 . Infigure 1.7 the values of seismic moment obtained in this study (M coda

0 ) are plottedagainst the values from waveform modeling (Mwm

0 ). On average there is a differenceof 36% between the waveform and the coda results.

Another reference that is interesting to compare with the displacement spec-tra is the local magnitude scale. Recently Reamer and Hinzen (2004) investigatedearthquakes in the Northern Rhine Area (around station BUG in figure 1.1) andestablished the relation

log10 MMl0 = 1.083Ml + 10.215 (1.12)

between the local magnitude Ml and the seismic moment MMl0 in Nm. The super-


1011

1013

1015

1017

1992_04_13 1995_06_20 1997_10_21 1997_11_29

1011

1013

1015

1017

2001_06_23

sour

ce d

ispl

acem

ent s

pect

rum

ω|M

(ω)|

[Nm

]

2002_07_22 2003_02_22

10−1 100 101

2003_03_22

10−1 100 1011011

1013

1015

1017

2004_05_22

10−1 100 101

2004_10_20

f [Hz]

10−1 100 101

2004_12_05

Figure 1.6: Double logarithmic plots of source displacement spectra of the 11 earth-quakes used in this study. Numbers within the charts indicate the dates of the eventsin yyyy mm dd format. Dots with continuous curve: source spectra estimated in thisstudy. Error bars correspond to 95% confidence limits of the resolution of the inversion(see text for explanation). The low frequency part shown in black is averaged to obtainM coda

0 . Thick horizontal line: seismic moment Mwm0 estimated independently using wave-

form modeling techniques (Braunmiller et al., 2002). Thin dashed curve: omega - squarespectra independently calculated from local magnitudes using empirical relations. Lowfrequency limit corresponds to MMl

0 .

script Ml indicates that the seismic moment MMl0 is estimated from local magni-

tudes. Figure 1.7 also shows the comparison of MMl0 with the results from waveform

modeling (Mwm0 ). On average the difference between the magnitude-derived and the

waveform-derived moments is 57%. The scatter in the magnitude derived estimatesis thus larger than that of the coda derived estimates. However, according to theF-test with 10 degrees of freedom the difference between this 57% and the 36% meandeviation of the coda-derived moment is not significant at a 95% confidence level.This means the coda derived estimates of the seismic moment are at least as goodas the estimates that can be obtained directly from the local magnitudes.

Assuming an omega-square source model according to ω|MMl(ω)| = MMl0 /(1+

f/fc)2 (Aki and Richards, 1980a, p. 823) one can use equation 1.12 to calculate

an approximate source displacement spectrum if the corner frequency fc is known.

1.5. RESULTS 21

1014 1016 10181014

1015

1016

1017

1018

M0 fr

om c

oda

/ mag

nitu

de [N

m]

M0wm from waveform inversion [Nm]

M0coda

M0Ml

Figure 1.7: Scatter plot of the seismic moments obtained from seismic coda (M coda0 )

in this study and moments estimated from local magnitudes (MMl0 ) against estimates of

seismic moments from long period waveform modeling (Mwm0 ). Straight line indicates the

locations of perfect agreement of the different estimates.

Further assuming a stress drop σ of 10 MPa, it can be obtained from fc = 3√σ/MMl

0 ·0.49v0 (Hough et al., 2000). These approximate source spectra are plotted as thindashed lines in figure 1.6. Note that these approximate source spectra are solelybased on the local magnitude estimated by the SZGRF.

Figure 1.6 shows that also the local magnitudes can give reasonable estimatesof the seismic moment. There are some differences between the magnitude- andcoda-derived spectra in the corner frequency and in the high frequency asymptote.ω|M coda(ω)| shows a steeper decline than the omega-square model.

1.5.3 Site response

Site amplification factors for velocity amplitudes range between 0.3 and 3. Mostlythey can easily be related to local geology (table 1.3). In the 0.375 Hz band for ex-ample we find amplifications of crystalline sites only below 1 (logarithmic mean: 0.8)whereas stations on sediments show amplifications exclusively above 2 (logarithmicmean: 3.0). Sites on sedimentary rocks occupy an intermediate range between 0.75and 1.75 (logarithmic mean: 1.1). As examples the site factors for the stations BUGand FUR in figure 1.8 are shown. BUG is situated on a clastic sedimentary rock.The site amplification is not frequency dependent with a mean of 0.9. In contrastFUR is placed on moraine over molasse (i.e. sediments). Here a clear amplificationof a factor greater than 2 is seen. Response factors for all stations and frequenciesare listed in table 1.3, together with information about geology.


10−1

100

101

BUG

site

am

plifi

catio

n fa

ctor

R

10−1 100 10110−1

100

101

FUR

f [Hz]

site

am

plifi

catio

n fa

ctor

R

individual events

mean with 95%confidence limits

Figure 1.8: Site response spectra for the stations BUG and FUR. BUG is placed onsedimentary rock whereas FUR is situated on sediments.

1.6 Discussion

The model used in the present study is a first-order approximation of S-wave scatter-ing in the earth. It neglects anisotropic scattering and mode conversions, does notpay attention to source radiation patterns and assumes a statistically homogeneoushalf-space. Some of these effects could be incorporated into the model by solving theradiative transfer equation with Monte-Carlo techniques. This allows to take intoaccount anisotropic scattering or mode conversion of elastic waves (Przybilla et al.,2006). Even anisotropic sources and deterministic earth structure, like a depth de-pendent velocity and scattering coefficient can be simulated (Margerin et al., 1998,chapter 2). The improvement of the model that can be achieved by incorporatingthese effects will be variable and the successful application of more complex modelsrequires detailed knowledge about the medium and the source (Hoshiba et al., 2001).The source radiation patterns for example will probably have a minor effect on thepresent results because in the late coda the effect vanishes due to the averaging bymultiple scattering and the direct S-wave is measured and averaged over a lot of

1.6. DISCUSSION 23

stations. Note that S-waves always have less pronounced radiation patterns thanP-waves irrespective of the station configuration (Sato and Fehler, 1998, p.150).

The acoustic approximation is justified here because the different magnitude ofthe conversion scattering coefficients for S to P and P to S conversion (Aki, 1992)leads to a stable partitioning of energy at large lapse times which is dominated byS-waves. In this regime P-waves carry only about 10% of total energy (Ryzhik et al.,1996, eq. 5.40). In order to model the early S- and P-wave coda the treatment of theelastic scattering problem seems to be the most promising modification to the presentapproach (c.f. chapter 2). The computational effort, however, is substantially higher.A still unsolved problem is the treatment of conversion and scattering of surfacemodes.

One consequence of the half-space model is the inability to handle Lg-waves.Lacombe et al. (2003) model Lg-waves using a Monte-Carlo technique to separateintrinsic absorption and scattering attenuation. Interestingly Lacombe et al. (2003)could not satisfactorily separate elastic and anelastic attenuation without the useof information from direct waves incorporated as total attenuation. This directlycorresponds to my observation that led to the constraint of the inversion with theballistic energy.

Nevertheless, the present results are encouraging. Though a very simple modelis used I can directly relate the shape and amplitude of seismogram envelopes tothe moment of the source. These moment estimates are in good agreement withindependent results from moment tensor inversions. I find reasonable values for thetotal scattering coefficient g0 and the intrinsic attenuation parameter IQ−1 and ob-tain site response factors that reflect local geology. I think that the approximationsare justified in order to obtain fast estimates of relevant parameters without a prioriknowledge about medium or source.

There are some technical parameters in the algorithm that might influence theresults. First, the choice of the misfit function might have an effect. I verified thatlog(Eobs/Emod) is approximately Gaussian distributed in the coda time window ifthe model fits well. For comparison I inverted the dataset using different misfitfunctions including the L1-norm and found no significant differences.

Another influence might come from the choice of the inversion time windows.The “seismic coda” time window starts at two times the S-wave travel time. Inthis part of the seismogram the distribution of energy is traditionally believed to beclose to the diffusion limit (Rautian and Khalturin, 1978; Wu and Aki, 1985). Thisis not required in the present method but also the assumption of isotropic scatteringis a better approximation in the later coda. Additionally direct surface waves areexcluded by fitting only the late coda. The “direct S-wave” time window starts atthe arrival of the direct S-wave (v = 3.5 km/s) and lasts until waves with apparentvelocities above 3.0 km/s have passed the station. This is somewhat arbitrary butthe objective was to collect energy that was radiated as S-wave from the source inthe direction of the receiver. Taking into account finite source duration and envelopebroadening due to forward scattering (Sato, 1989; Saito et al., 2002) this appearedto be a reasonable choice. By comparing the results with inversions done withdifferent lengths of the ”direct S-wave” time window I verified that the influence is


small compared to the uncertainties of the estimates.

A new aspect of this method is the frequency range used for inversion. Typicallyonly frequencies above 1 Hz have been used in studies of seismic scattering (Fehleret al., 1992; Mayeda et al., 1992; Bianco et al., 2002; Lacombe et al., 2003). On onehand this is due to the short period data used in these studies. On the other handlow frequency coda is believed to be influenced by surface waves (Wu and Aki, 1985;Abubakirov and Gusev, 1990). The latter argument appears to be not well founded.Mode conversions between body and surface waves as well as surface to surfacewave scattering are not bound to frequencies below 1 Hz. Since it is importantfor the purpose of estimating the low frequency level of the source spectrum, i.e.the seismic moment of the earthquakes, the broadband data is further exploited toobtain information below the corner frequencies.

The major advantage of this method is that I can easily obtain source parame-ters like the seismic moment without a priori knowledge about attenuation functionsthat are needed for magnitude estimations and without waveform modeled referenceevents that are needed if envelopes are approximated on the basis of empirical for-mulae. The approach is thus optimally suited for source parameter estimation withtemporal seismic networks.

Another possible application is related to the stability of the coda compared todirect waves. If the medium parameters g0 and IQ−1 and the site response factorsRi have been estimated once in a study like the present, the seismic coda can beused to obtain magnitude estimates that are more stable than estimates from directwaves because the coda is less sensitive to source radiation patterns. In fact Mayedaet al. (2003) showed that coda amplitudes have 3−5 times less interstation scatter ascompared to direct phases. Just like the method of Mayeda et al. (2003) my approachcan be used to significantly improve the magnitude estimates in permanent seismicnetworks.

Compared to waveform modeling techniques the present method has a funda-mental advantage for the estimation of the seismic moment. Modeling waveformsrequires knowledge about the Green’s function of the propagation medium. Dueto small scale heterogeneities in realistic earth models this is impossible for highfrequencies. Therefore waveform techniques can not be applied to small (Ml < 3.5)earthquakes. As I model energy densities rather than waveforms the approach onlyrequires knowledge of the energy density Green’s function. In this case the phaserelation between waves that interacted with these small scale structures can be ne-glected and it is possible to obtain the Green’s function in a statistical sense. It isthus possible with this method to obtain estimates of the seismic moment for smallevents with high frequency sources.

1.7 Conclusions

A new method is presented to extract information from the seismic coda using ra-diative transfer theory in the acoustic approximation. It is based on the modeling ofcoda envelopes as S-waves that were isotropically scattered at randomly distributed

1.7. CONCLUSIONS 25

heterogenieties in a statistically homogeneous half-space. I do not apply coda nor-malization but instead incorporate the effects of different sources and site responsesinto the inversion scheme which allows to explicitly resolve these parameters togetherwith medium parameters.

Source spectra and seismic moments can thus directly be obtained from seis-mogram envelopes and it is shown that there is excellent agreement between myestimates and values independently obtained in moment tensor inversions. Addi-tionally, site response factors are obtained and intrinsic and scattering attenuationare separated.

26C

HA

PT

ER

1.A

CO

UST

ICR

AD

IAT

IVE

TR

AN

SF

ER

TH

EO

RY

Table 1.1: List of earthquakes used in this study. Source time, location, depth, and local magnitudes are provided by the GermanCentral Seismological Observatory. Seismic moments estimated by modeling waveforms are listed in column 6 with respective referencesin column 7. The last column lists seismic moments estimated in the present study.

Date Time Location Ml Depth Mwm0 in Nm Reference M coda

0 in Nm(waveform modeling) (this study)

1992/04/13 01:20:03.1 51.14oN 6.05oE 6.0 15 9.2·1016 Braunmiller et al. (1994) 9.3·1016

1995/06/20 01:54:57.5 50.57oN 4.82oE 4.6 10 9.0·1014

1997/10/21 16:44:39.4 48.78oN 9.72oE 4.0 10 2.3·1014 Braunmiller (2002) 1.7·1014

1997/11/29 20:06:09.3 50.31oN 8.35oE 4.0 7 3.0·1014 Braunmiller (2002) 4.5·1014

2001/06/23 01:40:05.3 50.87oN 6.17oE 4.3 10 9.23·1014 www.seismo.ethz.ch/mt 5.4·1014







1.7.C

ON

CL

USIO

NS

27

Table 1.2: Comparison of structural parameters g0 and IQ−1 for different regions in frequency bands centered at 3 Hz.

Reference g0 in m−1 IQ−1 Region Distance range ModelAbubakirov and Gusev (1990) 6.7·10−6 0.0032 Kamchatka local half space

Fehler et al. (1992) 6.5·10−6 0.0026 Kanto-Tokai (Japan) local half spaceMayeda et al. (1992) 4.26·10−5 0.00317 Long Valley local half space

2.82·10−5 0.00336 Central California local half space3.19·10−5 0.00055 Hawaii local half space

Bianco et al. (2002) 2.7 · 10−6 0.0034 Italy local-regional half spaceLacombe et al. (2003) 4·10−6 0.0012 France regional layer over half-space

Abubakirov (2005) 5.3·10−6 0.0031 Kamchatka regional half spacethis study 1.45·10−6 0.0013 Germany regional half space

28C

HA

PT

ER

1.A

CO

UST

ICR

AD

IAT

IVE

TR

AN

SF

ER

TH

EO

RY

Table 1.3: Site response factors of the 25 seismic stations that were used in this study. The “÷” sign denotes the logarithmic uncertainty.The site factor is to be multiplied or divided by the uncertainty to obtain the limits of the 95% confidence interval. If no uncertainty isgiven only a single event could be inverted. The two stations marked with * are placed several 100 m below surface.

Station Geology 0.1875 Hz 0.375 Hz 0.75 Hz 1.5 Hz 3.0 Hz 6.0 Hz 12.0 Hz 24.0 HzBFO granite 0.71 ÷ 1.61 0.76 ÷ 1.26 0.72 ÷ 1.16 0.73 ÷ 1.20 0.64 ÷ 1.28 0.57 ÷ 1.18 0.94 ÷ 1.54 1.06 ÷ 2.85

BRNL cenozoic sediments 3.78 4.61 3.72 ÷ 7.52 2.03BRG slate 1.10 ÷ 2.09 0.81 ÷ 1.46 0.56 ÷ 1.15 0.60 ÷ 1.16 0.90 ÷ 1.23 1.21 ÷ 8.32BSEG anhydrit 2.11 ÷ 1.19 1.23 ÷ 2.04 0.58 ÷ 1.52 0.62 ÷ 2.05 1.13 0.53BUG sedimentary rocks 0.87 ÷ 1.37 0.96 ÷ 1.41 0.88 ÷ 1.25 0.86 ÷ 1.21 1.08 ÷ 1.39 1.11 ÷ 1.42 0.72 ÷ 1.39 0.85 ÷ 1.31CLL Greywacke rocks 1.12 ÷ 2.18 0.91 ÷ 1.27 0.74 ÷ 1.17 1.10 ÷ 1.21 1.20 ÷ 1.23 1.85 ÷ 1.40CLZ sedimentary rocks 0.87 ÷ 1.20 0.86 ÷ 1.22 0.91 ÷ 1.13 0.84 ÷ 1.09 0.91 ÷ 1.12 1.37 ÷ 1.12 2.13 ÷ 1.37FUR moraine over Molasse 2.80 ÷ 1.37 3.45 ÷ 1.14 2.49 ÷ 1.11 1.89 ÷ 1.14 1.89 ÷ 1.34 1.06 ÷ 2.99GEC2 granite 1.11 ÷ 1.27 0.77 ÷ 1.17 0.65 ÷ 1.20 0.68 ÷ 1.13 0.75 ÷ 1.45 1.30 ÷ 1.57GRA1 limestone 1.28 ÷ 1.18 1.58 ÷ 1.24 1.69 ÷ 1.08 1.26 ÷ 1.09 0.98 ÷ 1.11 0.63 ÷ 1.24GRB1 limestone 1.04 ÷ 1.19 1.32 ÷ 1.17 1.61 ÷ 1.09 1.73 ÷ 1.11 2.97 ÷ 1.18 2.51 ÷ 1.23GRC1 limestone 0.97 ÷ 1.33 0.86 ÷ 1.22 0.97 ÷ 1.10 1.18 ÷ 1.16 0.79 ÷ 1.14 0.79 ÷ 1.50GRFO limestone 1.31 ÷ 1.17 1.57 ÷ 1.24 1.55 ÷ 1.08 1.00 ÷ 1.07 0.51 ÷ 1.06 0.63 ÷ 1.20GSH slate 0.28 0.97 0.71 0.97 1.43 1.77 1.00HAM cenozoic sediments 7.26 5.03 ÷ 43.31 2.52 ÷ 8.11 1.40 1.46HLG sandstone 1.87 ÷ 3.51 1.37 ÷ 92.54 1.54 ÷ 2.35 1.68 ÷ 1.81 1.15 ÷ 10.82 1.44IBBN 1.21 ÷ 1.32 1.23 ÷ 1.60 1.17 ÷ 1.37 1.02 ÷ 1.26 1.22 ÷ 1.42 1.60 ÷ 1.23 2.40MOX slate 0.98 ÷ 1.35 0.87 ÷ 1.19 0.65 ÷ 1.09 0.69 ÷ 1.07 0.78 ÷ 1.08 0.94 ÷ 1.19 1.39 ÷ 4.65NRDL Zechstein sediments* 1.16 1.74 0.61 0.36 0.29 0.29RGN soft sediments 2.52 ÷ 5.06 2.60 ÷ 1.94 2.23 ÷ 2.61 5.14 5.83RUE limestone 1.83 ÷ 1.53 1.16 ÷ 1.52 0.93 ÷ 1.31 0.81 ÷ 1.61STU hard marls 0.92 ÷ 1.38 0.86 ÷ 1.27 1.05 ÷ 1.20 1.82 ÷ 1.14 1.47 ÷ 1.24 0.56 ÷ 2.70TNS quartzite 0.88 ÷ 2.00 0.70 ÷ 1.42 0.72 ÷ 1.28 1.09 ÷ 1.41 0.79 ÷ 1.36 0.74 ÷ 1.65 0.79 ÷ 1.77 1.11 ÷ 1.48

UBBA Zechstein sediments* 0.90 ÷ 17.77 0.76 ÷ 1.00 0.40 ÷ 1.76 0.40 ÷ 2.24 0.40 ÷ 1.33 0.42WET gneiss 0.89 ÷ 1.68 0.69 ÷ 1.17 0.54 ÷ 1.10 0.63 ÷ 1.17 0.78 ÷ 1.17 1.54 ÷ 1.18 2.16 2.09

Chapter 2

Monte-Carlo simulations of elasticwave energy in the presence ofdeterministic structure

29

30 CHAPTER 2. MONTE-CARLO SIMULATIONS

Abstract

In this chapter a short introduction into the transfer problem of elastic energy inheterogeneous media is given. A comparison of different strategies for Monte-Carlosimulations leads to the description of an algorithm that was developed to solvethe elastic radiative transfer equation in a model that is close to the situation ofthe earth crust. The algorithm treats anisotropic scattering including conversionsbetween compressional and shearing modes. Deterministic large scale structure isincorporated in terms of a free surface and an interface at depth that reflect andconvert and in the case of the interface also transmit seismic energy. The horizontalinterface separates an infinite half space which may have a vertical velocity gradientfrom the overlaying layer of constant velocity.

INTRODUCTION 31

2.1 Introduction

In chapter 1 radiative transfer theory is applied in the approximation of multipleacoustic and isotropic scattering. For the analysis of the S-wave coda and the es-timation of basic medium and source/site parameters this approximation is useful.The main shortcoming for practical applications is the inability of this approxima-tion to model the P-wave and its coda. From the theoretical point of view, thelimitation to acoustic and isotropic scattering means a severe simplification of thephysics of wave propagation through a heterogeneous elastic medium. In this chap-ter an algorithm is introduced that solves the radiative transfer equation with lessrestrictive approximations. The fundamental improvement is the simulation of elas-tic energy with conversion between P- and S-modes. This can be done in a modelthat contains velocity gradients and discontinuities.

The radiative transfer of elastic energy in a model with deterministic velocitystructure is perhaps the most accurate approach to the description of the scat-tering of seismic waves in the earth’s crust, that is available. Numerous approx-imation to this situation can be found in the literature. Most simple approachesassume isotropic single scattering of acoustic energy an infinite medium (e.g. Akiand Chouet, 1975). Isotropic means that the there is no preferred direction of scat-tering. The acoustic approximation limits the scattering to S-waves. P-waves andconversions are neglected. The single scattering approach assumes that waves arescattered only once on their way from the source to the receiver. Since there isstrong evidence for multiple scattering of seismic waves in the crust (Hennino et al.,2001; Larose et al., 2004), this approximation strongly limits the applicability of thisapproach. The diffusion approximation (Dainty and Toksoz, 1981; Wegler and Luhr,2001) that assumes strong scattering is the other extreme, that is only applicablefor late lapse times of crustal seismic records. Multiple scattering is suitable for theintermediate range of moderate scattering (e.g. Hoshiba, 1991).

The approximation of isotropic scattering was replaced by non-isotropic scat-tering by Hoshiba (1995). The modeling of seismogram envelopes using elastic scat-tering at distributed spheres is due to Margerin et al. (2000). This was extendedto random heterogeneity by Shearer and Earle (2004) and Przybilla et al. (2006).Deterministic structure in terms of an interface was first modeled by Margerin et al.(1998) for acoustic waves.

Here, everything is brought together and an algorithm is presented that is ableto synthesize seismogram envelopes for multiple non-isotropic scattering of elasticwaves at randomly distributed heterogeneities in deterministic structure. Before thedetails of the algorithm’s components are discussed, I present the radiative transferequation in the from that is used by the algorithm and briefly compare differentMonte-Carlo techniques for their solution.


2.2 Elastic radiative transfer theory

The central quantity of transfer theory is the radiance or specific intensity I(n, r).It measures the energy flux from a radiating surface into a unit solid angle arounddirection n. This angular spectrum of the energy flux may be a local quantity anddepend on position r. Alternatively to the energy flux described by the radianceone may work with the local and directionally resolved energy density a(n,x) thatis related to the radiance by

a(n,x) =1

vI(n,x) (2.1)

with v being the velocity of the ray. Here, and in all what follows I implicitlyassume dependence of the radiance or related quantities on time and frequency.This is justified because on one hand the energy transport in seismology is alwaysnon-stationary which implies an intrinsic time dependence. On the other hand thepropagation medium is stationary which ensures that frequencies do not change withtime and remain uncoupled.

In the classical book of Chandrasekhar (1960), the radiative transfer equationthat describes the radiance is introduced phenomenologically to describe light prop-agation through atmospheres. The phenomenological view of the transport problemis based on considerations of energy conservation. Regarding the energy flux indirection n at a fixed point in space, the energy balance can be expressed as:

variation of energyflux in direction n

= −decrement of energy flux due to ab-sorption and scattering from directionn into n′

+increment of energy flux in direction ndue to scattering from directions n′ inton

+increment of energy flux due to radia-tion of sources in direction n

These heuristic considerations are physically transparent and seem intuitivelyconvincing without further explanations. For several decades this field developedwithout complete physical foundations and independently of wave physics. A rigor-ous statistical substantiation of this concept was achieved when the direct connec-tion between the angular Fourier transform of the radiance and the spatial coherencefunction of the wave field was discovered. The coherence function defined as

Γ(x, ρ, t) = 〈u(x + ρ, t)u∗(x, t)〉 (2.2)

for a monochromatic wave field u(x, t) = A(x)ei(knx−ωt) is given by the angularFourier transform of the radiance (Apresyan and Kravtsov, 1996, p. 22):

Γ(x, ρ, t) ∝∫I(n,x, t)eiknρdΩn. (2.3)

ELASTIC RADIATIVE TRANSFER THEORY 33

This makes the radiance somewhat like the angular spectrum of the wave field.Finally this link between the radiance and its “wave-origin” lead to the derivationof the radiative transfer equations for elastic waves by Ryzhik et al. (1996) andWeaver (1990). This is a system of two coupled integro-differential equations for theP-wave Intensity and the 2 × 2 coherence matrix for the S-waves. The coherencematrix I takes into account the two S-polarizations with degenerate dispersion lawsand their phase relation. It is defined as (Apresyan and Kravtsov, 1996, p. 72)

I =

[Iθθ IθφIφθ Iφφ

]= c

[〈AθA∗θ〉〈AθA∗φ〉〈AφA∗θ〉〈AφA∗φ〉

](2.4)

where I define the directions of the components indicated by the subscripts θ and φas the directions of increasing angles θ and φ of the polar coordinate system shownin figure 2.2. Here Aθ and Aφ are the components of the S-wave amplitude and c isa proportionality coefficient. Alternatively the polarization can be described by theStokes vector P . Following Margerin et al. (2000) I define P including the P-waveintensityIP as

P =

IP

ISθISφUV

=

IP

IθθIφφ

Iφθ + Iθφi(Iφθ − Iθφ)

= c

〈|AP |2〉〈|Aθ|2〉〈|Aφ|2〉

2Re〈AθA∗φ〉2Im〈AθA∗φ〉

. (2.5)

with c being a vector of proportionality constants with components c1 = ρω2α/2and c2...5 = ρω2β/2. ISθ and ISφ correspond to the intensities of SV- and SH-wavesrespectively. The parameters U and V can be written in terms of the phase differenceδ between the two S-polarizations in θ and φ directions:

U = 2c4〈AθAφ cos(δ)〉 (2.6)

V = 2c5〈AθAφ sin(δ)〉 (2.7)

The Stokes vector provides a full description of the state of polarization in-cluding linear and elliptic polarizations. The character of a seismic shear source asdetailed in section 2.4.3 and the structure of the scattering coefficients (sec. 2.4.4)allows to simplify the treatment of polarizations. A seismic shear dislocation pointsource emits completely linearly polarized waves. For these waves the phase shift δis zero. Also during scattering events the two polarizations remain in phase. Thisleads me to an approximation of the Stokes vector. In the following I assume thatS-waves are completely linearly polarized and no elliptic polarization occurs. Thisallows to omit the V -component of the Stokes vector which is zero for linear polarizedwaves (cf. equation 2.6 for δ = 0). The U -component is determined by the followingrelation that holds for completely polarized waves (Apresyan and Kravtsov, 1996,p. 74)

(ISθ + ISφ )2 = (ISθ − ISφ )2 + U2 + V 2. (2.8)

For V = 0, U is defined by ISθ and ISφ . This approximation allows to write thetransport equations as three coupled scalar equations for IP , ISθ and ISφ . They fol-low directly from the energy balance in a multi-mode system (c.f. Apresyan and


Kravtsov, 1996, p. 84) and describe three propagating modes. The transport equa-tion for IP reads:

(∂

α0∂t+ n · ∇

)IP (n, r) = −

(gPP0 + gPS0 +

ω

α0IQP

)IP (n, r) (2.9)

+

∫

4π

gPP (Θ)IP (n′, r)dΩn′

+

∫

4π

gSP (Θ,Φθ)ISθ (n′, r)dΩn′

+

∫

4π

gSP (Θ,Φφ)ISφ (n′, r)dΩn′.

Here I used the following notation: α0 is the mean P-wave velocity, gAB0 are thetotal scattering coefficients given by equations 2.47, and IQP is the quality factor ofintrinsic P-wave attenuation (section 2.4.5). The scattering coefficients gAB(Θ,Φ)are discussed in section 2.4.4. These coefficients are assumed to be given in the localscattering coordinate system, i.e. Θ is the scattering angle between n and n′ and Φ ismeasured between the polarization direction of the incident wave and the scatteringplane spanned by n and n′. Φ is different for ISφ and ISθ . Here dΩn′ = sin(Θ)dΘdΦ.

In a similar notation the transport equations for ISθ and ISφ read

(∂

α0∂t+ n · ∇

)ISθ (n, r) = −

(gSS0 + gSP0 +

ω

β0IQS

)ISθ (n, r) (2.10)

+

∫

4π

gSSl (Θ,Φθ)ISθ (n′, r) cos2(ψ)dΩn′

+

∫

4π

gSSr (Θ,Φθ)ISθ (n′, r) sin2(ψ)dΩn′

+

∫

4π

gSSl (Θ,Φφ)ISφ (n′, r) cos2(ψ)dΩn′

+

∫

4π

gSSr (Θ,Φφ)ISφ (n′, r) sin2(ψ)dΩn′

+

∫

4π

gPS(Θ)IP (n′, r) cos2(ψ)dΩn′.

ELASTIC RADIATIVE TRANSFER THEORY 35

XZ

Y

ΘΦθΦφISlISr

ψ

Figure 2.1: Illustration of the S-polarization directions in the radiative transfer equations.Angles are identified by color of the arcs. The red and blue arrows indicate directions n ′

and n respectively. The bluish lines with the ISl and ISr labels indicate the directions ofpolarization of the parallel and perpendicular intensity components respectively

and(

∂

α0∂t+ n · ∇

)ISφ (n, r) = −

(gSS0 + gSP0 +

ω

β0IQS

)ISφ (n, r) (2.11)

+

∫

4π

gSSl (Θ,Φθ)ISθ (n′, r) sin2(ψ)dΩn′

+

∫

4π

gSSr (Θ,Φθ)ISθ (n′, r) cos2(ψ)dΩn′

+

∫

4π

gSSl (Θ,Φφ)ISφ (n′, r) sin2(ψ)dΩn′

+

∫

4π

gSSr (Θ,Φφ)ISφ (n′, r) cos2(ψ)dΩn′

+

∫

4π

gPS(Θ)IP (n′, r) sin2(ψ)dΩn′.

The angle ψ is measured between the scattering plane and the vertical plane thatcontains the vector n. This form of the transfer equations is unusual but it bestrefers to the way I solve the equations with the Monte-Carlo technique. There aretwo rotations hidden in the integrands of these equations. The first corresponds tothe notation of the scattering coefficients in the local frame with orientation of thepolar axis along n′ and Φ measured from the direction of polarization. This resultsin the use of the two angles Φθ and Φφ for the two polarizations which of courseobey the relation Φθ = π/2 + Φφ. Figure 2.1 illustrates this situation.

The second rotation is done by the sin2(ψ) and cos2(ψ) terms which projectthe intensities ISr (n) and ISl (n) given in the directions perpendicular and parallelto the scattering plane respectively onto the θ and φ directions. The sin2(ψ) and


cos2(ψ) terms correspond to the rotation matrix for the Stokes vector (Ishimaru,1978). These two rotations transfer the scattering matrix for intensities (given by itscomponents gAB) into the Muller matrix (Margerin et al., 2000). I do not explicatethe Muller matrix because I evaluate the scattering coefficients in the local frame inthe Monte-Carlo simulations.

2.3 Monte-Carlo techniques to solve the radiative

transfer equation

Monte-Carlo techniques are widely used to solve transport problems. Applicationsrange from light propagation through atmospheres, and neutron scattering to globalillumination computations in computer graphics. In seismology there are two differ-ent strategies to solve the transfer equation with Monte-Carlo methods. The firstmethod that one could refer to as “probability summation method” was first usedby Hoshiba (1991). The idea of this approach is to simulate the propagation ofparticles through the medium and evaluate the probability that a particle is scat-tered towards a receiver at every scattering event. This probability which dependson the scattering coefficient (section 2.4.4), the geometry of the incident ray, andthe direction to the receiver, is multiplied by a factor that accounts for attenuationand geometrical spreading. The resulting quantity is the contribution to the energydensity at the site of the receiver at the potential arrival time of the particle. Soeach scattering event contributes to the energy density estimate at the receiver. Thismethod allows to measure the energy density at point-like receivers which is usefulif there is weak symmetry in the model.

The second method for Monte-Carlo simulations in seismology was first usedby Gusev and Abubakirov (1987) and could be called “particle counting method”.Here the particles are also propagated through the medium but do only contribute tothe energy density estimate at the receiver if they pass through a finite size volumeplaced around the receiver. This is very efficient if there is some symmetry in themodel that allows to use a for instance a spherical shell as receiver volume. Forweak symmetry that occurs if realistic shear sources are simulated, the efficiencydecreases drastically because many particles will delay the simulation that neverpass through the small receiver volume.

In a first approach I implemented the probability summation method. It turnedout to have a major disadvantage.

The principle source of the fluctuations in the energy density estimates is therandomness of the particles’ trajectories. In the probability summation method theprobability of the particle to be scattered towards the receiver and the geometricalspreading between the scattering point and the receiver are additional sources offluctuations. For strongly anisotropic scattering these fluctuations increase dramat-ically and appear to be non-Gaussian. The envelopes have a relatively smooth lowerlevel but irregular spikes which strongly exceed this surrounding level. Looking ata small time window in the coda part of the envelope for increasingly larger num-bers of particles one observes the following. Not only the standard deviation of the

OUTLINE OF THE RadTrans ALGORITHM 37

fluctuations decreases which is expected for stronger averaging but also the meanof the envelope increases. This behavior is caused by the unbounded variance ofthe next-point estimator (Lux and Kolbinger, 1991, p. 386). Despite the fact thatthe energy density estimate converges towards the correct expectation (Lux andKolbinger, 1991, p. 386) this bad convergence characteristic renders the probabilitysummation method inapplicable for strongly non-isotropic scattering.

2.4 Outline of the RadTrans algorithm

With the RadTrans algorithm I follow the approach of Gusev and Abubakirov (1996)and perform a particle counting. In this case the radiance I(n, r) is modeled by thenumber density of particles N(n, r) that are located at position r and move intodirection n. Let us define ∆N(n, r) as the change of particle density along the rayin direction n, i.e. ∆N(n, r) = N(n, r+n∆l)−N(n, r). Noting that the differentialoperator on the left hand side of transfer equations 2.9, 2.10 and 2.11 is the totalderivative along the ray (Margerin, 2005):

d

dl≡ d

vdt=

∂

v∂t+ n · ∇ (2.12)

I can reformulate these equations in terms of ∆N(n, r). This clarifies the relationbetween the transfer equations and the Monte-Carlo modeling. I demonstrate thisfor the case of P-intensity IP . Substituting ∆NP (n, r)/∆l for dIP (n.r)/dl (accord-ing to equation 2.12) in equation 2.9 and NP for IP I obtain after multiplicationwith ∆l and division by NP (c.f. Gusev and Abubakirov, 1996)

∆NP (n, r)

NP (n, r)= −∆l

(gPP0 + gPS0

)−∆l

ω

α0IQP

(2.13)

+∆l

N(n, r)

[∫

4π

gPP (Θ)NP (n′, r)dΩn′

+

∫

4π

gSP (Θ,Φ)NSθ (n′, r)dΩn′

+

∫

4π

gSP (Θ,Φ)NSφ (n′, r)dΩn′

].

This is a probabilistic interpretation of equation 2.9. The left hand side of equa-tion 2.13 is the fractional change of the number of particles (that move in directionn) on a path element ∆l along a ray with direction n. The right hand side has twoterms representing loss of particles. The first describes the fraction of particles thatchange the direction and leave the ray. It depends on the total scattering coeffi-cients gAB0 . The probability of the new direction in which these particles continueto propagate is determined by the scattering coefficient of the accordant scatteringtype:

ΠAB(Θ,Φ) =gAB(Θ,Φ)

4πgAB0

. (2.14)


The fraction of particles that die off because of intrinsic attenuation is described bythe second term containing ω/(α0

IQP ).

The last term of equation 2.13 represents the gain of particles. It depends onall particles at position r independent of their directions n′ and polarizations. Thisgain is balanced by the scattering-loss of rays with direction n′. Thus, the first andlast term of equation 2.13 are closely related and ensure the energy conservation ofthe scattering process.

During the Monte-Carlo simulations a large number of particles has to be sim-ulated to obtain a smooth estimate of N(n, r). But the particles are simulatedsuccessively. For a single particle only the first two terms of equation 2.13 are im-portant. The third term will be automatically modeled during the simulation ofother particles. This leads to the core of the Monte-Carlo modeling that is detailedin the next sections.

The coordinate system in which the simulations are performed is explained insection 2.4.1. Section 2.4.2 describes how the particles are characterized in the sim-ulation. After the particles are initialized i.e. emitted from a source (section 2.4.3)they propagate through the medium in small time steps ∆t in which they move apath ∆l depending on the velocity. With a certain probability depending on the totalscattering coefficients the particles are scattered into new directions (section 2.4.4)and with another probability they are eliminated or their weight is reduced becauseof intrinsic attenuation (section 2.4.5). If the particles are neither scattered noreliminated they continue to propagate in their initial directions. The motion of theparticles on the small path elements ∆l is governed by ray theory which includesreflection/conversion and transmission at interfaces and bending of paths in thepresence of velocity changes (section 2.4.6). The recording process is detailed insection 2.4.7. An auxiliary section describes how random numbers are generatedfrom the given probability distributions for scattering angles and source radiationdirections (section 2.4.8).

2.4.1 The coordinate system

For the recording and the propagation of the particles it is convenient to describetheir position and velocity in Cartesian coordinates. The treatment of reflection,transmission and conversion at horizontal interfaces and the treatment of the scat-tering processes requires knowledge of the polar and azimuthal angles of the prop-agation direction. The orientation of the sources and the position of the seismicstation in turn requires a relation to geographical coordinates. In figure 2.2 the re-lations between the Cartesian-, the polar- and the geographical coordinates as theyare used in the algorithm are illustrated. The Cartesian system is oriented with theZ coordinate increasing downwards. X and Y point north- and eastwards respec-tively. The azimuthal angle φ is measured eastwards from north and the polar angleθ is measured from the vertical downwards.


Y (east)

X (north)

Z (depth)

θφ

Figure 2.2: Illustration of the relations between the geographical, the Cartesian and thepolar coordinate systems as they are used in the algorithm. Angles are identified by color.

2.4.2 Description of a particle

For a complete description of a particle the following parameters are needed:

• Coordinates in space and time (x, y, z, t)

• Velocity vector (vx, vy, vz)

• Stokes vector (IP , ISθ , ISφ , U, V ) describing the polarization.

For convenience during the simulation the following parameters are stored addition-ally:

• Mode of the particle

• Direction of velocity vector in polar coordinates (vθ, vφ).

• Polarization angle ψ

Knowledge of the propagation direction in polar coordinates and of the po-larization angle, defined as the angle between the polarization axis of a linearlypolarized S-wave and the vertical plane that contains the ray path, is useful for thetreatment of interface reflections/conversions/transmissions and the transformationof scattering angles from the local into the global coordinate system.

2.4.3 Particle emission from a point shear dislocation

For the emission of particles, a description of the source in terms of probabilitydensities for radiation directions and particle type with polarization is needed.


A seismic source can be represented as a point shear dislocation if the sourcedimensions are small compared to wavelength and distance. Mathematically thiscan be described by a double couple of forces (Aki and Richards, 1980b, p. 113).The far field displacement vector of such a source is

u(x, t) = BPr er

2√15

M(t− r/α)

4πρα3r︸︷︷︸P−radiation

+[BSθ eθ +BS

φeφ]√2

5

M(t− r/β)

4πρβ3r︸︷︷︸S−radiation

(2.15)

where x = (r, θ, φ) is the position vector in spherical coordinates, M is the timederivative of the moment time function and er, eθ and eφ are unit vectors in thedirections of the spherical coordinates. For a shear source with fault plane normalto the Cartesian (1, 0, 0) direction and slip in the direction (0, 1, 0) the radiationpattern in 2.15 are given by (Sato and Fehler, 1998):

BPr =

√15

2sin2(θ) sin(2φ)

BSθ =

1

2

√5

2sin(2θ) sin(2φ) (2.16)

BSφ =

√5

2sin(θ) cos(2φ).

Here subscripts denote polarization direction and superscripts indicate wave type.The probability distributions of the radiation directions are given by the squares ofthe radiation pattern:

ΠPr = sin4(θ) sin2(2φ)

ΠSθ = sin2(2θ) sin2(2φ) (2.17)

ΠSφ = sin2(θ) cos2(2φ)

The radiation pattern are normalized such that

∮ ∣∣BPr

∣∣2 dΩ =

∮ [∣∣BSθ

∣∣2 +∣∣BS

φ

∣∣2]dΩ = 4π. (2.18)

By integrating the energy flux density vρ|u|2 with v being either P-wave velocityα or S-wave velocity β over a spherical shell and time one can obtain the radiated P-and S-energy respectively (c.f. Sato and Fehler, 1998, p. 149). The ratio of radiatedP- and S-energy is given by

W S

W P=

3

2γ5. (2.19)

The ratio between the energies radiated in the two different S polarizations can


be obtained by integrating e.g.∣∣BS

φ

∣∣2 over the whole solid angle:

∫ 2π

0

∫ π

0

∣∣BSφ

∣∣2 sin(θ)dθdφ =5

2

∫ 2π

0

∫ π

0

sin3(θ) cos2(2φ)dθdφ

=5

2

[(−1

3sin2(θ)− 3

2

)cos(θ)

]π

0

[π

2+

1

8sin(4φ)

]2π

0

=10

3π (2.20)

Using this together with the normalization condition eq. 2.18 I find the desired ratio:

W Sθ

W Sφ

=1

5(2.21)

Equations 2.17, 2.19 and 2.21 allow to simulate the energy radiation of a shearsource with fault normal in (1,0,0) direction and slip direction (0,1,0). For an ar-bitrarily oriented source a rotation matrix is calculated from the given strike, dipand rake angles. The rotation is first performed for the rake around the axis along(1,0,0) then for dip around (0,1,0) and finally for strike around (0,0,1). Both thevelocity and polarization vectors are multiplied with this rotation matrix to accountfor the specific source orientation.

2.4.4 Scattering

The most interesting and complicated part in this algorithm is the simulation of scat-tering of elastic vector waves by distributed heterogeneities. I follow the derivationby Sato and Fehler (1998).

Random medium

Let’s assume a certain property v of the medium is not constant in space but dependson position x. This property may be the wave velocity, the density or a Lamecoefficient. The field v(x) can be written as the sum of a constant mean field v0 anda perturbation δv:

v(x) = v0 + δv(x) = v0[1 + ξ(x)]. (2.22)

ξ(x) is the fractional fluctuation. v0 is chosen such that the average, denoted by〈. . .〉, of an ensemble of heterogeneous media is:

v0 = 〈v(x)〉 and 〈ξ(x)〉 = 0 (2.23)

The random function ξ(x) is best described by its spectral properties. For thisI introduce the auto-correlation function (ACF) as

R(x) ≡ 〈ξ(y)ξ(y + x)〉. (2.24)

In the following it is assumed that the random function ξ is isotropic and homo-geneous which means that R does only depend on the distance r between y and x(r ≡ |y − x|).


The strength of the fluctuation is given by the mean square fractional fluctua-tion

ε2 ≡ R(0) = 〈ξ(x)2〉 (2.25)

and the spatial variation of the fluctuations is characterized by the correlation lengtha. The Fourier transform of the ACF R gives the power spectral density function(PSD). For isotropic and homogeneous random media the PSD is a function of themagnitude of the wavenumber vector m ≡ |m|. The most commonly used ACFsand their corresponding PSDs are given below.

• Gaussian ACF

R(r) = ε2e−r2/a2

and PSD(m) = ε2√π3a3em

2a2/4 (2.26)

This ACF can be used for media that are poor in short wavelength components.

• Exponential ACF

R(r) = ε2e−r/a and PSD(m) =8πε2a3

(1 + a2m2)2(2.27)

The exponential ACF is much richer in short wavelength heterogeneities as itdecays for large wavenumbers as a power low with exponent −4.

• von Karman ACF

R(r) =ε221−κ

Γ(κ)

(ra

)κKκ

(ra

)for κ = 0 ∼ 0.5 (2.28)

PSD(m) =8π3/2ε2a3Γ(κ+ 3/2)

Γ(κ)(1 + a2m2)κ+ 32

(2.29)

Here Γ is the gamma function and Kκ is the modified Bessel function of secondkind and order κ. The von Karman ACF with κ = 0.5 is equivalent to theexponential ACF. For smaller κ the amount of short wavelength fluctuationsincreases even further.

Scattering coefficients

Scattering of distributed heterogeneities in a random medium can be described bythe scattering coefficients. The formal derivation starting from the wave equationin a random medium can be found in Sato and Fehler (1998, p. 95-105). Beginningwith a wave of unit amplitude and angular frequency ω incident on an inhomogeneitywhich is localized at the origin one can define the scattering amplitude F as:

usc(x, t) =ei(kr−ωt)

rF. (2.30)

It relates the amplitude of a wave usc(x, t) that is scattered in a certain directionto the amplitude of the incident wave. There is one scattering amplitudes for each


type of scattering. Comparing the scattered energy flux Jsc to the incident one (J0)leads to the definition of the scattering cross section:

dσAB

dΩ≡ Jscr

2

J0=vBvA

∣∣FAB∣∣2 (2.31)

Here A and B indicate the wave mode which may be either P or S. vA and vB arethe mean velocities of the respective modes. If the inhomogeneities are not localizedbut distributed over a larger volume the scattering coefficient may be defined as 4πtimes the ensemble average of the scattering cross section of a region of size L.

gAB ≡ 4π1

L3

⟨dσAB

dΩ

⟩= 4π

1

L3

⟨∣∣FAB∣∣2⟩

(2.32)

Thus the scattering coefficients relate the incident energy flux to the angular de-pendent energy flux that is scattered from distributed heterogeneities in a region ofunit size around the origin.

The main assumptions that are used for the derivation of the scattering coeffi-cients are the following:

• The heterogeneities can be described by eq. 2.22 with |ξ| 1 indicating thatthe fluctuation must be small compared to the average value of the quantity.

• In order to be able to use a first order perturbation method to decomposethe wave field into a scattered and a homogeneous unscattered componentit is assumed that the scattered field is much smaller in amplitude than theincident field. This approximation means that the scattering must be weakand is known as the Born approximation.

• To solve the wave equation for the scattered component of the wave field theGreen function in homogeneous space is used.

• It is further assumed that observations are made in the far field (r a)and the Fraunhofer zone (r a2k/π) of the scatterer. This means that thedistance between scatterer and receiver on one hand and the wavelength onthe other hand have to be large compared to the spatial length scale of thefluctuations.

• The fluctuations of the wave velocities are equal. They are correlated with themass density of rock as stated by Birch’s law (Birch, 1961). These correlationsreduce the fluctuations to one single fractional fluctuation that is common tothe wave velocities and the mass density.

ξ(x) ≡ δα(x)

α0=δβ(x)

β0=

1

ν

δρ(x)

ρ0(2.33)

α and β are the velocities of P- and S-waves respectively and ρ is the massdensity. The parameter ν ranges between 0.78 < ν < 0.8 in typical lithosphere(Sato and Fehler, 1998, p. 101).


These are the formal assumptions and approximations used to derive the spe-cific form of the scattering coefficients that is used here. The range of applicabilityof the radiometric concept is much wider. Przybilla et al. (2006) showed that evenwith scattering coefficients derived under these assumptions the correspondence be-tween the radiative transfer equation and the full wave equation goes far beyondwhat would be expected from the above conditions.

In spherical coordinates (Θ,Φ) the scattering coefficients are given by Sato andFehler (1998, p. 104):

gPP (Θ) =l4

4π

∣∣XPP (Θ)∣∣2 PSD

(2l

γ0sin

(Θ

2

))(2.34)

gPS(Θ) =1

γ0

l4

4π

∣∣XPS(Θ)∣∣2 PSD

(1

γ0

√1 + γ2

0 − 2γ0 cos(Θ)

)(2.35)

gSP (Θ,Φ) = γ0l4

4π

∣∣XSP (Θ,Φ)∣∣2 PSD

(1

γ0

√1 + γ2

0 − 2γ0 cos(Θ)

)(2.36)

gSSl (Θ,Φ) =l4

4π

∣∣XSSl (Θ,Φ)

∣∣2 PSD(

2l sin

(Θ

2

))(2.37)

gSSr (Θ,Φ) =l4

4π

∣∣XSSr (Θ,Φ)

∣∣2 PSD(

2l sin

(Θ

2

))(2.38)

Θ measures the angle between the initial direction and the direction of the scat-tered particle. This angle is usually called scattering angle. Φ measures the anglebetween the polarization direction of the incident particle and the scattering planespanned by the incident and scattered directions. Superscripts indicate the typeof scattering. SP for example stands for the scattering of a S-wave into a P-wave.Subscripts r and l indicate the polarization of the scattered particle with respect tothe scattering plane. r and l stand for perpendicular (normal to scattering plane)and parallel (in the scattering plane) respectively. This distinction is only neces-sary for SS-scattering since the polarization in the other cases is non-ambiguous.The polarization of P-waves is radial (in the direction of propagation) and S-wavesthat were generated by scattering from P-waves have parallel polarization in thescattering plane. The argument of the power spectral density function is the socalled exchange wavenumber expressed in terms of the scattering angle and the S-wavenumber l. The terms denoted by X are the basic scattering patterns. They arethe same for all types of media. Modifications of the angular patterns of the scatter-ing coefficients are solely introduced by the power spectrum of the heterogeneities.


The scattering patterns are given by Sato and Fehler (1998, p. 102):

XPP (Θ) =1

γ0

[ν

(−1 + cos(Θ) +

2

γ20

sin2(Θ)

)− 2 +

4

γ20

sin2(Θ)

](2.39)

XPS(Θ) = − sin(Θ)

[ν

(1− 2

γ0cos(Θ)

)− 4

γ0cos(Θ)

](2.40)

XSP (Θ,Φ) =1

γ20

sin(Θ) cos(Φ)

[ν

(1− 2

γ0cos(Θ)

)− 4

γ0cos(Θ)

](2.41)

XSSl (Θ,Φ) = cos(Φ) [ν (cos(Θ)− cos(2Θ))− 2 cos(2Θ)] (2.42)

XSSr (Θ,Φ) = sin(Φ) [ν(cos(Θ)− 1) + 2 cos(Θ)] (2.43)

The scattering coefficients 2.34 to 2.38 define the probability distributions forthe directions of the scattered particles in a local frame with respect to the initialdirection and polarization.

Transformation of angles

The new propagation and polarization directions of a scattered particle are selectedin a local coordinate system defined by the incident direction and the direction ofpolarization. For further computations the new directions have to be transformedinto the global coordinate system described in sec. 2.4.1. The geometry of thistransformation is depicted in figure 2.3. Let the direction of the incident particle begiven by the polar coordinates (θi, φi). The polarization of the incident particle isgiven by ψi which is the angle between the polarization axis and the vertical planethat contains the incident ray. The two angles that describe the scattering processare the scattering angle Θ and Φ, the angle between the scattering plane and theinitial polarization direction. To express the polar angle θs and azimuth φs of thescattered particle in terms of the direction of the incoming particle (θi, φi) and thescattering and polarization angles (Θ,Φ, ψi) one can use spherical trigonometry. Thepolar angle of the scattered particle θs can be expressed as

cos(θs) = cos(θi) cos(Θ) + sin(θi) sin(Θ) cos(ψi + Φ). (2.44)

The azimuth φs of the scattered particle is given by

cos(φi − φs) =cos(Θ)− cos(θi) cos(θs)

sin(θi) sin(θs). (2.45)

To find the polarization angle of the scattered particle the angle ψs between thescattering plane and the vertical plane that contains the scattered ray is calculatedusing

sin(ψs) =sin(φi − φs)

sin(Θ)sin(θi). (2.46)

The polarization direction of the scattered particle is either ψs+π as given by eq. 2.46if the polarization is parallel to the scattering plane or ψs + π/2 if the polarizationis perpendicular.


XZ

Y

θi

φi

θs

φs

ΘΦψ

iψ

spol

Figure 2.3: Illustration of the scattering and the polarization angle and their relation tothe global polar coordinate system. Red and blue arrows indicates the directions of theincoming and scattered particle respectively. The polarization direction of the incomingS-particle is indicated by the bold black line.

Total scattering coefficients

If a particle undergoes scattering, the probability distribution of the direction inwhich the particle continues to propagate is given by the scattering coefficients 2.34to 2.38. The integral of the scattering coefficients over the solid angle gives the totalscattering coefficient gAB0 , which defines the probability that a certain scatteringevent from mode A to B occurs while the particle moves one unit length (Sato andFehler, 1998, p. 43):

gAB0 ≡ 1

4π

∮gABdΩ =

1

lAB(2.47)

gSS0 is the sum of the total scattering coefficients for the two different polarizations.The inverse of gAB0 is the path lAB , i.e, the average length that a particle of type Apropagates until it is scattered into mode B. The length lA defined as

lA =

(∑

B

gAB0

)−1

(2.48)

is the mean free path of wave mode A. The length LA that a particle of mode Apropagates between two scattering events is determined by the exponential proba-bility distribution

Π(LA) = (1/lA) exp(−LA/lA) (2.49)

.

Random walk

In the previous sections I have outlined how the probability distributions for thescattering events are related to the properties of the medium. Here I combine these


information and explain how it is used to simulate the random walk of the particles.This is largely analogous to Margerin et al. (2000).

The random walk is performed in incremental time steps ∆t. Depending oftheir velocity the particles move paths of length ∆l in this time. The probabilityof an incident particle of mode Mi to be scattered while moving ∆l is according toequation 2.49

Π(Mi) = 1− exp

(−∆l

lMi

). (2.50)

A first uniform random number ε1 ∈ (0, 1) is used to decide whether a particle isscattered by comparing it to Π(Mi). If ε1 < Π(Mi) scattering occurs.

In this case the mode of the outgoing particle Mo is selected according to theconditional probability distribution

Π(Mi|Mo) =gMiMo

0∑Mo

gMiMo0

. (2.51)

It is the probability that a particle is scattered into mode Mo given that the incidentmode is Mi. For the selection of the outgoing mode a second uniform random numberε2 ∈ (0, 1) is used. Here the two S-polarizations are considered different modes,which means that in the case of an incoming S-particle it is decided at this pointwhether the outgoing particle is polarized parallel or perpendicular to the scatteringplane, i.e. whether the new direction is selected using gSSl or gSSr . This choice alsodetermines the polarization of the outgoing particle.

Due to the special form of the scattering coefficients (equations 2.34 to 2.38)that are used to select the new directions of the particle, and the assumption oflinear polarization, the probability distributions of the scattering angle Θ and angleΦ are independent. The probability of the angle Θ under the condition that theincoming particle of mode Mi is scattered into mode Mo is

Π(Mi,Mo|Θ) =

∫gMiMo(Θ,Φ)dΦ

4πgMiMo0

. (2.52)

The probability of the angle Φ is

Π(Mi,Mo|Φ) =

∫gMiMo(Θ,Φ) sin(Θ)dΘ

4πgMiMo0

. (2.53)

For fast evaluation of the probability distributions Π(Mi,Mo|Θ) and Π(Mi,Mo|Φ)during the simulation I generate tables with random numbers drawn from thesedistributions prior to the simulation (section 2.4.8). This means that tables withvalues of Θ are generated for the five types of scattering that can occur, namelyP → P , P → S, S → P S → Sl, and S → Sr. For the angle Φ only a singletable has to be generated with numbers distributed according to sin2(Φ). Numberswith a cos2(Φ) distribution are obtained by adding π/2 to the sin2(Φ) distributednumbers. In the cases of P → P and P → S scattering the values for the angle


Φ are uniformly distributed in [0, 2π). Two further uniformly distributed randomnumbers ε3, ε4 ∈ (0, 1) are used to select values for Θ and Φ from the appropriateprobability tables.

Since the radiative transfer equations in section 2.2 are written as a 3-modesystem, the ISθ and ISφ intensities could be modeled with different particles. ButI decided to allow combinations of ISθ and ISφ , i.e. linear polarization in arbitrarydirections to reduce the need of averaging. This implies that for an incoming S-particle the angle Φ in equation 2.53 is measured from the direction of its polarization(ψi in figure 2.3).

The separate treatment of parallel and perpendicular S → S scattering is donefor the probabilities of Θ and Φ to be independent. This greatly simplifies the simu-lation. Otherwise the probability distribution of Φ would depend on the choice of Θ.The effect is that the outgoing particle has either parallel or perpendicular polariza-tion. The correct polarization of the S → S scattering is obtained after averagingover a number of particles. The correct convergence of this process is assured by thesin2(Φ) and cos2(Φ) terms appearing in the scattering coefficients gSSl and gSSr whichare the same as in the rotation matrix for the Stokes vector (Ishimaru, 1978). Thismeans the Stokes vector that would result from simultaneously treated parallel andperpendicular S → S scattering is, except for a rotation of the coordinates, equalto the average of a number of separately treated S → Sl and S → Sr scatterings inthe same direction (Θ,Φ).

2.4.5 Intrinsic attenuation

The modeling of intrinsic attenuation seems straight forward, but is algorithmicallysomewhat tricky. If the intrinsic attenuation is uniform in space and equal for allmodes it can be accounted for by the multiplication of the envelopes with an ex-ponentially decaying function after the simulation. In the elastic case where twodifferent wave modes with different attenuation properties are modeled, the attenu-ation has to be taken into account during the simulation. There are two possibilitiesfor this. The particles can either have constant energy but die off after a certaintime with a probability depending on the attenuation or the energy of the particlesis decreased continuously as they propagate.

The parameter that is usually used to describe the amount of intrinsic atten-uation is the intrinsic quality factor IQ or its reciprocal quantity the attenuationparameter IQ−1. The attenuation parameter is proportional to the exponential de-cay rate of seismic wave amplitude with distance. In a non-scattering medium theenergy density of a plane wave decays like

EP (r, f) ∝ e−2πrf

α0IQP and ES(r, f) ∝ e

−2πrf

β0IQS . (2.54)

This decay is already contained in the transfer equations 2.9 through 2.11. Inequation 2.13 this decay is only described correctly for infinitesimal ∆l. For finitesize path length ∆l the probability that a particle dies off can be calculated from 2.54to be

Π(∆l) = 1− e−ω∆l

α0IQP and Π(∆l) = 1− e

−ω∆l

β0IQS (2.55)


for P- and S-waves respectively. These are also the factors that the Stokes vector hasto be multiplied with if the particles keep propagating. This is the strategy chosenin the present algorithm because it accelerates the convergence of the energy densityestimates in the later part of the envelopes. However, it requires some considerationsof the numerical accuracy.

Total attenuation parameters Q−1 =I Q−1 +ScQ−1, with ScQ−1 being the scat-tering attenuation parameter, of the earth’s crust were estimated between 0.0005and 0.05 (Sato and Fehler, 1998, fig. 5.1 and 5.2). Assuming that IQ−1 and ScQ−1

are of equal size, one can roughly approximate the maximum decay that can occurduring a simulation. For a large simulation length of 1000 s and an intrinsic qualityfactor of 100 the energy can decrease by factor of 10−41 at 3 Hz whereas it decreasesonly by a factor of 10−4 in a medium with IQ = 1000. Obviously there is a largedynamic range required for the simulation of intrinsic attenuation. In a simulationthat consists of a low IQ crustal layer and an underlying mantle with high IQ theenergy of two particles of which one propagated in the mantle and the other inthe crust can differ by a factor of 1037 after 1000 s. When these two particles arerecorded in the same cell the numerical accuracy of the floating point representationcan lead to a phenomenon known as absorption in computer science. The valueof the larger number does not change if the smaller is added. Because of this thealgorithm works with 64-bit integer numbers that have the highest dynamical rangewith constant precision.

2.4.6 Ray tracing

Once the mode and direction of a particle are fixed, it propagates through themedium according to ray theory (e.g. Aki and Richards, 1980b). This includes theinteraction with interfaces at which the particles can be reflected or transmittedjust like seismic waves. They can also change their mode if conversion occurs. In amedium with constant mean velocity the particles move on straight lines. If thereare gradual changes of the mean velocity in the medium the paths are bent. Thealgorithm is capable of dealing with any depth dependent velocity structure.

Interfaces

If a particle hits an interface between two media characterized by different massdensities and wave velocities it can be reflected, transmitted or converted. Theprobabilities for the particle to undergo either of these events are given by therespective energy reflection coefficients. These will be expressed in the notation ofAki and Richards (1980b). Assuming horizontal interfaces acute and grave accentsare used to indicate up-going and down-going particles respectively. The generalform of an energy reflection coefficient is

Π(A→ B) =ρovo cos(γo)

ρivi cos(γi)

(AB)2

. (2.56)

Here Π(A → B) is the energy reflection coefficient i.e. the probability that anincoming particle of type A with direction · continues in direction · as type B. AB


is the reflection coefficient for displacement amplitudes. A and B stand for the wavetypes P or S and · and · stand for the accents · and · indicating up- or down-goingdirections. For an interface between two solids all 16 combinations of P , S, · and · aredefined. ρ and v indicate mass density and wave velocity respectively. γ is the anglebetween the direction of the ray and the normal of the interface. Subscripts i and oindicate the medium on the side of the incoming or outgoing particle respectively.The velocity, however, does not only depend on the medium but also on the wavetype.

For horizontal interfaces the two different S-wave polarizations can be treatedseparately. For horizontally polarized waves there is no conversion into other modes.The reflection coefficients for displacement amplitudes for this case are given by Akiand Richards (1980b) equation 5.33 p. 139. For vertically polarized S-waves andP-waves which are coupled at the interface the amplitude reflection coefficients aregiven by Aki and Richards (1980b) equations 5.40 p. 144-145. The interaction withthe free surface is a special case of these coefficients for the case of zero velocity inthe upper layer. This causes total reflection for horizontally polarized S-waves. Thedisplacement amplitude coefficients for horizontally polarized S- and P-waves aregiven in equations 5.27, 5.28 (p. 134), 5.31 and 5.32 (p. 136) by Aki and Richards(1980b).

I treat supercritical reflection i.e. S-energy reflection for which the P-reflectionangle would be greater than π/2 as total reflection similar to horizontally polar-ized waves. Doing so I neglect phase shifts between the two S-polarizations uponreflection which could change the polarization type (e.g. from linear to elliptical).

Ray paths in a medium with vertical velocity variations

Ray theory predicts that in a laterally homogeneous medium the particles propagateon curved rays in a vertical plane. In spherical coordinates the azimuthal angle φof the velocity vector remains constant and only the polar angle varies along thepath. The curvature of the ray is determined by the velocity gradient. In order touse the usual notation of Snell’s law I introduce the incidence angle α = π− θ for θas defined in figure 2.2. In this notation Snell’s law states that

sin(α)

v= q = const. (2.57)

Here v is the wave velocity of either P- or S-waves and the constant q is called theray parameter. In order to keep q constant along a ray that passes through regionswith changing velocity the relative change of sin(α) must balance the relative changeof v. I thus have

d sin(α)ds

sin(α)=

dvds

v(2.58)

where s parameterizes the length of the path. Considering the geometry of the raypath in a medium with a vertical velocity gradient this leads to

dαds

cos(α)

sin(α)=

dvdz

cos(α)

v(2.59)


and finally todα

ds=dv

dz

sin(α)

v(2.60)

where the last factor can be recognized as the ray parameter (equation 2.57). Equa-tion 2.60 determines the change of the incidence angle along the ray. For a constantvelocity gradient dv/dz the ray path is a circle.

With equation 2.60 I can approximately model the particle propagation inany depth dependent velocity structure by simply changing the θ component of thevelocity vector after the propagation of a path element ∆l an amount of − dv

dzsin(θ)v

∆l.

2.4.7 Recording of particles

The principal result of a simulation is the number density of particles which gives anestimate of the energy density (equation 2.13). Since seismograms do not containinformation about the direction of seismic waves I do not accumulate informationabout the direction of particles in the simulation. The total number of particlesis counted in a 4-dimensional array. This array represents rectangular boxes in 3-dimensional space with a 4-th dimension for time. Snapshots of the energy field aretaken in time intervals ∆t. Like in a histogram the position of a particle is recordedby adding its energy to the cell of the array in which the particle resides. Thisway I obtain a spatio-temporal estimate of the local energy density as defined byequation 2.1 for directions integrated over the whole solid angle.

2.4.8 Evaluation of probability distributions

A key point in Monte-Carlo simulations is the evaluation of probability distributions.In the present algorithm random numbers have to be generated with probabilitydistributions for the direction of radiation from the source (eq. 2.17) and for thescattering angles (eq. 2.34-2.38). Since, during a simulation a large quantity ofrandom numbers is needed, it is useful for the computational efficiency to generatetables with random numbers prior to the simulation.

I use an iterative trial and reject method to generate random numbers witharbitrary probability distributions. This method does a random sampling of a givenprobability distribution P (x). A set X of N random numbers xi ∈ [xmin, xmax],i = 1, . . . , N having the distribution P (x) is selected as a subset of a uniformlydistributed set Y of random numbers yj ∈ [ymin, ymax], j = 1, . . . ,M . The randomnumber yj belongs to X if zj ≤ P (yj)/Pmax with zj ∈ [0, 1] being another randomnumber from a uniformly distributed set. P (x)/Pmax is the normalized probabilitydistribution where Pmax is the maximum of P (x). As an approximate estimate ofPmax I use Pmax = maxj=1,...,M(P (yj)). This method is similar to the Metropolis-Hastings algorithm (Hastings, 1970) with uniform proposal density.

The probability distributions encountered here depend on two angular vari-ables. But as can be seen in equations 2.17 and 2.34-2.38 the two variables areindependent and the joint two-dimensional probability distribution is the product


of the two one-dimensional distributions. However the trial and reject method asexplained above can be generalized to multi-dimensional probability distributions.

2.5 Conclusions

The algorithm is implemented in C programming language. The simple overall struc-ture of the code allows to parallelize it by splitting the loop for the particles. UsingMPI it can be run in parallel on clusters of any size. The speed of the computationstrongly depends on the parameters of the random medium and it is impossible togive general statements about the performance. The system requirements are de-fined by the size of the 4-dimensional recording grid and the desired smoothness ofthe energy density estimates in space and time. Models for regional wave propa-gation over several hundred kilometers with a simulation length of several hundredseconds require memory of a couple of GB. In combination with inversion proceduresthat depend on the simulation of a large number of models, the application of thisalgorithm certainly requires the use of large computer clusters.

The correctness of the implementation is verified in three different way. Themost simple property that is tested is energy conservation. For the Monte-Carloalgorithm this means that the number of particles in the simulation remains con-stant over time. In the case of zero intrinsic absorption this is confirmed also in thepresence of discontinuities. The second property that more deeply depends on thephysics of the radiative transfer is equipartition. It is verified that the energy ratiobetween S- and P-energy converges to the ratio of 2α3

0/β30 predicted by radiative

transfer theory (Ryzhik et al., 1996). The third test that is performed is a compar-ison with finite difference calculations. No significant or systematic differences areobserved between ensemble averaged energy density envelopes produced with finitedifference calculations in full space and the envelopes produced with my algorithm.

Chapter 3

Elastic radiative transfer theory toexplain the Lg-blockage in thePyrenees

53

54 CHAPTER 3. LG-BLOCKAGE IN THE PYRENEES

Abstract

The RadTrans algorithm that is presented in chapter 2 is well suited to simulateenergy propagation in a model that contains both, deterministic structure and sta-tistically described small scale heterogeneities. With this algorithm the propagationof seismic waves through the Pyrenees is studied. Using a series of seismogramsthe variations of the attenuation properties of the crust along the Pyrenean axis isdemonstrated. The western Pyrenees exhibit strong attenuation of crustal phases.This phenomenon is known as Lg-blockage. Since the velocity structure of the west-ern Pyrenees cannot account for the attenuation, the observation was unexplained.Here I present a model that is able to explain the Lg-blockage. I propose that theLg-blockage is caused by a perturbed body in the crust that differs in the scatteringand attenuation properties from the surrounding crust. The velocities are identical.The results show that it is possible to model the observed Lg-blockage with increasedscattering in the perturbed body. Stronger intrinsic attenuation than the surround-ing crust improves the model further but intrinsic attenuation alone, cannot explainthe observations.

STRUCTURE AND EVOLUTION OF THE PYRENEES 55

3.1 Introduction

The Monte-Carlo algorithm that is introduced in chapter 2 allows to simulate fullseismogram envelopes of regional earthquakes in continental crust. It can be usedto improve existing methods for the estimation of medium or source/site parame-ters like in chapter 1. Compared to such simple methods, the RadTrans algorithmprovides a much more realistic model, but needs more parameters to describe it. Todescribe the elastic scattering in the way discussed in chapter 2, four parametersare needed: the correlation length a, the strength of the fluctuations ε, the qualityfactor of intrinsic absorption Q, and to characterize the type of the auto-correlationfunction in the case of a von Karman medium, the parameter κ. This results ineight unknown parameters to describe the energy propagation in crust and man-tle. Due to this multitude of free parameters their precise estimation with sufficientconfidence appears very difficult.

Here a study is presented that has a different focus. The chapter starts withan introduction to the geology and tectonics of the Pyrenees and describes the ob-servation of anomalous attenuation of seismic waves in some parts of this mountainrange. The special capability of the RadTrans algorithm to simulate energy propa-gation in a realistic crust-mantle structure is then used to test different models fortheir ability to explain the observed attenuation.

3.2 Structure and evolution of the Pyrenees

The geological setting and tectonics of the Pyrenean range is characterized in thereview paper by Choukroune (1992) with the words: “the Pyrenean belt cannot beeasily compared to other belts and is probably unique in its unusual structural styleand evolution”. I therefore restrict to a simplified description which assumes aneast-west orientation of the tectonic units resulting from the collision between theIberian and European plates.

The central part of the Pyrenees is formed by the Paleozoic Axial Zone (PAZ).This zone coincides with the main topographic expression of the Pyrenees. It iscomposed of Hercynian rocks that were reactivated during the Alpine orogeny. Thiszone is confined by the North Pyrenees Fault which marks the boundary to theNorth Pyrenean Zone (NPZ) to the north. The NPZ is composed of Mezosoic flyschdeposits that are locally highly deformed and metamorphosed. Some large Paleozoicoutcrops form the North Pyrenean- and the Basque Massifs. Further north the NPZoverrides the Aquitaine molasse basin along the North Pyrenean Frontal Thrust.

Almost symmetrically the South Pyrenean Zone (SPZ) confines the PAZ to thesouth. The SPZ is composed of Mesozoic and Cenozoic sediments. At the SouthPyrenean Thrust the SPZ overrides the Ebro Basin.

Several seismic refraction and wide-angle reflection profiles revealed the Mohostructure beneath the Pyrenees. There is a jump in Moho depth between the 50 kmthick Iberian crust and the European crust which has a thickness between 29 kmand 35 km.


Irrespective of this general north south structure that is common to the wholerange there are variations along the axis of the Pyrenees. The compressional motionduring the last 65 Ma is related to subduction of the Iberian lower crust beneaththe Aquitaine basin in the eastern and central parts of the Pyrenees (Souriau andGranet, 1995). This observation, made with seismic tomography was confirmed bymagnetotelluric measurements (Pous et al., 1995) and gravity modeling (Vacher andSouriau, 2001). In the western part the shortening has been accommodated in thecrust.

The lower crust in the western part of the NPZ has a much higher P-wavevelocity than the PAZ. Daignieres et al. (1981) speculate that this difference is causedby intrusion of mantle material into the lower crust. This difference is not observed inthe eastern part of the range. The emplacement of lower crustal blocks into the uppercrust of the western Pyrenees has been detected by body wave tomography (Souriauand Granet, 1995) and gravity modeling combined with petrological information(Vacher and Souriau, 2001).

3.3 Data and observations

The geological variations in the Pyrenees along strike have an interesting manifesta-tion in the attenuation properties of crustal waves. There is a pronounced differencebetween the western part and the eastern and central parts of the Pyrenees. Thewestern part strongly attenuates crustal waves (Chazalon et al., 1993). I illustratethis observation using data from the seismic network run by the French AtomicEnergy Commission (CEA) and from the Spanish station PAB run by the InstitutoGeografico National, Spain. Figure 3.1 shows a map of the earthquakes and stations.The bold red line with black dots along the Pyrenees marks the axis used to sort thetraces according to the place where the rays (blue lines) crossed the Pyrenees. Thedots indicate 50 km distance intervals. The corresponding waveforms are shown infigure 3.2. They are band-passed between 2 and 4 Hz, scaled to group velocity, nor-malized for their maximum and sorted for place of intersection with the axis of thePyrenees. The horizontal lines in figure 3.2 indicate the velocities of the dominantregional seismic phases. These are the two crustal phases Pg (vPg = 6.03km/s) andSg (vSg = 3.56) and the two mantle phases Pn (vPn = 8.16) and Sn (vSn = 4.65).The velocities are taken from the velocity model of the Laboratoire de Detection etde Geophysique (LDG).

For wave paths that cross the eastern Pyrenees (east of the 250 km mark infigure 3.2) seismograms are dominated by the crustal phases. The reverberationsin the crust form an extended Sg-arrival of high amplitude, the Lg-wave. The Sn-arrival is barely visible and the Pn amplitude is small compared to the Pg phase.In the western part of the Pyrenees the situation is different. The traces that crossthe Pyrenees west of the 200 km mark show almost no crustal phases. A strongPn-arrival is followed by Sn which usually has the highest amplitude. The Lg-wavecan hardly be recognized in the coda of the Sn-phase. This strong attenuation ofcrustal phases, which is most obvious for the Lg-waves is referred to as Lg-blockage.

DATA AND OBSERVATIONS 57

6oW 4oW 2oW 0o 2oE 4oE 6oE

39oN

42oN

45oN

48oN

51oN

Figure 3.1: Map of the study area with wave paths used in figure 3.2. Black trianglesindicate stations and red stars indicate earthquakes. Red line with black dots shows theaxis of the Pyrenees used to sort the traces in figure 3.2. The black dots indicate 50 kmintervals.

Another illustration of the Lg-blockage is presented in figure 3.3. Part A ofthis figure shows records of six earthquakes that occurred north and south of thePyrenees made at the station SJPF inside the western Pyrenees. Part B showsrecords of the same earthquakes made after the waves passed through the westernPyrenees. A map with the events and the stations is shown in figure 3.4. The clearPg and Lg waves in figure 3.3A demonstrates that crustal phases are generatedby the sources and travel from both sides towards the station inside the Pyrenees.Mantle phases are very weak. After passing through the Pyrenees the shape of theseismograms is again completely different (figure 3.3B). Except for events 1 and6 where a small amplitude increase is visible at the arrival time of the Lg phase,there are no crustal phases. Using the interchangeability of source and receiver, thisexactly reflects an observation made by the analysts at the LDG. Some earthquakesin the western Pyrenees generate Lg-waves that can be recorded at both sides of therange. But Lg-waves that where generated outside the Pyrenees do not propagatethrough.

To generate traces that represent the main characteristics of the undisturbedpropagation through the easter Pyrenees and the propagation through the westernPyrenees with blocked Lg-waves I use the following approach.

I split the collection of traces in figure 3.2 into groups according to the positionof intersection between the ray path and the axis of the Pyrenees. For the first groupI select traces that cross the Pyrenean axis between the 100 km and 200 km marks.I exclude wave paths west of the 100 km mark to not be biased by paths throughthe Bay of Biscay. This group represents the propagation through the westernPyrenees. The second group that represents the eastern part is made up of pathseast of the 300 km mark. The traces are band-passed between 2 and 4 Hz before theenvelopes are computed with a Hilbert transform. The envelopes are normalized fortheir maximum, scaled to group velocity and stacked within each group. Figure 3.5


0 50 100 150 200 250 300 350 400

8.16

6.03

4.65

3.56

3

2

grou

p ve

loci

ty [k

m/s

]

position along axis of Pyrenees [km]

Figure 3.2: Collection of seismograms for wave paths through the Pyrenees. Traces aresorted according to the position of the intersection between their paths and the axis of thePyrenees (red line in figure 3.1). The traces are band pass filtered between 2 and 4 Hz,normalized for their maximum and scales to group velocity.

A8.16 6.03 4.65 3.56 3 2

Event 6

Event 5

Event 4

Event 3

Event 2

Event 1

v [km/s]B

8.16 6.03 4.65 3.56 3 2v [km/s]

Figure 3.3: Seismograms of six earthquakes recorded within (A) the zone that blocksLg-waves and recorded after the waves passed through this zone (B). A map showing thelocations of the events and the stations is shown in figure 3.4.

shows the resulting envelopes. These clearly show the characteristics of Lg-blockage.

There are three mechanism known in the literature that could possibly explainsuch an observation of variable amplitude ratio between crustal and mantle phases.

• The mantle phases propagate in a high Q medium with little attenuationwhereas the crustal phases are stronger attenuated (Bormann, 2002). Thisgenerally causes a distance dependence of the amplitude ratio. Indeed thewave paths through the western Pyrenees are generally slightly longer thanthe paths in the east but there are records with equal or reversed distancerelation. Also these traces with the same source and equal epicentral distance

DATA AND OBSERVATIONS 59

6oW 4oW 2oW 0o 2oE 4oE 6oE 36oN

39oN

42oN

45oN

48oN

51oN

2

1

MFF3

SJPF

4 5

6PAB

Figure 3.4: Map of the study area with wave paths used in figure 3.3. Black trianglesindicate stations and red stars indicate earthquakes.

8.16 6.03 4.65 3.56 3 2v [km/s]

wes

tern

Pyr

enee

s

e

aste

rn P

yren

ees

Figure 3.5: Reference traces from the eastern and western Pyrenees Vertical lines indicatearrival times of seismic phases. The upper black trace represents undisturbed propagationthrough the eastern Pyrenees. It shows clear arrivals of crustal phases. The bottomtrace represents wave paths through the western Pyrenees. Here the crustal phases areattenuated and only the mantle phases show clear arrivals.

show the Lg-blockage.

• The radiation pattern of the source excites strong mantle phases and weakcrustal phases (Bormann, 2002). There are examples that confirm this situ-ation but they seem rather rare. I use several earthquake sources to averagethe observation over different sources in order to minimize this effect.

• The velocity structure in the western Pyrenees causes the Lg-blockage.

The last point is worth a more detailed discussion.The sensitivity of Lg-waves to crustal structure is known from the extinction of

this phase when passing through oceanic crust. Lg-blockage has also been been re-ported from the Barents Sea (Baumgardt, 2001), eastern Mediterranean, the Caspian


Sea, the Black Sea, the Red Sea (McNamara and Walter, 2001; Rodgers et al., 1997)and the Lingurian Sea (Shapiro et al., 1996) and the Alpine range (Campillo et al.,1993).

In most of these cases the attenuation is caused by soft marine sediments. Theyform a layer of very low shear velocity and high attenuation that traps the energywhich propagates in the crust. This phenomenon of Lg-blockage by sedimentarybasins was modeled by Shapiro et al. (1996). Apart from minor details the velocitystructure explains the Lg-blockage. The Lg-blockage in the Alps and the Pyreneesis different. There are sedimentary basins north and south of the Pyrenees but theyare not comparable in depth, extent and velocity to the marine situations. For thespecific situation of the western Pyrenean crust Chazalon et al. (1993) modeled thewave propagation with a realistic velocity model including the Moho jump and highvelocity obstacles in the upper crust. Chazalon et al. (1993) conclude from theirmodeling that neither the Moho jump nor the detailed crustal velocity structure cancause the observed attenuation. The inability of the Moho jump to cause the Lg-blockage is supported by the observed Lg propagation through the eastern Pyreneeswhere the Moho jump is also present.

It can be summarized that that the observations reported in figure 3.2 reflectintrinsic properties of the propagation medium and not of the sources or the velocitystructure. There is no satisfactory explanation for the Lg-blockage by the Pyreneesor the Alps. It was speculated by Chazalon et al. (1993) that the increased scatteringby small scale heterogeneities in the western Pyrenean crust causes the attenuationof crustal phases.

3.4 Modeling

To test this hypothesis, the energy propagation through the Pyrenees is modeledwith elastic radiative transfer theory (RTT). This theory in the acoustic approxi-mation was applied in a number of studies focusing on the separation of intrinsicand scattering attenuation in simple half space models (e.g Hoshiba, 1991; Fehleret al., 1992; Gusev and Abubakirov, 1996). Lacombe et al. (2003) use a model thatconsists of a layer overlaying an infinite half space for this purpose. Here the elasticversion of RTT is used in a complex geometrical model that consists of a scatteringlayer over a scattering half space with a vertical velocity gradient. Additionallya cuboid body with different scattering and attenuation properties is included inthe crustal layer to model the local effect of the Pyrenean crust. Because of thiscomplex model the main focus is not on the estimation of precise parameters of themodel. Instead I try to present a conceptual model to explain the phenomenon ofLg-blockage in the Pyrenees.

The velocity model that is used for the simulation is based on the model of theLDG. The crust is a constant velocity layer with density ρcrust = 3.0 g/cm3, P-wavevelocity αcrust = 6.03 km/s and S-wave velocity βcrust = 3.56 km/s. The mantleis represented by a half space with constant vertical velocity gradient and constantdensity ρmantle = 3.5 g/cm3. At the Moho the P-wave velocity is αmantle = 8.16 km/s

INVERSION 61

and S-wave speed βmantle = 4.65 km/s. The velocity gradient for P-waves is 1300

km/skm

which is a approximate average of the velocity profile in the uppermost 600 km man-tle in the Preliminary Reference Earth Model (PREM) (Dziewonski and Anderson,1981). The velocity ratio γ is constant throughout the mantle. Small scale hetero-geneities are described by an exponential autocorrelation function (equation 2.27)with the parameters fluctuation strength ε and correlation length a. The scatteringproperties in crust and mantle are different but constant, i.e. the scattering proper-ties do not vary with depth in the mantle. Intrinsic attenuation is described by thequality factor. I fix the ratio between the quality factor of P- and S-waves to the theo-retical value that assumes that attenuation is limited to shearing: Q = QP = 9/4QS

(Shearer, 1999, p. 114).

In this model I include an additional cuboid body in the crust that has thesame velocities and density as the surrounding crust but differs in the scatteringand attenuation properties. The body extends from the surface down to the Mohoand has a width in propagation direction of 100 km. This width corresponds tothe topographic expression of the Pyrenees. In the following I try to show whetheror not the attenuation and scattering properties of this body can account for theLg-wave blockage.

Since the reference traces are the result of averaging over several differentsources I use sources with isotropic radiation pattern. According to equation 2.19the ratio between radiated P and S-energy is fixed at 96% S-energy and 4% P-energy.

3.5 Inversion

The goal of this non-linear inversion process is to find scattering and attenuationparameters for crust, mantle and Pyrenean body that result in energy density traceswhich fit the eastern and western reference traces. Though some of the parametersmay be estimated independently the whole model has 9 free parameters, the fluctu-ation strength ε the correlation length a and the quality factor Q in three differentbodies. Here the type of the medium is already fixed as I decided to restrict theinversion to the exponential auto-correlation function. Because of the high numberof free parameters I apply a genetic algorithm for the inversion. The principal el-ements of the algorithm are recombination of the parameters of the fittest modelsand mutation to maintain diversity. Fitness of the models is linked to the misfitbetween the model and the reference trace. The misfit is the sample wise averageof the absolute value of the difference between the logarithm of the model envelopeand the logarithm of the reference. This is the L1-norm of the difference betweenthe logarithmic traces. It ensures that the low amplitude coda has similar weight forthe misfit calculation as the high amplitude direct arrivals and is not too sensitiveto fluctuations in the data and the simulated envelopes.


3.5.1 The undisturbed case

In a first inversion I estimate the parameters of the background model. For this, thecrustal phases of the eastern reference trace are fitted by adjusting the parametersa, ε and Q of the crust in a model that has a transparent (no scattering) andnon-attenuating mantle and no additional body. The lapse time windows used tocalculate the misfit in this inversion step are the Pg window between the Pg arrivaland the Sn arrival and the S-coda from the Sg arrival until a group velocity of2 km/s. Similar to the generation of the reference traces I stack envelopes in acertain distance range after scaling to equal group velocity. For comparison withthe eastern reference trace I stack records between 500 km and 800 km distancefrom the source.

With this model of the crust I invert the western reference trace for scatteringand attenuation properties of the mantle. Because of the blockage of the crustalphases the wave paths through the western Pyrenees offer the interesting possibilityto estimate mantle properties relatively independent of the crust. In this inversionstep I use time windows that contain the mantle phases, i.e. the Pn window be-tween the Pn and the Pg arrivals and the Sn window between the Sn and the Sgarrivals. Records between 600 km and 900 km are stacked in this case to accountfor the slightly longer wave paths through the western Pyrenees. The results of thisinversion step are not well determined and a relatively wide range of values for theparameters ε, a and Q can be used to fit the reference trace. I therefore fixed the pa-rameter for the fluctuation strength and correlation length of the mantle to ε = 0.02and a = 2 km. This resulted in the quality factors QP = 1070 and QS = 475.

With these new parameters for the mantle I re-inverted the eastern refer-ence trace for the crustal properties because the increased scattering in the mantlechanged the amount of energy that returns from the mantle into the crust. Thisprocess raised the level of the coda and required changes of the parameters in thecrust. The final parameters of the crustal material are ε = 0.021, a = 0.77 km,QP = 1400 and QS = 623.

Figure 3.6 shows the result of this inversion. The general fit of the model isfairly good. Especially the decay of the S-coda is well fitted. The kink in the S-codaenvelope shortly after 3 km/s originates from the reverberations in the crust. Itmarks the end of the Lg-wave train which is well modeled by the algorithm. Thefit is less good for the P-coda and the Pn-phase. The deviation from the shape ofthe Pn-arrival is due to the fact that I estimated the mantle properties from thewestern reference trace. I can better fit the Pn-phase with stronger scattering in themantle but then the impulsive shape of the mantle phases in the western referencetrace cannot be reproduced. In fact this is an indication for slightly different mantleproperties under the eastern and western Pyrenees.

3.5.2 The disturbed case

To model the Lg-blockage I use the background model described in the previoussection and include an extra body with different scattering and attenuation prop-

INVERSION 63

8.16 6.03 4.65 3.56 3 2v [km/s]

Figure 3.6: Simulation Results for the undisturbed propagation in through the easternPyrenees. Bold red and black traces represent data and model respectively. Thin tracesindicate types of energy. Cyan trace shows radiated P-energy that arrives at the receiveras P-energy. The difference between the cyan and blue traces shows radiated P-energythat arrives as S-energy. The difference between blue and green traces shows radiated S-energy that arrives as P. Radiated S-energy that arrives as S is indicated by the differencebetween the green and the black trace. Vertical lines indicate arrival times of seismicphases

erties. To speed up the Monte-Carlo simulations I slightly modify this model toenforce cylindrical symmetry which allows better averaging. In the modified modelthe cuboid attenuating body is replaced by a cylinder ring centered at the source.The radius of the inner boundary is 300 km and the width is 100 km. I verified thatthe effect of this modification on the inversion is very weak.

Three different variants of the model are used for the inversion. In the firstvariant the fluctuation strength of the small scale heterogeneities is adjusted to fitthe western reference trace. The intrinsic attenuation and the correlation lengthin the inclusion are equal to the surrounding crust. With ε = 11% the best fit isobtained (figure 3.7A). The general shape of the data is reasonably well fitted. Thestrong scattering in the inclusion efficiently attenuates the crustal phases. There isno Pg-arrival and the Sg-phase is only visible as a small bump in the Sn-coda. Thisshows that scattering is a good candidate to cause the Lg-blockage in the Pyrenees.But there is some mismatch between this model and the data in the S-coda whichhas too much energy in this part of the envelope.

In a second variant of the model I assume that the Lg-blockage is caused byintrinsic attenuation. The inclusion has the same scattering properties as the normalcrust but variable Q values. Again the ration between QP and QS is fixed at 9/4.The best fit is obtained for QP = 86 and QS = 38. These values are extremely low.Figure 3.7B shows the results of the model. The fit is worse compared to the fitachieved by changing the scattering strength. The S-coda of the model is much toolow but there is still Pg-energy that traverses the inclusion. It is impossible to bothattenuate the Pg-phase strong enough but simultaneously increase the energy in the


A8.16 6.03 4.65 3.56 3 2

v [km/s]B

8.16 6.03 4.65 3.56 3 2v [km/s]

Figure 3.7: Simulation results for propagation through the western Pyrenees withstronger scattering (A) and increased attenuation (B). For explanations see figure 3.6.

8.16 6.03 4.65 3.56 3 2v [km/s]

Figure 3.8: Simulation results for propagation through the western Pyrenees with bothstronger scattering increased attenuation. For explanations see figure 3.6.

S-coda. Leaving the QP/QS-ratio free would allow to obtain a much better fit. Butthis would require much stronger P- than S-attenuation which is unrealistic. Thisshows the attenuation alone does not suffice to explain the Lg-blockage.

In the third model I allow for variable scattering and attenuation in the in-clusion. With this model the best fit is obtained (figure 3.8). The parameters areQP = 402, QS = 179 and ε = 7.7%. With these parameters the crustal phasescan both be appropriately attenuated. The Pg-phase is weak enough not to thestick out of the Pn-coda and simultaneously the small signature of the Lg-wave isreproduced. Also the amplitude ratio between the Pn and Sn phases is reasonable.For this model the parameters are listed in table 3.1.

I conclude that strong attenuation alone cannot explain the observation of Lg-blockage in the Pyrenees. Increased scattering seems an essential ingredient of thisphenomenon but is not able to explain all details of the observation. In contrast theobservation can be well explained with increased scattering and stronger attenuation

DISCUSSION 65

Table 3.1: Parameters of the final model used to explain the phenomenon of Lg-blockageby the Pyrenees. The columns show values for the quality factor of P-waves (QP ) andS-waves (QS), the fluctuation strength (ε) and the correlation length (a) of the small scalevelocity fluctuations and the transport mean free path l∗S (Turner, 1998). The values ofε and a in the mantle are fixed by hand. The correlation length a of the Pyrenean bodyis assumed to equal the correlation length of the surrounding crustal material. Values inbrackets show representative numbers of the PREM (Dziewonski and Anderson, 1981).The 750 km transport mean free path can be compared to the mean free path of 690 kmobtained in chapter 1.5 with the isotropic scattering model for the German crust.

QP QS ε a l∗SCrust 1400 (1350) 623 (600) 2.1% 0.77 km 750 km

Mantle 1070 (321) 475 (143) 2.0% 2.0 km 1530 kmPyrenean body 402 179 7.2% 0.77 km 63 km

together.

3.6 Discussion

Compared to the waveform modeling by Chazalon et al. (1993), I use a simple struc-tural model that consists of a flat layer over an infinite half space. It does neithercontain any deterministic velocity structure inside the crust nor the Moho topogra-phy. Despite these limitations the shape of the seismogram envelopes in the easternPyrenees can be reproduced. Including an extra body with stronger scattering andattenuation the energy propagation can be influenced to reproduce the seismogramenvelopes that are observed in the western Pyrenees. For the final parameters Iillustrate the energy field with two snapshots of the simulations in figures 3.9A andB. These figures show the shadowing effect of the inclusion (marked by the blackbox) in the western Pyrenees for crustal waves. The epicenter is indicated by thestar. The linear color-scale represents the square root of seismogram envelopes, i.e.the fourth root of energy. Figure 3.9A displays the energy distribution at a lapsetime of 105 s. At this time the P-phases already passed the Pyrenees. The Pnmantle phase which propagated north to about 47 is unaffected by the Pyrenees.Its energy is low compared to the other phases and the speckle reflect the individualparticles of the Monte-Carlo simulations. P-energy that propagated in the crustforms the Pg-phase bounded towards north at about 45. This phase experiencedstrong attenuation in the western Pyrenean crust. The Sn wave front just leaves thebox and the Lg-wave enters it. Figure 3.9B shows the situation at 180 s lapse time.The P wave fronts already left the map and the Sn and Lg wave front extend northto about 47 and 45 respectively. Again the mantle phase is continuous behind thePyrenees whereas the crustal Lg shows a gap. The almost homogeneous color be-hind the Lg-wave illustrates the spatial equalization of energy density in the seismiccoda.

The simulations show a consistent picture of the phenomenon of Lg-blockage.


A

-6˚

-6˚

-5˚

-5˚

-4˚

-4˚

-3˚

-3˚

-2˚

-2˚

-1˚

-1˚

0˚

0˚

1˚

1˚

2˚

2˚

3˚

3˚

4˚

4˚

5˚

5˚

6˚

6˚

37˚ 37˚

38˚ 38˚

39˚ 39˚

40˚ 40˚

41˚ 41˚

42˚ 42˚

43˚ 43˚

44˚ 44˚

45˚ 45˚

46˚ 46˚

47˚ 47˚

48˚ 48˚

49˚ 49˚

200 km

low heigh

amplitudeB

-6˚

-6˚

-5˚

-5˚

-4˚

-4˚

-3˚

-3˚

-2˚

-2˚

-1˚

-1˚

0˚

0˚

1˚

1˚

2˚

2˚

3˚

3˚

4˚

4˚

5˚

5˚

6˚

6˚

37˚ 37˚

38˚ 38˚

39˚ 39˚

40˚ 40˚

41˚ 41˚

42˚ 42˚

43˚ 43˚

44˚ 44˚

45˚ 45˚

46˚ 46˚

47˚ 47˚

48˚ 48˚

49˚ 49˚

200 km

low heigh

amplitude

Figure 3.9: Snapshots of the energy density field at 105 s (A) and 180 s (B) lapsetime for an earthquake in Spain. Green star and black box indicate the epicenter andthe disturbed body beneath the western Pyrenees. respectively. The linear color-scaleindicates the square root of seismogram envelopes, i.e. fourth root of energy density.Relief is indicated by brightness. The shadowing effect of the medium in the box can beseen for the Pg-wave in part A and for the Lg-wave in part B.

They are obtained with model parameters (see table 3.1) that appear to be reason-able. The parameters for the crust match existing estimates quite well. The qualityfactor of intrinsic S-wave attenuation (623) is close to the value of 600 given in thePREM (Dziewonski and Anderson, 1981). It is also comparable to the estimateobtained with the isotropic half space model in chapter 1 of 770. The scatteringproperties are best described by the S-wave transport mean free path of 750 km.Again this is close to the 690 km found in chapter 1. The parameters for the mantleare not well constrained. In fact the scattering parameters are fixed manually butcertainly the transport mean free path has to be larger than in the crust to fit theimpulsive Pn-arrival of the western reference trace. The medium that causes the isdescribed by a transport mean free path of 63 km which is in the range of parame-ters estimated in other regions of the world. The attenuation with a S-wave qualityfactor of 179 is rather strong for crustal material (c.f. table 1.2). But this is incorrespondence with implication of the seismicity distribution. The intrusive bodiesof lower crustal material of which two were identified in the seismically active uppercrust of the wester Pyrenees are seismically quiet (Souriau et al., 2001). Souriauet al. (2001) proposed that this is due to the more ductile lower crust material whichmost likely shows stronger attenuation.

However, no confidence intervals are given for the model parameters and I do

CONCLUSIONS 67

not put any emphasis on these values because their exact determination is not thescope of this study. Also the geometrical parameters of the model like the extendof the perturbed zone, the depth of the Moho or the velocity gradient in the mantleshould not be regarded precise estimates. It is assumed that the mantle phases Pnand Sn are made up of energy that left the crust through the Moho and returnedafter traveling through the mantle according to ray-theory. Energy that travels asinterface wave along the discontinuity is neglected just like surface waves. It mightthus be possible to find models with different parameters which similarly fit thedata. However, these effects are small compared to the Lg-blockage and will thusnot influence the general results and conclusions of this study. I present a conceptualmodel with realistic parameters that explains the so far unexplained phenomenonof Lg-blockage in the western Pyrenees.

3.7 Conclusions

With the elastic Monte-Carlo simulations in a model with deterministic velocitystructure and statistically described small scale heterogeneities it is possible to cal-culate seismogram envelopes of crustal earthquakes including mantle phases andcoda. Including an additional body with different scattering and attenuation prop-erties this model provides a possible explanation for the Lg-blockage by the Pyrenees.To coincide with the geographical observation this body should be located under thewestern Pyrenees. This region is characterized by strong deformation. The crustwas squeezed during the compressional phase that led to the formation of the Pyre-nees. There is geological and geophysical evidence for material exchange betweenthe different depths. Mantle material is lifted to lower crustal depths and lowercrust material can be found in the upper crust. These observations appear to beconsistent with increased heterogeneity as suggested by my model. In the easternand central Pyrenees the shortening led the the subduction of the Iberian lower crustbeneath the European plate. This process reduced the deformation in the crust.

I speculate that the strong deformation in the western Pyrenees caused somekind of mixing of mantle and crustal material that increased the heterogeneity andthe attenuation in the western Pyrenean crust. These material changes cause theattenuation of crustal phases known as Lg-blockage. The model might also beapplicable to the situation in the western Alps (Campillo et al., 1993).

Part II

Phase Correlations

69

Chapter 4

Passive image interferometry andseasonal variations of seismicvelocities at Merapi volcano,Indonesia

71

72 CHAPTER 4. PASSIVE IMAGE INTERFEROMETRY AT MT. MERAPI

Abstract

Passive image interferometry is proposed as a novel technique for seismology thatallows to continuously monitor small temporal changes of seismic velocities in thesubsurface. The technique is independent of sources in the classical sense and re-quires just one or two permanent seismic stations. This is possible because thesignals that are used for interferometry are retrieved from ambient seismic noise.Applying passive image interferometry to data from Merapi volcano, it is shownthat velocity variations can be measured with an accuracy of 0.1% with a temporalresolution of a single day. At Mt. Merapi the velocity variations show a strong sea-sonal influence and a depth dependent hydrological model is presented that describesthe observations solely based on precipitation.

4.1. INTRODUCTION 73

4.1 Introduction

Imaging techniques with elastic waves like seismic tomography, reflection seismicor ultrasonic imaging proved their great abilities for the determination of spatialstructure. Numerous applications range from medicine to hydrocarbon exploration.Temporal changes are often thought to be of little importance and repeated mea-surements designed for structural investigations are used to detect such changes.However, in some fields of geophysics and nondestructive testing, monitoring tem-poral changes is of greater interest than details of the spatial structure. Especiallyfor monitoring of volcanoes, fault zones, dams, and hydro-carbon or geothermalreservoirs it is valuable to observe changes of elastic properties like seismic velocity,even if the precise spatial structure is unknown. A technique named coda waveinterferometry (CWI) has been proposed (Snieder et al., 2002; Snieder, 2006) formeasuring weak changes of seismic velocities using multiple scattered waves. Thesame approach was previously used by Poupinet et al. (1984) and Ratdomopurboand Poupinet (1995). The basic concept of CWI is that despite their complexitythe waveform of the seismic coda is constant for constant source receiver configu-rations. On the other hand the delay time, i.e. the travel time variation due toa velocity variations in the medium increases with travel time. So the changes inthe medium are inferred by comparison of the coda obtained from the same sourcereceiver configuration at different times. This way CWI strongly depends on thephase information of the seismic coda which is neglected in the radiative transfertheory applied in part I.

CWI has been applied to measure velocity changes associated with heat, pres-sure, and water saturation in a rock sample (Gret et al., 2006), with temperaturechanges in a building (Larose et al., 2006) and with earthquakes (Nishimura et al.,2000; Peng and Ben-Zion, 2006). Velocity changes were also reported to precedeeruptions of Merapi volcano (Ratdomopurbo and Poupinet, 1995; Wegler et al.,2006b). The major handicap of CWI, which hinders continuous monitoring in thefield, is the requirement of repeatable sources. This results in temporally irregularmeasurements when repeating earthquakes (multiplets) are used (Ratdomopurboand Poupinet, 1995; Peng and Ben-Zion, 2006) or necessitates repeatable activesources that are expensive on shore (Wegler et al., 2006b; Nishimura et al., 2000).

Here I suggest to use a technique called passive imaging to obtain the signalsfor CWI. This technique requires only one or two permanent stations and workswithout active sources. The basic principle of passive imaging is that the Green’sfunction (GF, or impulse response) between two points A and B can be retrievedby cross-correlation of a random isotropic wave field sensed at A and B (Derodeet al., 2003; Sanchez-Sesma and Campillo, 2006; Wapenaar, 2004; Snieder, 2004).In real applications the cross-correlation of a random wave field will only yield anapproximation of the real GF. Consequently a term like field correlation functionappears more appropriate. However, to stress the relation to the impulse response,that the field correlations functions converge to for adequate wave fields or sufficientaveraging, the term Green’s function is used here.

Applications of passive imaging were reported from helioseismology (Duvall


110.4˚ 110.45˚ 110.5˚

−7.6˚

−7.55˚

−7.5˚

3 km

110.4˚ 110.45˚ 110.5˚−7.6˚

−7.55˚

−7.5˚

3 kmBEB

BAT

KEN

GRW

1000

15002000

0.1 km

GRW

0

1

23

Figure 4.1: Location of seismic stations on Mt. Merapi. Triangles denote seismometersand stars indicate sources of active CWI measurements (Wegler et al., 2006b).

et al., 1993; Rickett and Clearbout, 1999) and from acoustics (Lobkis and Weaver,2001; Weaver and Lobkis, 2001). In seismology GFs were retrieved from seismicnoise and from coda waves of earthquakes (Campillo and Paul, 2003; Shapiro andCampillo, 2004; Paul et al., 2005; Larose et al., 2005; Roux et al., 2005; Shapiroet al., 2005). So far the recovery of the GF with passive imaging has only beendemonstrated for ballistic surface waves and P-waves (Roux et al., 2005). Evidencefor the emergence of scattered waves from the GFs in seismology has been missingso far. In this study it is shown that scattered waves can be retrieved by correlationof seismic noise.

GFs are studied at Mt. Merapi (figure 4.1) which is one of the most activestrato volcanoes on earth. Its activity is characterized by dome growth and partialdome collapse threatening the surroundings with pyroclastic flows. Most recentactivity started in middle of April 2006 and caused evacuations of several thousandinhabitants in May and June. Monitoring the state of Mt. Merapi is of greatsocio-economic importance. Since 1997 three seismic arrays have been operated onMt. Merapi (Wassermann and Ohrnberger, 2001). I focus on data from stationsGRW0 and GRW1 (figure 4.1) because these sensors were connected to a commondigitizer which provides precise relative timing. Data are almost continuous betweenAugust 1997 and June 1999. This study is published in Sens-Schonfelder and Wegler(2006a).

4.2 Green’s function retrieval

Green’s functions are retrieved by cross-correlating two one-day seismic records thatare high passed at 0.5 Hz. To down-weight the contribution of coherent phases suchas teleseismic body waves in the correlation, the records are clipped at one standarddeviation of the recorded seismic noise. Additionally, the causal and acausal parts ofthe GFs which represent waves that travel from A to B and vice versa, are averaged.Since I use three-component receivers I can reconstruct the full Green’s tensor (GT,

4.3. SCATTERED WAVES IN THE GREEN’S FUNCTION 75

Paul et al. (2005)) consisting of nine GFs, by correlating all combinations of thecomponents.

4.3 Scattered Waves in the Green’s Function

In order to apply CWI on the GFs retrieved with passive imaging I first establishthat multiply scattered or reflected waves are contained in the GF. I present twoarguments. Firstly, GFs retrieved from noise at different days show coherent phasesat late lapse times denoted τ (figure 4.2A). Given a distance of about 170 m betweenstations GRW0 and GRW1 the late arrivals (e.g. the one marked by the arrow infigure 4.2A) clearly correspond to indirect wave paths, i.e. reflected or scatteredwaves. Coherent phases can also be observed in the auto-correlation function that isan approximation of the zero offset GF (figure 4.3). Secondly, I study the envelopesof the GT in rotated coordinates where L is oriented along the connecting line ofthe stations, Q is perpendicular to the free surface and T is perpendicular to Land Q. Figure 4.4 shows envelopes of the Green’s tensor components in which theoff-diagonal components are averaged and the traces have been normalized in thetime window between 10 s≤ τ ≤ 12 s. Two different segments can be observed. Theearly part (τ < 2 s) contains surface waves that travel in a low velocity layer at agroup velocity of about 200 m/s. After two seconds lapse time all GT componentsexhibit a common decay. This shape is reminiscent of the typical decay of codawave envelopes from active sources. The decay can be fitted with a diffusion model(Dainty and Toksoz, 1981) expressed by: W (τ) ∝ τ−3/2 exp[−bτ ]. Here W is theseismogram envelope, τ denotes lapse time. b is the absorption parameter. Theestimates of parameter b using GFs reconstructed from seismic noise agree wellwith measurements on Mt. Merapi performed with active sources (Wegler andLuhr, 2001). This supports the hypothesis that the late part of the GFs consists ofscattered coda waves.

4.4 Passive image interferometry

I have presented evidence for the fact that the GFs obtained from cross-correlationsof seismic noise do not only contain the ballistic wave part of the GF as previ-ously assumed in the seismological literature (Campillo and Paul, 2003; Shapiro andCampillo, 2004; Shapiro et al., 2005; Paul et al., 2005; Roux et al., 2005) but alsoreflected and scattered waves. This opens the possibility to apply interferometricmethods to the late part of the passively retrieved Green’s functions (images). Icall this technique passive image interferometry. It allows to precisely infer changesin the medium by comparing the GFs retrieved from the noise records at differenttime periods.

Standard CWI (Snieder et al., 2002) measures the delay times (δτi) in varioustime windows i centered around lapse times τi as the time shifts that result inthe best correlation of the segments in these windows. Alternatively the delay times(δτi) can be estimated from the derivative of the of difference of the phase spectra in


0 2 4 6 8 1005/27/9905/25/9905/23/9905/21/9905/19/9905/17/9905/15/9905/13/9905/11/9905/09/9905/07/9905/05/99

τ [s]7.1 7.15 7.2 7.25 7.3

τ [s]A B

Figure 4.2: Coherent phases in the Green’s function coda obtained from cross-correlation.A: cross-correlation between Z-components from stations GRW0 and GRW1 for severaldays. Note the coherent phase marked by the arrow. B: Closeup around the phase markedby the arrow in A. This phase which is seen in all GFs arrives earlier in late May 1999than in early May. This change in the delay time corresponds exactly to the apparentvelocity change of 0.015% per day as seen in figure 4.6A.

0 2 4 6 8 1005/28/9905/26/9905/24/9905/22/9905/20/9905/18/9905/16/9905/14/9905/12/9905/10/9905/08/9905/06/99

t [s]3.85 3.9 3.95

t [s]A B

Figure 4.3: Coherent phases in the Green’s function coda obtained from auto-correlation.Similar to figure 4.2. A: auto-correlation of the Z-component of station GRW1. Note thecoherent phase marked by the arrow. B: Closeup around the phase marked by the arrowin A.

these time windows with respect to frequency. This approach is used by the MovingWindow Cross Spectral Analysis (MWCSA) that was applied by Ratdomopurboand Poupinet (1995). The relative delay time (δτ/τ) is then the mean of all δτi/τivalues.

In contrast to CWI and MWCSA the relative delay time is determined here asthe factor by which the time axis of one trace has to be stretched or compressedto obtain the best correlation with the other trace. This technique is illustratedin figure 4.5 where the relative delay time is to be measured between the blue andthe red trace shown in the upper part. The red trace is successively stretched byincrements of 1%. The correlation coefficient between the blue and the stretched redtrace is depicted in the lower part of figure 4.5 for each increment of the stretching,which is directly related to the relative delay time δτ/τ . The maximum correlation

4.4. PASSIVE IMAGE INTERFEROMETRY 77

0 5 10 15

100

102

τ [s]

Figure 4.4: Envelopes of the GT components between stations GRW0 and GRW1 inrotated coordinates (see text for orientation of coordinates). Non-diagonal components ofthe GT are averaged with their counterparts. Dashed lines show LL, QQ, TT and LQcomponents of the GT which can be excited by ballistic waves. Thin lines show the LTand QT components which can not be excited by ballistic waves. The best fitting diffusionmodel is shown by a bold line.

0 0.02 0.04 0.06 0.08−0.5

0

0.5

1

corr

elat

ion

δτ / τ

Figure 4.5: Illustration of the technique that is used to estimate the relative delay timebetween two traces. Upper graph: the two traces that are to be compared in red andblue. Black trace is a copy of the red trace stretched by 8%. Lower graph: correlationcoefficient between the blue and the stretched red trace for various amounts of stretching.

is assumed for a stretching of 5%, corresponding to a relative delay time of 0.05.

The advantage of this approach becomes evident by comparison with standardCWI. Both techniques measure the similarly between two traces by means of thecorrelation coefficient. But standard CWI estimates the correlation in a number ofsmall time windows. The size of these windows has to be chosen such that theyare long enough to calculate a meaningful correlation coefficient but they have tobe short enough to ensure that the distortion caused by the medium change issmall within these windows. Also the number of these time windows has to be largeenough for a reliable least squares fit of the δτi−τi curve. In the approach presentedhere there is no need to adjust such algorithmic parameters. Since the distortionalong the trace is accounted for by the stretching or compression, the similarity


between the two traces can be evaluated along the whole trace. This means thatthe correlation coefficient can be calculated in a single time window ranging fromthe beginning of the trace until the signal reaches the noise level. This makes ourapproach more robust and sensitive.

If the time delay is caused by a spatially homogeneous relative velocity changeδv/v, the relative delay time, is independent of the lapse time at which it is measuredand δv/v = −δτ/τ . In this case error estimates are obtained by choosing consecutivenon-overlapping time windows in which the correlations are calculated separately.This gives independent estimates of the relative delay time and mean value andstandard deviation can be calculated.

4.5 Application to Merapi volcano

Passive image interferometry is applied to almost two years of nearly continuousseismic records from Merapi volcano, Indonesia. Velocity changes are estimated foreach one-day GF (between stations GRW0 and GRW1) with respect to a referencetrace, assuming that δv/v is spatially homogeneous. As reference trace a stack of allone-day GF in January 1998 is used. Figure 4.6A shows the daily relative velocityvariation obtained from the LL component of the GT between 2 and 8 s lapse time.I do not use lapse times smaller than 2 s because the surface waves contained inthe GFs prior to 2 s have a different spatial sensitivity to velocity variations. Grayshading in figure 4.6A marks the standard deviation of independent measurementsin three consecutive 2 s time windows starting at 2, 4, and 6 seconds lapse time.

Seismic velocity shows temporal variations of a few percent with a clear seasonaltrend. Except for a constant offset, the curve only marginally depends on the periodused to generate the reference trace. However stacking the GFs from a longer periodof time may result in a poorer correlation of the reference trace with the dailycorrelations because the velocity variations that occurred in the longer period woulddegrade the reference trace. January 1998 is chosen because the velocity appearedto be stable during this period.

Evidence for the fact that the inferred velocity variations are related to phys-ical properties of the subsurface and not artefacts introduced by the sensors or therecording system is provided by the analysis of records from different componentsand stations. In figure 4.7A the velocity changes measured from the transverse com-ponents are shown. The results are virtually identical to the ones obtained from theLL-component. An interesting property of our new method for practical applica-tions in geophysics is that it can be applied to records from a single station usingthe auto-correlation function. Figure 4.7B shows the velocity variations measuredon the auto-correlation of the vertical component of station GRW0. Again the curveis very similar to the one in figure 4.6A.

The blue curve in figure 4.6A shows the results of CWI measurements (Wegleret al., 2006b) made at the same stations (GRW) but with an active source (BATin figure 4.1). Results from the southern flank of Merapi volcano (station KEN,source BEB) are shown by the green curve. Considering the daily variability of the

4.6. HYDROLOGICAL MODELING 79

measurements and the source receiver distance of more than 3 km used for the activemeasurement, I conclude that the present measurements are in agreement with theresults of the active experiment. The advantage of my new method is that changescan be monitored continuously without an expensive active source.

Actually the velocity changes can directly be seen in the GFs. Figure 4.2B and4.3B show close-ups around the phases marked by the arrows in figures 4.2A and4.3A respectively. The coherent phases tend to arrive earlier in late May than inearly May 1999. This corresponds to the almost linear velocity increase of about0.015% per day in May 1999 that can be observed in figure 4.6A.

4.6 Depth dependence of the velocity changes and

hydrological modeling

The periodicity of the velocity fluctuations of approximately one year in this tropicalenvironment suggests a climatic influence, most likely by precipitation. Anotherobservation also supports a connection with precipitation: The amplitude of thevelocity variations depends on the lapse time at which it is measured. This hasnever been reported for CWI measurements before and indicates that the velocityperturbations are spatially inhomogeneous. I can explain this observation togetherwith details of the temporal variations with a hydrological model of the groundwater level (GWL).

An illustration of the hydrological model and its influence on the seismic veloc-ity structure is shown in figure 4.8. I assume that drainage of ground water occursthrough a stationary aquifer that can approximately be described by Darcy’s law.Then the drainage is proportional to the height of the ground water table whichresults in an exponential decrease of the water level after rain events. A convolutionof the precipitation rates with this exponential function thus gives the ground waterlevel GWL (below surface) at time ti:

GWL(ti) = GWL0 −i∑

n=0

p(tn)

φe−a(ti−tn). (4.1)

Here φ is porosity, a is the parameter describing the decay, GWL0 is the asymptoticwater level, and p(ti) denotes the daily precipitation. Good agreement of such amodel was found with water level measurements in a well (Akasaka and Nakanishi,2000). Daily precipitation data were provided by the NASA/Goddard Space FlightCenter’s Laboratory for Atmospheres (Huffman et al., 2001).

To relate the ground water level with delay times, I define the depth (z) andtime (t) dependent relative slowness perturbation S(ti, z) and a reference water levelGWLref that is chosen equal to the mean level of January 1998. Then S(ti, z) = δsfor GWL(ti) < z < GWLref , S(ti, z) = −δs for GWLref < z < GWL(ti) andS(ti, z) = 0 elsewhere. δs is the relative slowness difference between the states


saturated and dry. The delay time δτ at time ti measured at lapse time τ is then

δτ(ti, τ) =

0∫

z=−∞

S(ti, z)K(z, τ)dz (4.2)

where K(z, τ) is the sensitivity kernel (Pacheco and Snieder, 2005). I use the kernelin the approximation of coincident source and receiver with a diffusivity constant ofseismic energy of D = 0.05 km/s2 as estimated at Mt. Merapi (Wegler and Luhr,2001). The kernel decays with distance r from the station as r−1 exp(−r2/Dτ).For the one dimensional structure of the velocity perturbation resulting from ourhydrological model, figure 4.9 shows the exact depth dependent sensitivity kernelsfor lapse times τ = 3s and τ = 7s.

With a genetic algorithm I found the best model with the following parameters:GWL0 = 51 m, φ = 0.03, a = 0.008 d−1, and δs = 0.17. Figure 4.6B shows the dailyrainfall, the inferred GWL and the modeled apparent homogeneous velocity varia-tion for the 2–4 s and 6–8 s lapse time windows together with the measured values.Considering that the precipitation data are averages over an area of approximately100 by 100 km and the simplicity of our hydrological model the velocity variationscan be explained in remarkable detail. The model also explains the difference be-tween the blue and the red curves in figure 4.6B that represent the measurements atdifferent lapse times. I am thus confident that the observed velocity variations arerelated to precipitation via changes in the ground water level. The successful expla-nation of the lapse time dependence of the relative delay times by depth dependentvelocity perturbations indicates that the coda in the frequency range above 0.5 Hzis made up of scattered body waves rather than surface waves. In contrast to bodywaves the depth profile of the sensitivity kernel of surface waves does not changewith lapse time.

For spatially inhomogeneous velocity perturbations, the accuracy of the methodcan be characterized by the scatter of measurements around a common trend. Inrain-free periods like the first three month in the study period I obtain a standarddeviation of 0.1%.

Assuming that the annual rainfall is similar each year the observation canbe compared with measurements from 1991 (Ratdomopurbo and Poupinet, 1995).Indeed a velocity increase of a little more than 1% was observed during the periodfrom April until September 1991. During April to October 1998 and April to June1999 we observe a similar increase. According to my interpretation these periodsof increasing velocities are caused by a falling ground water level in the dry season.I conclude that the effect of stress changes inside the volcano on velocity changes(Ratdomopurbo and Poupinet, 1995; Wegler et al., 2006b) is secondary.

4.7 Conclusions

The present study shows that the hydrological conditions can change the seismicvelocities by more than 10% in a depth range of several tens of meters. This effect

4.7. CONCLUSIONS 81

should be taken into account in studies of the local velocity structure. With thehydrological model an independent check for PII is provided. It proves that de-terministic information about the propagation medium is contained in the GFs atlarge lapse times, i.e. one can retrieve scattered and reflected waves by correlation ofseismic noise. This provides the basis for passive seismic imaging with non-ballisticwaves. I anticipate this observation to stimulate a variety of studies in reflectionseismology using correlations of seismic noise.

The temporal velocity variations detected at Mt. Merapi in this study aredominated by shallow perturbations and demonstrate the potential of PII. If thescattering is weaker or the changes close to the stations are smaller, the sensitivitykernel allows the detection of more distant changes. In a lower frequency range theGF are dominated by surface waves. The depth sensitivity of surface waves is verydifferent from body waves and opens new possibilities for passive image interfer-ometry. The is a variety of possible applications for PII (Snieder and Page, 2007).Monitoring of the exploitation process of geothermal or hydrocarbon reservoirs withPII seems possible as well as the monitoring of nuclear waste disposal sites. I thinkthat PII can be used for routine monitoring in various engineering, geotechnical,and geophysical applications.

In a separate study a drop of seismic velocity associated with the Mid-Niigataearthquake in Japan is reported that was detected with PII (Wegler and Sens-Schonfelder, 2007). Crucial for the results presented here is the precise relativetiming of the two sensors. But it is shown that similar results can be obtained fromjust a single station and in chapter 5.2 introduces a technique to measure clockerrors of seismic stations. This means that PII can be applied to data from anyexisting permanent seismograph station and we expect other interesting discoveriesfrom the application of PII in different environments.


−0.04

−0.03

−0.02

−0.01

0

0.01

0.02

0.03

0.04

dv/v Jun Jul Aug

−0.01

−0.005

0

Sep Nov 1998 Mar May Jul Sep Nov 1999 Mar May Jul

40

30

20

10

0

−0.02

−0.01

0

0.01

0.02

GW

L [m

]

δv

/v

0

30

60

prec

ipita

tion

[mm

/d]

A

B

Figure 4.6: Measured and modeled velocity variations at Merapi volcano. A: Red dotsmark the measurements in the 2–8 s time window. Individual measurements in the 2–4 s,4–6 s, and 6–8 s seconds lapse time windows and their standard deviation are indicatedby light gray dots and gray shading respectively. Results of the active CWI experiment(Wegler et al., 2006b) are indicated by green and blue lines. Inset shows a close up aroundthese active measurements. The blue line corresponds to the measurements at the GRWstation (used here). Green line shows measurements from station KEN on the southernflank of Mt. Merapi. B: Bottom part: daily precipitation rate (blue) and modeled groundwater level (black). Top part: measured (light red dots) and modeled (dark red line)velocity variations in the 2–4 s lapse time window. Same for the 6–8 s lapse time windowin blue.

4.7. CONCLUSIONS 83

−0.03

−0.02

−0.01

0

0.01

0.02

0.03

dv/v

Sep Nov 1998 Mar May Jul Sep Nov 1999 Mar May Jul−0.03

−0.02

−0.01

0

0.01

0.02

0.03

dv/v

A

B

Figure 4.7: A: Velocity variations measured from the TT-component of the GT betweenthe stations GRW0 and GRW1 in the 2-8 s lapse time window. B: Velocity variationmeasured on the auto-correlation of the Z-component from station GRW1.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

precipitation p(t)

drainage d(t)

vs

dept

h0

GW

L

δv

Figure 4.8: Illustration of the hydrological model that is used to explain the observedvelocity variations. On the left the velocity structure is sketched that shows a discontinuityat the depth of the ground water level. The hight of the ground water level on the righthand side is determined by the interplay between precipitation and drainage though aDarcy-like aquifer.


0 2 4

x 10−4

0

0.5

1

K(z)

dept

h [k

m]

τ = 3sτ = 7s

Figure 4.9: Amplitudes of the depth dependent sensitivity kernels for τ = 3 s and τ = 7 slapse times.

Chapter 5

Daily velocity variations andinstrumental errors at Merapivolcano

85


Abstract

Two case studies are presented that demonstrate the use of passive imaging for dif-ferent seismological purposes. In the first example the temporal resolution of PassiveImage Interferometry is improved to obtain measurements of temporal changes ona sub one day scale. Periodic changes of the Green’s functions are found that canbe related to temperature induced changes of the seismic velocities. In the sec-ond example Green’s functions on paths across Merapi volcano are retrieved. It isshown how these Green’s functions can be used to identify and correct instrumentalerrors in a network of four stations. A polarity reversal of the vertical componentat station GRW0 between the 25th September and the 25th October is identified.Measurements of the clock differences between station pairs are performed and usedto obtain detailed information about the timing errors of individual stations. Theseshow short offsets of about 1 s as well as periods of constant clock drifts. Relativetiming with this method can exceed the precision of an internal digitizer clock al-ready after the clock drift of a single day. The results are extendable to networks ofarbitrary size.

IMPROVING PASSIVE IMAGE INTERFEROMETRY 87

5.1 Improving the temporal resolution of passive

image interferometry

The accuracy of techniques that infer information from passive imaging depends onthe convergence rate of the field correlations toward a stable approximation of theGreen’s function (GF). This rate is determined by the characteristics of the wavefield that is used. The ideal field would be equipartitioned in phase space meaningthat all modes and directions are equally present (e.g. Sanchez-Sesma and Campillo,2006). This is seldom fulfilled. Especially the low frequency ambient seismic noise,that is often used for passive imaging in seismology, is excited predominantly in theoceans (Sabra et al., 2005; Stehly et al., 2006). Thus the sources have preferredlocations and frequencies. To compensate for this shortcoming the correlations areusually averaged over a long period of time. For structural investigations this doesnot pose a problem because it is assumed that the medium does not change withtime. In monitoring applications, however, the time needed for the noise correlationsto converge defines the temporal resolution of the method. It is thus desirable todecrease the necessary averaging time in order to increase the temporal resolutionof the monitoring.

5.1.1 Hourly Green’s functions

In chapter 4 I have shown that a temporal resolution of one day can be achieved withPassive Image Interferometry. The GFs obtained from 24 hours of continuous noiserecords at stations GRW0 and GRW1 (figure 4.1), as they were used in chapter 4, areshown in detail in figure 5.1. In contrast to the processing in chapter 4 I apply a one-bit normalization (Larose et al., 2004) before the correlation. The noise records arecorrelated in segments of one hour length. After correlation the 24 traces obtainedfor each day are stacked. The section in figure 5.1 consists of more than 600 traces(one for each day) with amplitudes mapped on a color scale. Several observationscan be made in this figure.

• The shape of the traces is very similar during the whole period of time. Promi-nent phases such as the ones at 3 s and 7 s lapse time can easily be tracedover the 2 years. This indicates that the medium remains stationary, i.e. thescatterers that generate the coda do not move.

• The velocity change reported in figure 4.6 can be observed in detail in theGFs. The Phase at 7 s lapse time for instance, exactly resembles the shape ofthe velocity variations shown in figure 4.6A.

• The delays for positive and negative delay times are symmetric. This indicatesthat the delays are indeed caused by medium changes and not by instrumentalerrors.

• The amplitudes of some phases in the GFs changes in July 1998. This changeprobably reflects variations in the characteristic of the seismic noise, that was

88 CHAPTER 5. DAILY VELOCITY VARIATIONS / INSTRUMENT ERRORS

laps

e tim

e [s

]

Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun

−8

−6

−4

−2

0

2

4

6

8

Figure 5.1: Time section of 24-hour Green’s functions between stations GRW0 andGRW1. There is one Green’s function for every day from beginning of August 1997 untilmid June 1999.

influenced by rock falls and accelerated seismic activity during a volcanic crisis.This change gradually disappears until beginning of January.

Figure 5.2 shows GFs obtained from correlations of one hour seismic noise.This is similar to figure 5.1 except that no averaging of the 1-hour GFs is made.There are more than 16000 GFs plotted in figure 5.2, one for each hour in which datais available. The similarity between figure 5.2 and figure 5.1 is remarkable. Thisindicates that even with averaging over only one hour of seismic noise, the correlationfunction appears to be a reasonable approximation of the GF. I therefore performthe same measurements as presented in chapter 4 on the one-hour GFs.


laps

e tim

e [s

]

Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun

−8

−6

−4

−2

0

2

4

6

8

Figure 5.2: Time section of 1-hour Green’s functions between stations GRW0 and GRW1.There is one Green’s function for every hour from beginning of August 1997 until mid June1999.

5.1.2 Hourly velocity variations

In figure 5.3 the velocity variations are shown for each one-hour GF together withthe velocity variations obtained from GF that are averaged over one day. The bluecurve is similar to figure 4.6A, except that the measurements are performed for theZ-component. The offset between end September and mid October 1997 is due toa polarity reversal at station GRW0 (cf. chapter 5.2). Compared to these dailymeasurements the red dots that represent the hourly measurements appear to showstronger fluctuations. A closer look at these fluctuations reveals some periodicity.Figure 5.4A shows an enlarged segment of the curve in figure 5.3. For comparisonthe GF for this period are shown in figure 5.4B. The similarity between the themeasured velocity variations and the shape of the Green’s functions in figure 5.4 isnot as striking as for the annual changes. This is due to the less averaging of the


Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun Jul−0.02

−0.015

−0.01

−0.005

0

0.005

0.01

0.015

0.02

δ v/

v

Figure 5.3: Hourly velocity variations measured between stations GRW0 and GRW1from August 1997 until mid June 1999. Red dots are measured from 1-hour GFs and blueline indicates the measurements obtained from GFs averaged over 24 hours.

A

03/15 03/16 03/17 03/18 03/19 03/20 03/21 03/22−15

−10

−5

0

5x 10−3

δ v/

v

B

laps

e tim

e [s

]

03/15 03/16 03/17 03/18 03/19 03/20 03/21

−8

−6

−4

−2

0

2

4

6

8

Figure 5.4: A: Enlarged period (15-22 March 1998) of the hourly velocity variationsmeasured between stations GRW0 and GRW1. Red line is measured from 1-hour GFsand blue line indicates the measurements obtained from GFs averaged over 24 hours. B:1-hour Green’s functions for this period.


0 2 4 6 8 10 12

velocity variations

0 2 4 6 8 10 12

air pressure variations

0 2 4 6 8 10 12

air temperature variations

0 2 4 6 8 10 12

precipitation rate

frequency [1/day]

Figure 5.5: Power spectra of hourly velocity variations and environmental processes thatmight influence the Green’s functions. The vertical scale is logarithmic.

one-hour GFs. In fact also the amplitude of the GRs seems to fluctuate periodically.This indicates that the noise sources are not constant.

To identify the cause of these fluctuations a spectral analysis of the velocityvariations is performed. For this analysis, velocity measurements are only used whenthe correlation of the GF with the reference trace exceeds 0.4. Values with lowercorrelation are substituted by linear interpolation of adjacent vales. In figure 5.5 thepower spectrum of the velocity variations is compared with different environmentalprocesses that might cause the changes in the GFs. These are rain, air temperature,and air pressure for which data form Merapi was provided by Malte Westerhaus.

The spectrum of the velocity variations shows a dominant peak at one cycle perday and 4 higher harmonics. There is a significant difference between the amplitudeof the 24 h period peak and the 12 h peak. The spectrum of the air pressurefluctuations has peaks at 24 h period and higher harmonics. In contrast to thevelocity variations the amplitudes of the 24 h and 12 h peaks are similar. It is a wellknown characteristic of air pressure variations that the semi-diurnal oscillation hasequal or higher amplitude then the diurnal. Air temperature variations again showdominant fluctuations with 24 h period and higher harmonics. The 12 h periodpeak is of considerably lower amplitude than the 24 h peak. This is similar to thespectrum of the velocity variations. The precipitation rate does only show one peakat 24 h period. Higher harmonics are missing. Since the 24 h periodicity seemsto be the fundamental oscillation in all these time series I stack each of the series


0 5 10 15 20t [h]

velocity variationair temperatureair pressurerain

Figure 5.6: Stack of the diurnal velocity variations, air-pressure, air temperature andprecipitation rate. Y-axis shows scaled amplitudes. Time axis is Coordinated UniversalTime (UTC) which is 7 h behind local time at Merapi.

with a 24 h periodicity to compare the phases. Therefor each time series is cut intoconsecutive segments of 24 h length. The stacks of these segments are shown infigure 5.6. Taking into account the 7 h difference between UTC used in the figureand local time at Merapi, on can recognize the expected shape of the air-temperatureand -pressure variations. The temperature reaches its maximum at about noon andits minimum before sunrise. The frequency doubling of the pressure variations isalso evident. The precipitation assumes its maximum in the afternoon. The velocityvariations have a shape similar to the temperature variations and the oscillationsare approximately in phase. The exact phase difference between the air-temperaturevariations and the velocity variations is 1 h with the temperature variations beingin advance.

5.1.3 Conclusions

With the increased temporal resolution of PII I identified periodic variations in theGFs. These variations concern the amplitude and the shape of the GFs. Applyingcoda wave interferometry to these GFs reveals periodic fluctuations of the seismicvelocity. The frequency spectra of the velocity variations and precipitation areclearly different. The 12 h oscillation and the higher harmonics are missing in theprecipitation data. I conclude that the diurnal velocity variations are not caused byor related to precipitation.

The spectra of air-pressure and -temperature variations are more similar tothe velocity variations. In particular the amplitude ratio between the 24 h and 12 hpeaks argues for a relation between the velocity variations and air-temperature. Thisis confirmed by the shape of the diurnal stack of the data which shows similaritybetween temperature and velocity data.

A possible mechanism that could cause variations of the seismic velocity asresponse to temperature variations is thermal stressing. The associated temperaturechanges in the subsurface are certainly limited to a very shallow layer but the stress


changes affect a larger volume. A combination of different causes like air pressureand temperature can not be excluded with the present analysis. A detailed modelinglike in chapter 4 for the annual changes is necessary to clearly identify the cause ofthe diurnal changes.


5.2 Detecting instrumental errors with passive

imaging

One of the most difficult tasks of operating a seismic network is the verification ofdata quality. Especially the data that is taken home from temporal network doesnot always have the expected characteristics. An interesting example is presented byMiller et al. (1998) who deployed a temporal network for the estimation of momenttensors and tomography. They had to detonate an explosion to check the polaritiesof the sensors of which 25% were wrong. Leaving these 25% uncorrected would bedisastrous for moment tensor estimation. Though tomography and event locationexperiments require high timing precision, the timing problems especially of tempo-ral stations are well known. Here a method is presented that uses passive imagingto verify the timing and polarity of a seismic station directly from the data.

5.2.1 Cross-volcano Green’s functions

To estimate GFs between between distant stations a lower frequency band is to beused. The strong attenuation of high frequency waves destroys the correlation ofthe wavefield on a larger distance. For the distance of 5 km between the Merapiarrays GRW and KEN, the frequency band between 0.1 Hz and 10 Hz proved to beuseful for the short period stations and between 0.01 Hz 10 Hz for the broadbandstations. Since this frequency band contains the microseismic peaks there is somevariability in the power spectrum that is smoothed by spectral whitening (Bensenet al., 2007) before the correlation. Figure 5.7 shows the color coded GF betweenthe stations GRW0 and KEN1 for the vertical components.

The following observations can be made.

1. The GFs show a good coherence over the whole period of time.

2. The GFs are not symmetric about zero.

3. Between end the 25th of September and the 15th of October 1997 the polarityof one of the instruments is reversed.

4. Between March 1998 and October 1999 the individual phases show a saw-toothpattern.

Item one shows that there is a good reproducibility of the GFs. Their shapedoes only marginally depend on the time of the year. This in agreement with Stehlyet al. (2006) who found stable noise sources in the frequency band 0.1 - 0.2 Hz.GFs obtained from noise with this characteristics are well suited for the study oftemporal changes. The asymmetry of the GFs mentioned in item two above indicatesthat either noise sources are unequally distributed or that the scatters which act assecondary sources have preferred locations (Paul et al., 2005; Stehly et al., 2006).The sign reversal mentioned in item three was already observed in chapter 5.1.2 forcorrelations between stations GRW0 and GRW1. This means that the reversal most

DETECTING INSTRUMENTAL ERRORS WITH PASSIVE IMAGING 95

laps

e tim

e [s

]

Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr May

−40

−30

−20

−10

0

10

20

30

40

Figure 5.7: Time section of the GF between stations GRW0 and KEN1 for the periodAugust 1997 until May 1999. The amplitudes are linearly amplified from zero lapse timetowards the ends.

probably occurred at GRW0. In contrast to the symmetric changes observed insection 5.1.1 the changes in figure 5.7 have the same sense for positive and negativelapse times. They are thus of different origin.

All these observations can be made even at late lapse times such as 100 s.Figure 5.8 shows a time section of a late phase in the GFs. This phase is found inall GFs between a GRW-station and a KEN-station.

5.2.2 Instrumental errors

Temporal changes that affect positive and negative times in the same way are due toerrors in one of the station’s clock (Stehly et al., 2007). To estimate the clock errorsof the individual stations I try to optimally use the configuration of the network. Thenetwork consists of two arrays: GRW and KEN (figure 4.1). Each array consists offour three-component stations (see figure 4.1 for the GRW array), that are pairwiseconnected to a 6-channel digitizer. Thus there are four digitizers with independentclocks and GPS timing system. Digitizers are named A with sensors KEN0 andKEN1, B with sensors KEN2 and KEN3, C with sensors GRW0 and GRW1 andD with sensors GRW2 and GRW3. To estimate to clock difference between twodigitizers, e.g. A and C, I use the GFs between the four sensor pairs that areconnected to the two digitizers, i.e. A1-C1, A1-C2, A2-C1, A2-C2, where A1 andA2 symbolize the two sensors connected with digitizer A.

I estimate the clock differences in a multi step procedure.

1. Calculation of a reference trace by stacking all GF for each sensor pair


laps

e tim

e [s

]

Aug Sep Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr May

−120

−115

−110

−105

−100

−95

−90

−85

−80

Figure 5.8: Time section of the GF between stations GRW0 and KEN1 for the periodAugust 1997 until May 1999.

2. Measurement of the clock difference as the position of the maximum of thecross-correlation function between the individual GFs and the reference tracesfor each sensor pair

3. Combination of the clock difference measurements from the four sensor pairs

4. Calculation of a new reference trace for each sensor pair like in 1 but shiftedto account for the clock difference

5. Measurement of the clock difference for each sensor pair like in 2 but with thenew reference trace.

6. Combination of the new clock difference measurements from the four sensorpairs

Steps four up to six are necessary because the reference trace generated in step 1 isaltered by the clock errors. These steps can be repeated several times to improve thereference trace, but after one iteration the improvements are marginal. For the cross-correlation in items 2 and 4 I use the complete GFs between -100 and 100 s lapse timeincluding direct and scattered waves. Because of the larger amplitudes the cross-correlation functions are dominated by the direct waves. For the combination ofthe clock difference measurements from the different sensor pairs a simple average iscalculated between measurements that fulfill the following criterion: The correlationbetween the reference trace and the GF, that was shifted to account for the measuredclock difference exceeds 0.4. The number of traces that are used in this averagingcan be increased by using more components of the Green’s tensors between the


station pairs. In fact in the present case one could use 36 component pairs for thisaveraging. Figure 5.9 shows the clock difference measurements for all digitizer pairs.

Similar measurements were obtained for stations in Southern California byStehly et al. (2007). This information is useful for analysis of the station pairs, likeambient noise tomography (Shapiro et al., 2005). But it does not help to improveevent localization that is strongly dependent on precise timing information of in-dividual stations. To provide this information the following approach is used here.The relation between the clock differences shown in figure 5.9 and the clock errorsof individual stations can be written in matrix form:

1 −1 0 01 0 −1 01 0 0 −10 1 −1 00 1 0 −10 0 1 −1

∆A

∆B

∆C

∆D

=

δABδACδADδBCδBDδCD

(5.1)

This equation is valid for each point in time where the clock differences are known.Here ∆X represents the clock error of digitizer X and δXY represents the clock differ-ence between digitizers X and Y. Equation 5.1 is a classical least squares inversionproblem of the form Gm = d with m being the unknown model vector and d be-ing the known data. The system is solved in the least squares sense to obtain theunknown clock offsets. These are shown in figure 5.10. Since it is still possible toadd an equal offset to all clocks, the curves are tied to the timing of digitizer C. Itsclock error is set to zero. Actually figure 5.9 shows no signature of a clock correctionthat is common to all clock difference measurements that contain digitizer C. Thisindicates that digitizer C either had a perfect clock or was, which is more likely,continuously locked to the GPS signal.

The curves in figure 5.10 show the characteristics of the individual clocks. Dig-itizer C probably has continuous GPS reception. Digitizer D also has a good GPSbut the worst clock. The GPS only lost the connection to the satellites in August1998 but the clock has a drift of 6.71e-7 or of 1.74 s/month. Digitizers A and Bhave clock drifts of 1.5e-7 and 9.33e-8 respectively. They lost their GPS connectionbetween March and October 1998. Between December 1997 and March 1998, thedigitizers A and B of the array KEN show frequent spikes of about 1 s amplitudethat can be attributed to instrumental errors. The two spikes of digitizer A arecommon with digitizer B. This indicates that these errors do not occur randomlybut are externally triggered. Probably the meteorological conditions caused poorGPS reception which led to misinterpretation of the–pulse–per second GPS signalat both digitizers at KEN. The accuracy of the method depends on the distancebetween the receivers, because of the damping of high frequency components overlonger distances. For stations pairs within one of the small arrays, i.e. with dis-tances of about 200 m, the standard deviation of the clock difference is about 10 ms.For longer distances of about 5 km across the volcano the standard deviation of theclock difference measurements is of the order of 50 ms.


Sep Nov 1998 Mar May Jul Sep Nov 1999 Mar May−1

−0.5

0

0.5

1cl

ock

diffe

renc

e [s

]

C1−A1C1−A2C2−A1C2−A2average


−0.5

0

0.5

1

cloc

k di

ffere

nce

[s]

C1−B1C1−B2C2−B1C2−B2average


−0.5

0

0.5

1

cloc

k di

ffere

nce

[s]

C1−D1C1−D2C2−D1C2−D2average


−0.5

0

0.5

1

cloc

k di

ffere

nce

[s]

D1−A1D1−A2D2−A1D2−A2average


−0.5

0

0.5

1

cloc

k di

ffere

nce

[s]

D1−B1D1−B2D2−B1D2−B2average


−0.5

0

0.5

1

cloc

k di

ffere

nce

[s]

A1−B1A1−B2A2−B1A2−B2average

Figure 5.9: Clock differences between the digitizers of the GRW and KEN stations. Dotsindicate the measurements between the different sensor pairs of a given digitizer pair. Redline shows the combination of the measurements from the corresponding sensor pairs.


Sep Nov 1998 Mar May Jul Sep Nov 1999 Mar May

−0.5

0

0.5

−0.5

0

0.5

−0.5

0

0.5

−0.5

0

0.5

Digitizer D

Digitizer C

Digitizer B

Digitizer A

cloc

k of

fset

[s]

Figure 5.10: Clock errors of the four digitizer clocks

5.2.3 Conclusions

The method presented here provides information about the polarity of the seismicrecords and about the clock offsets. This information is very valuable for the verifi-cation and analysis of seismic data but it can also be obtained in a more traditionalway. The polarity can for example be checked by comparison of first motions ofteleseismic arrivals or by explosions as done by Miller et al. (1998). In principle alsothe curves that show the clock offset in figure 5.10 can be reproduced by analysis ofGPS-logfiles that usually record GPS reception and clock corrections. Nevertheless Ithink it is very useful to be able to get all these information directly from the data,without the need for external information from logfiles or seismological analysis.Especially for the operation of temporal networks with difficult GPS conditions likein high latitudes, mountainous areas or on the seafloor the possibility to correct thetiming may prove very valuable. In OBS networks one can obtain precise relativetiming and average the clock drifts of the whole network which will increase theaccuracy of the absolute timing. Onshore one can relate the need for GPS receptionat all stations since a single station without GPS can easily be synchronized to theGPS-locked clocks of other stations. The accuracy of 10 ms for small receiver dis-tances (200 m) and 50 ms for distances of 5 km can already outstrip the accuracyof internal digitizer clocks after a single day.

100CHAPTER 5. DAILY VELOCITY VARIATIONS / INSTRUMENT ERRORS

Summary

This work contains a number of studies that extract information from scattered seis-mic waves. In part I the amplitude of the seismic coda is analyzed with radiativetransfer theory. The first study introduces an approach to simultaneously invertthe seismic coda for intrinsic attenuation, scattering strength, source spectrum, andsite amplification factors. This analysis is based on radiative transfer in the acous-tic approximation with isotropic scattering. It is shown that the estimated seismicmoments are in good agreement with results from waveform modeling, that site am-plifications correspond to geology, that the transport mean free path in the Germanis 690 km and that the quality factor of intrinsic attenuation is of the order of 500.

A Monte-Carlo algorithm is presented that was developed to solve the radiativetransfer equation for elastic waves. It accounts for non-isotropic scattering andconversion between P- and S-modes. The large scale velocity structure of the crustis taken into account in terms of different velocities in crust and mantle, and thediscontinuity in-between, that leads to conversion, reflection, and transmission ofseismic energy. Reflection and conversion at the free surface is also treated. Avertical velocity gradient in the mantle can be used to generate mantle phases likePn and Sn.

The study on Lg-blockage applies this algorithm to explain the reason of thestrong attenuation of crustal phases in the western Pyrenees. In simulations of theradiative transfer process in a 3-dimensional structure it is shown that increasedheterogeneity in a localized region below the western Pyrenees is able to explain theLg-blockage. The fit of the model can be improved by additionally increasing theintrinsic attenuation in the western Pyrenean body. Intrinsic attenuation alone cannot account for the observed attenuation.

The second part of this work uses the phase information of the seismic coda.Additionally the seismic coda is not simply measured in a classical source–receiverconfiguration, but retrieved from the ambient seismic noise by cross-correlation ofthe signals from two receivers. This process of Green’s function retrieval is basedon the phase information of the seismic noise.

In the first study of data from Merapi Volcano it is shown that the passivelyretrieved Green’s functions can be retrieved from noise records of 24 h length. Thisallows to obtain a seismic record independent of active sources once a day. Basedon this finding a monitoring technique named passive image interferometry is intro-duced. It allows to precisely monitor velocity changes in the subsurface. Applicationto noise records from Merapi volcano shows velocity variations of a few percent ona time scale of weeks. A strong seasonal trend is observed. In a sophisticated in-

101

102 SUMMARY

version process, the parameters of a hydrological model are estimated that relatesthe measured velocity variations with precipitation induced changes of the groundwater level. With the final parameters the hydrological model can explain the shapeof the velocity variations and their lapse time dependent amplitudes solely based onprecipitation data.

The temporal resolution of passive image interferometry is increased in anotherstudy from one day to one hour. Based on spectral comparison with various mete-orological parameters, this study reports evidence for a relation of the fast velocityvariations with atmospheric temperature.

The last study introduces a method to assess the accuracy of instrument timingand setup directly from records of seismic noise. The signature of timing errors onthe passively retrieved Green’s functions is different from the signature of velocityvariations, that were analyzed in the studies above. The method above is appliedto data from four digitizers from the Merapi network. It is shown how the clockdifference of each instrument pair can be estimated from the noise records. With aleast squares inversion the clock offset of each digitizer is obtained.

Bibliography

Abubakirov, I. and Gusev, A. (1990). Estimation of scattering properties of thelithosphere of Kamchatka based on Monte-Carlo simulation of record envelope ofa near earthquake. Phys. Earth Planet. Int., 64(1):52–67.

Abubakirov, I. R. (2005). Attenuation characteristics of transverse waves in thelithoshere of Kamchatka estiamted from observations at the Petropavlovsk digitalbroadband station. Izvestiya, Physics of the Solid Earth, 41(10):813–824.

Akasaka, C. and Nakanishi, S. (2000). Correction of background gravity changes dueto precipitation: Oguni geothermal field, Japan. Proceedings World GeothermalCongress, pages 2471–2475.

Aki, K. (1980). Attenuation of shear-waves in the lithosphere for frequencies from0.05 to 25 Hz. Phys. Earth Planet. Inter., 21(1):50–60.

Aki, K. (1992). Scattering conversions P to S versus S to P. Bull. Seismol. Soc.Am., 82(4):1969–1972.

Aki, K. and Chouet, B. (1975). Origin of coda waves: source, attenuation, andscattering effects. J. Geophys. Res., 80(23):3322–3342.

Aki, K. and Richards, P. G. (1980a). Quantitative Seismology, Theory and Methods,volume II. W. H. Freeman and Company, San Francisco.

Aki, K. and Richards, P. G. (1980b). Quantitative Seismology, Theory and Methods,volume I. W. H. Freeman and Company, San Francisco.

Apresyan, L. and Kravtsov, Y. A. (1996). Radiation Transfer: Statistical and WaveAspects. Gordon and Breach, Amsterdam.

Baumgardt, D. R. (2001). Sedimentary basins and the blockage of lg wave propa-gation in the continents. Pure Appl. Geophys., 158:1207–1250.

Bensen, G. D., Ritzwoller, M. H., Barmin, M. P., Levshin, A. L., Lin, F., Moschetti,M. P., Shapiro, N. M., and Yang, Y. (2007). Processing seismic ambient noise datato obtain reliable broad-band surface wave dispersion measurements. Geophys. J.Int., 169(3):1239–1260. doi: 10.1111/j.1365-246X.2007.03374.x.

103

104 BIBLIOGRAPHY

Bianco, F., Del Pezzo, E., Castellano, M., Ibanez, J., and Di Luccio, F. (2002). Sep-aration of intrinsic and scattering seismic attenuation in the Southern Apenninezone, Italy. Geophys. J. Int., 150:10–22.

Birch, F. (1961). The velocity of compressional waves in rocks to 10 kilobars, part2. J. Geophys. Res., 66:2199–2224.

Bormann, P., editor (2002). New Manual of Seismological Observatory Practice.GFZ, IASPAI, Potsdam.

Braunmiller, J. (2002). Moment tensor solutions of stronger earthquakes in Ger-many with GRSN data. In Korn, M., editor, Ten Years of German RegionalSeismic Network (GRSN), number Report 25, pages 227–235. Senate Commissonfor Geoscience of the Deutsche Forschungsgemeinschaft.

Braunmiller, J., Dahm, T., and Bonjer, K. (1994). Source mechanism of the 1992Roermond earthquake from surface wave inversion of regional data. Geophys. J.Int., 116:663–672.

Braunmiller, J., Kradolfer, U., Baer, M., and Giardini, D. (2002). Regional mo-ment tensor determination in the European-Mediterranean area - initial results.Tectonophysics, 356:5–22.

Buttkus, B. (2000). Spectral Analysis and Filter Theory in Applied Geophysics.Springer-Verlag Berlin, Heidelberg.

Campillo, M., Feigner, B., Bouchon, M., and Bethoux, N. (1993). Attenuation ofcrustal waves across the Alpine range. J. Geophys. Res., 98(B2):1987–1996.

Campillo, M. and Paul, A. (2003). Long-range correlations in the diffuse seismiccoda. Science, 299:547 – 549. doi:10.1126/science.1078551.

Chandrasekhar, S. (1960). Radiative Transfer. Dover, New York.

Chazalon, A., Campillo, M., Gibson, R., and Carreno, E. (1993). Crustal wavepropagation anomaly across the pyrenean range. comparison between observationsand numerical simulations. Geophys. J. Int., 115:829–838.

Choukroune, P. (1992). Tectonic Evolution of the Pyrenees. Ann. Rev. Earth Planet.Sci., 20:143–158.

Daignieres, M., Gallart, J., and Banda, E. (1981). Lateral variations of the crust inthe North Pyrenean Zone. Ann. Geophys.

Dainty, A. M. and Toksoz, M. N. (1981). Seismic codas on the earth and the moon:a comparison. Phys. Earth Planet. Int., 26:250–260.

Derode, A., Larose, E., Campillo, M., and Fink, M. (2003). How to esti-mate the Green’s function of a heterogeneous medium between two passivesensors? Application to acoustic waves. Appl. Phys. Lett., 83(15):3054–3056.doi:10.1063/1.1617373.

BIBLIOGRAPHY 105

Dewberry, R. S. and Crosson, R. S. (1995). Source scaling and moment estimationfor the Pacific northwest network using S-coda amplitudes. Bull. Seismol. Soc.Am., 85(5):1309–1326.

Duvall, T. L., Jeffferies, S. M., Harvey, J. W., and Pomerantz, M. A. (1993). Time-distance helioseismology. Nature, 362:430 – 432. doi:10.1038/362430a0.

Dziewonski, A. M. and Anderson, D. L. (1981). Preliminary reference earth model.Phys. Earth Planet. Int., 25:297–356.

Fehler, M., Hoshiba, M., Sato, H., and Obara, K. (1992). Separation of scatteringand intrinsic attenuation for the Kanto-Tokai region, Japan, using measurementsof S-wave energy versus hypocentral distance. Geophys. J. Int., 108(3):787–800.

Gret, A., Snieder, R., and Scales, J. (2006). Time-lapse monitoring of rockproperties with coda wave interferometry. J. Geophys. Res., 111(B03305).doi:10.1029/2004JB003354.

Gusev, A. A. and Abubakirov, I. R. (1987). Monte-carlo simulation of record enve-lope of a near earthquake. Phys. Earth Planet. Int., 49:30–36.

Gusev, A. A. and Abubakirov, I. R. (1996). Simulated envelopes of non-isotropicallyscattered body waves as compared to observed ones: another manifestation offractal heterogeneity. Geophys. J. Int., 127(1):49–60.

Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains andtheir applications. Biometrika, 57:97–109.

Hennino, R.and Tr´A c©goures, N., Shapiro, N. M., Margerin, L., Campillo, M., vanTiggelen, B., and Weaver, R. L. (2001). Observation of equipartition of seismicwaves. Phys. Rev. Lett., 86:3447–3450. doi:10.1103/PhysRevLett.86.3447.

Hoshiba, M. (1991). Simulation of multiple-scattered coda wave excitation based onthe energy conservation law. Phys. Earth Planet. Inter., 67(1-2):123–136.

Hoshiba, M. (1995). Estimation of nonisotropic scattering in western Japan usingcoda wave envelopes: Application of a multiple nonisotropic scattering model. J.Geophys. Res., 100(B1):645–657.

Hoshiba, M., Rietbrock, A., Scherbaum, F., Nakahara, H., and Haberland, C. (2001).Scattering attenuation and intrinsic absorption using uniform and depth depen-dent model - application to full seismogram envelope recorded in northern Chile.J. Seismology, 5(2):157–179.

Hoshiba, M., Sato, H., and Fehler, M. (1991). Numerical basis of the separation ofscattring and intrinsic absorption from full seismogram envelope - A Monte-Carlosimulation of multiple isotropic scattering. Pa. Meteorol. Geophys., Meteorol. Res.Inst., 42:65–91.

106 BIBLIOGRAPHY

Hough, S., Dollar, R., and Johnson, P. (2000). The 1998 earthquake sequence southof Long Valley caldera, California: hints of magmatic involvement. Bull. Seismol.Soc. Am., 90(3):752–763.

Huffman, G., Adler, R., Morrissey, M., Curtis, S., Joyce, R., McGavock, B., andSusskind, J. (2001). Global precipitation at one-degree daily resolution frommulti-satellite observations. J. Hydrometeor., 2:36–50.

Ishimaru, A. (1978). Wave Propagation and Scattering in Random Media, volume1 and 2. Academic Press, New York.

Korn, M., editor (2002). Ten Years of German Regional Seismic Network (GRSN).Number Report 25. Senate Commisson for Geoscience of the Deutsche Forschungs-gemeinschaft.

Lacombe, C., Campillo, M., Paul, A., and Margerin, L. (2003). Separation of intrin-sic absorption and scattering attenuation from Lg coda decay in central Franceusing acoustic radiative transfer theory. Geophys. J. Int., 154(2):417–425.

Larose, E., De Rosny, J., Margerin, L., Anache, D., Gouedard, P., Campillo, M.,and Van Tiggelen, B. (2006). Observation of multiple scattering of kHz vibrationsin a concrete structure and application to monitoring weak changes. Phys. Rev.E, 73:016609. doi:10.1103/PhysRevE.73.016609.

Larose, E., Derode, A., Campillo, M., and Fink, M. (2004). Imaging from one-bit correlations of wideband diffusive wave fields. Journal of Applied Physics,95(12):8393–8399. doi:10.1063/1.1739529.

Larose, E., Khan, A., Nakamura, Y., and Campillo, M. (2005). Lunar subsurfaceinvestigated from correlation of seismic noise. Geophys. Res. Lett., 32:L16201.doi:10.1029/2005GL023518.

Lobkis, O. I. and Weaver, R. L. (2001). On the emergence of the Green function inthe correlation of a diffuse field. J. Acoust. Soc. Am., 110:3011–3017.

Lux, I. and Kolbinger, L. (1991). Monte Carlo particle transport and methods:neutron and photon calculations. CRC Press.

Margerin, L. (2005). Introduction to radiative transfer of seismic waves. In Levan-der, A. and Nolet, G., editors, Seismic data analysis with global and local arrays,volume 157 of AGU Monograph Series, chapter 14, pages 229–252. AGU, Wash-ington.

Margerin, L., Campillo, M., and Van Tiggelen, B. (2000). Monte Carlo simulationof multiple scattering of elastic waves. J. Geophys. Res., 105(B4):7873–7892.

Margerin, L., Campillo, M., and Van Tiggelen, B. A. (1998). Radiative transfer anddiffusion of waves in a layered medium: new insight into coda Q. Geophys. J.Int., 134:596–612.

BIBLIOGRAPHY 107

Mayeda, K., Hofstetter, A., O’Boyle, J., and Walter, W. (2003). Stable and trans-portable regional magnitudes based on coda-derived moment-rate spectra. Bull.Seismol. Soc. Am., 93(1):224–239.

Mayeda, K., Koyanagi, S., Hoshiba, M., Aki, K., and Zeng, Y. (1992). A comparativestudy of scattering, intrinsic, and coda Q−1 for Hawaii, Long Valley Caldera, andcentral California between 1.5 and 15.0 hz. J. Geophys. Res., 97(B5):6643–6659.

Mayeda, K. and Walter, W. (1996). Moment, energy, stress drop, and source spectraof western United States earthquakes from regional coda envelopes. J. Geophys.Res., 101(B5):11195–11208.

McNamara, D. E. and Walter, W. R. (2001). Mapping Crustal Heterogeneity UsingLg Propagation Efficiency Throughout the Middle East, Mediterranean, SouthernEurope and Northern Africa. Pure and Applied Geophysics, 158:1165–1188.

Miller, A. D., Julian, B. R., and Foulger, G. R. (1998). Three-dimensional seismicstructure and moment tensors of non-double-couple earthquakes at the hengill-grensdalur volcanic cmoplex, iceland. Geophys. J. Int., 133:309–325.

Morasca, P., Mayeda, K., Malagnini, L., and Walter, W. (2005). Coda-derivedsource spectra, moment magnitudes and energy-moment scaling in the westernAlps. Geophys. J. Int., 160:263–275.

Nakahara, H., Nishimura, T., Sato, H., and Ohtake, M. (1998). Seismogram envelopeinversion for the spatial distribution of high-frequency energy radiation from theearthquake fault: Application to the 1994 far east off Sanriku earthquake, Japan.J. Geophys. Res., 103(B1):855–867.

Nishimura, T., Uchida, N., Sato, H., Ohtake, M., Tanaka, S., and Hamaguchi,H. (2000). Temporal changes of the crustal structure associated with the M6.1earthquake on September 3, 1998 and the volcanic activity of mount Iwate, Japan.Geophys. Res. Lett., 27(2):269–272.

Paasschens, J. (1997). Solution of the time-dependent Boltzmann equation. Phys.Rev. E, 56(1):1135–1141.

Pacheco, C. and Snieder, R. (2005). Time-lapse travel time change of mul-tiply scattered acoustic waves. J. Acoust. Soc. Am., 118(3):1300–1310.doi:10.1121/1.2000827.

Paul, A., Campillo, M., Margerin, L., and Larose, E. (2005). Empirical synthesisof time-asymmetrical Green functions from the correlation of coda waves. J.Geophys. Res., 110:B08302. doi:10.1029/2004JB003521.

Peng, Z. and Ben-Zion, Y. (2006). Temporal changes of shallow seismic veloc-ity around the Karadere-Duzce branch of the North Anatolian Fault and strongground motion. Pure Appl. Geophys., 163:567–600. doi:10.1007/s00024-005-0034-6.

108 BIBLIOGRAPHY

Poupinet, G., Ellsworth, W. L., and Frechet, J. (1984). Monitoring velocity vari-ations in the crust using earthquake doublets: An application to the Calaverasfault, California. J. Geophys. Res., 89(B7):5719–5731.

Pous, J., Ledo, J., Marcuello, A., and Daignieres, M. (1995). Electrical resistivetymodel of the crust and upper mantle from manetotelluric survey through thecentral Pyrenees. Geophys. J. Int., 121:750–762.

Przybilla, J., Korn, M., and Wegler, U. (2006). Radiative transfer of elastic wavesversus finite difference simulations in two-dimensional random media. J. Geophys.Res., 111:B04305. doi:10.1029/2005JB003952.

Ratdomopurbo, A. and Poupinet, G. (1995). Monitoring a temporal change ofseismic velocity in a volcano: application to the 1992 eruption of Mt. Merapi(Indonesia). Geophys. Res. Lett., 22(7):775–778.

Rautian, T. and Khalturin, V. (1978). The use of the coda for determination of theearthquake source spectrum. Bull. Seismol. Soc. Am., 68(4):923–948.

Reamer, S. and Hinzen, K.-G. (2004). An earthquake catalog for the Northern Rhinearea, Central Europe (1975-2002). Seismol. Res. Lett., 75(6):713–725.

Rickett, J. and Clearbout, J. (1999). Acoustic daylight imaging via spectral factor-ization: Helioseismology and reservoir monitoring. The Leading Edge, 18(2):957–960.

Rodgers, A. J., Ni, J. F., and Hearn, T. M. (1997). Propagation characteristics ofshort-period sn and lg in the middle east. Bull. Seismol. Soc. Am., 87(2):396–413.

Roux, P., Sabra, K. G., Gerstoft, P., Kuperman, W. A., and Fehler, M. C. (2005).P-waves from cross-correlation of seismic noise. Geophys. Res. Lett., 32:L19303.doi:10.1029/2005GL023803.

Ryzhik, L., Papanicolaou, G., and Keller, J. B. (1996). Transport equations forelastic and other waves in random media. Wave Motion, 24:327–370.

Sabra, K. G., Gerstoft, P., Roux, P., Kuperman, W. A., and Fehler, M. (2005).Surface wave tomography from microseisms in Southern California. Geophys.Res. Lett., 32:L14311. doi:10.1029/2005GL023155.

Saito, T., Sato, H., and Ohtake, M. (2002). Envelope broadening of sphericallyoutgoing waves in three-dimensional random media having power-law spectra. J.Geophys. Res., 107(B5):3–1 – 3–16. doi:10.1029/2001JB000264.

Sanchez-Sesma, F. J. and Campillo, M. (2006). Retrieval of the Green functionfrom cross-correlation: The canonical elastic problem. Bull. Seismol. Soc. Am.,96:1182–1191.

BIBLIOGRAPHY 109

Sato, H. (1989). Broadening of seismogram envelopes in the randomly inhomoge-neous lithosphere based on the parabolic approximation: southeastern Honshu,Japan. J. Geophys. Res., 94(B12):17,735–17,747.

Sato, H. and Fehler, M. C. (1998). Seismic Wave Propagation and Scattering in theHeterogeneous Earth. Springer-Verlag, New York.

Sato, H., Nakahara, H., and Ohtake, M. (1997). Synthesis of scattered energydensity for non-spherical radiation from a point shear dislocation source based onthe radiative transfer theory. Phys. Earth Planet. Inter., 104:1–13.

Sens-Schonfelder, C. and Wegler, U. (2006a). Passive image interferomety and sea-sonal variations of seismic velocities at Merapi volcano, Indonesia. Geophys. Res.Lett., 33:L21302. doi: 10.1029/2006GL027797.

Sens-Schonfelder, C. and Wegler, U. (2006b). Radiative transfer theory for es-timation of the seismic moment. Geophys. J. Int., 167(3):1363–1372. doi:10.1111/j.1365-246X.2006.03139.x.

Shapiro, N., Bethoux, N., Campillo, M., and Paul, A. (1996). Regional seismicphases across the Ligurian Sea: Lg blockage and oceanic propagation. Phys.Earth Planet. Int., 93:257–268.

Shapiro, N. M. and Campillo, M. (2004). Emergence of broadband Rayleigh wavesfrom correlations of the ambient seismic noise. Geophys. Res. Lett., 31:L07614,.doi:10.1029/2004GL019491.

Shapiro, N. M., Campillo, M., Stehly, L., and Ritzwoller, M. H. (2005). High-resolution surface-wave tomography from ambient seismic noise. Science,307:1615–1618. doi:10.1126/science.1108339.

Shearer, P. and Earle, P. (2004). The global short-period wavefield modelledwith a monte carlo seismic phonon method. Geophys. J. Int., 158:1103–1117.doi:10.1111/j.1365-246X.2004.02378.x.

Shearer, P. M. (1999). Introduction to Seismology. Cambridge University Press,Cambridge.

Snieder, R. (2004). Extracting the Green’s function from the correlation of codawaves: A derivation based on stationary phase. Phys. Rev. E, 69:046610.doi:10.1103/PhysRevE.69.046610.

Snieder, R. (2006). The theory of coda wave interferometry. Pure Appl. Geophys.,163:455–473. doi:10.1007/s00024-005-0026-6.

Snieder, R., Gret, A., and Scales, J. (2002). Coda wave interferometry for estimatingnonlinear behavior in seismic velocity. Science, 295:2253–2255.

Snieder, R. and Page, J. (2007). Multiple scattering in evolving media. Phys. Today,60(4):49–55.

110 BIBLIOGRAPHY

Souriau, A. and Granet, M. (1995). A tomographic study of the lithosphere beneaththe Pyrenees from local and teleseismic data. J. Geophys. Res., 100:18117–18134.

Souriau, A., Sylvander, M., Rigo, A., Fels, J., Douchain, J., and Ponsolles, C. (2001).Sismotectonique des Pyrenees: principales contraintes sismologiques. Bull. Soc.Geol. Fr., 172(1):25–39.

Stehly, L., Campillo, M., and Shapiro, N. M. (2006). A study of the seismicnoise from its long-range correlation properties. J. Geophys. Res., 111:B10306.10.1029/2005JB004237.

Stehly, L., Campillo, M., and Shapiro, N. M. (2007). Travel time measurementsfrom noise correlations: stability and detection of instrumental errors. Geophys.J. Int. doi: 10.1111/j.1365-246X.2007.03492.x.

Turner, J. A. (1998). Scattering and diffusion of seismic waves. Bull. Seismol. Soc.Am., 88(1):276–283.

Vacher, P. and Souriau, A. (2001). A three-dimensional model of the Pyrenean deepstructure based on gravity modelling, seismic images and petrological constraints.Geophys. J. Int., 145:460–470.

Wapenaar, K. (2004). Retrieving the elastodynamic Green’s function of an arbi-trary inhomogeneous medium by cross correlation. Phys. Rev. Lett., 93:254301.doi:10.1103/PhysRevLett.93.254301.

Wassermann, J. and Ohrnberger, M. (2001). Automatic hypocenter determinationof volcano induced transients based on wavefield coherence - an application to the1998 eruption of Mt. Merapi, indonesia. J. Volcanol. Geotherm. Res., 110:57–77.

Weaver, R. L. (1990). Diffusivity of ultrasound in polycrystals. J. Mech. Phys.Solids, 38:55–86.

Weaver, R. L. and Lobkis, O. I. (2001). Ultrasonics without a source: Ther-mal fluctuation correlation at MHz frequencies. Phys. Rev. Lett., 87(13):134301.doi:10.1103/PhysRevLett.87.134301.

Wegler, U. (2004). Diffusion of seismic waves in a thick layer: Theory and applicationto Vesuvius volcano. J. Geophys. Res., 109(B7):10.1029/2004JB003048.

Wegler, U., Korn, M., and Przybilla, J. (2006a). Modelling full seismogram en-velopes using radiative transfer theory with Born scattering coefficients. PureAppl. Geophys., 163:503–531.

Wegler, U. and Luhr, B. G. (2001). Scattering behaviour at Merapi volcano (Java)revealed from an active seismic experiment. Geophys. J. Int., 145(3):579–592.doi:10.1046/j.1365-246x.2001.01390.x.

BIBLIOGRAPHY 111

Wegler, U., Luhr, B. G., Snieder, R., and Ratdomopurbo, A. (2006b). Increaseof shear wave velocity before the 1998 eruption of Merapi volcano (Indonesia).Geophys. Res. Lett., 33:L09303. doi:10.1029/2006GL025928.

Wegler, U. and Sens-Schonfelder, C. (2007). Fault zone monitoring with passiveimage interferometry. Geophys. J. Int., 168(3):1029–1033. doi: 10.1111/j.1365-246X.2006.03284.x.

Wu, R. and Aki, K. (1985). The fractal nature of the inhomogeneities in the litho-sphere evidenced from seismic wave scattering. Pure Appl. Geophys., 123(6):805–818.

Zeng, Y., Su, F., and Aki, K. (1991). Scattering wave energy propagation in arandom isotropic scattering medium: 1. Theory. J. Geophys. Res., 96(B1):607–619.

112 BIBLIOGRAPHY

Deutsche Zusammenfassung

113

DEUTSCHE ZUSAMMENFASSUNG 115

Gestreute Wellenfelder werden in der Seismologie im Allgemeinen als Storsignale an-gesehen. Aufgrund der Komplexitat dieser Wellenfelder sind Strahl- und Wellentheo-rie, die Werkzeuge der klassischen Seismologie, fur deren Beschreibung ungeeignet(Aki and Richards, 1980). Damit ist der uberwiegende Teil der seismologischen Da-ten, namlich sowohl die Coda von Erdbeben als auch das seismische Umgebungsrau-schen, auf traditionallem Wege nicht zur Informationsgewinnung nutzbar. Die vor-liegende Arbeit beschaftigt sich mit zwei alternativen Ansatzen, die es ermoglichen,diesen enormen Datenbestand fur geophysikalische Zwecke zu nutzen.

Energietransfertheorie

Im ersten Teil der Arbeit werden Erdbebenaufzeichnungen auf Grundlage der Ener-gietransfertheorie (Ryzhik et al., 1996) ausgewertet. Diese Theorie betrachtet nichtdie Ausbreitung von elastischen Wellen, sondern vernachlassigt deren Phasenbe-ziehungen in der Annahme, dass diese unkorreliert sind. Damit kann die Energie-transfertheorie nur Aussagen zur Ausbreitung der Energie des Wellenfeldes treffen.Diesem Nachteil steht der Vorteil gegenuber, dass durch die Vernachlassigung derPhasenbeziehungen auf eine deterministische Beschreibung des Ausbreitungsmedi-ums verzichtet werden kann. Damit konnen die kleinskaligen Heterogenitaten, diefur die Ausbildung der seismischen Coda verantwortlich sind, statistisch beschriebenund die Form der Seismogrammeinhullenden (Envelope) einschließlich Coda model-liert werden.

Akustische Naherung der Energietransfertheorie

Hier wird die Energietransfertheorie in der akustischen Naherung benutzt. Die Aus-breitung der Scherwellenenergie lasst sich in diesem Fall durch eine skalare Energie-transfergleichung beschreiben. Die Energie der Longitudinalwelle wird vernachlassigt.Es wird eine Methode entwickelt, die es ermoglicht, aus den Aufzeichnungen derScherwellen und deren Coda von mehreren Erdbeben an mehreren Stationen dasseismische Moment der Erdbeben, die Verstarkungsfaktoren des Untergrundes dereinzelnen Stationen sowie die Dampfungsparameter der Erdkruste zu bestimmen.Insbesondere lassen sich die Effekte von intrinsischer Dampfung und Streudampfungtrennen.In einer Anwendung auf Breitbanddaten des Deutschen Regionalnetzes wird die Me-thodik demonstriert und es werden Informationen uber die Dampfungseigenschaf-ten der Kruste sowie die Verstarkungsfaktoren des Untergrundes der Regionalnetz-Stationen gewonnen. Fur die mittlere freie Weglange, die die Starke der Streuungbeschreibt, wurde ein Wert von 690 km bestimmt. Der Qualitatsfaktor der intrin-sischen Dampfung liegt im Frequenzbereich um 3 Hz bei 500. Die Verstarkungs-faktoren zeigen einen eindeutigen Bezug zur Geologie. Durch einen Vergleich dergemessenen seismischen Momente mit unabhangigen Messungen durch Wellenform-inversion werden die Inversionsergebnisse verifiziert.Wesentliche Bestandteile dieser Untersuchung wurden in Sens-Schonfelder and Weg-ler (2006b) veroffentlicht.

116 SUMMARY

Monte-Carlo Simulationen fur die Ausbreitung elastischer Energie

Um die Einschrankungen zu uberwinden, die sich aus der Anwendung der akusti-schen Naherung ergeben, wird in diesem Kapitel ein Algorithmus entwickelt, der dieEnergietransfergleichung im elastischen Fall lost. Das Energietransferproblem wirddafur zunachst in Form von drei gekoppelten skalaren Gleichungen formuliert, diesich fur die Implementierung eines Monte-Carlo-Algorithmus’ eignen. Die einzelnenKomponenten des Algorithmus’ werden im Anschluss diskutiert. Der Algorithmusermoglicht die Berechnung von vollstandigen Seismogrammeinhullenden in einemModell, das durch statistische Eigenschaften der kleinskaligen Heterogenitaten so-wie durch deterministische Struktur beschrieben wird. Die Berechnung erfolgt unterBeachtung von Konversion zwischen Longitudinal- und Transversalwellen und vonReflexion/Transmission und Konversion an der freien Oberflache und einer Diskon-tinuitat in der Tiefe.

Elastische Energietransfertheorie und die Ausloschung von Lg-Wellen inden Pyrenaen

Der im vorangegangenen Kapitel beschriebene Algorithmus kann Seismogramm-einhullende in einem Modell berechnen, das neben den kleinskaligen Heterogenitatenauch deterministische großskalige Struktur enthalt. Diese Fahigkeit erlaubt eine sehrrealistische Beschreibung von Aufzeichnungen regionaler Erdbeben und wird hierausgenutzt, um außergewohnliche Dampfungseigenschaften der Erdkruste unterhalbder westlichen Pyrenaen zu untersuchen. In dem Bereich werden durch die Krustelaufende Wellen ungewohnlich stark gedampft. Dieses, nach einer krustalen PhaseLg-Ausloschung genannte, Phanomen ist von marinen Sedimentbecken bekannt, indenen es durch das Einfangen von Energie in den Sedimenten mit sehr geringenGeschwindigkeiten und hoher Dampfung verursacht wird. Auf die Situation in denPyrenaen ist diese Erklarung allerdings nicht anwendbar und es wurde gezeigt, dassdie lokale makroskopische Geschwindigkeitsstruktur nur einen sehr geringen Ein-fluss auf die Amplituden der krustalen Phasen hat. Damit war der Mechanismus derLg-Ausloschung in den Pyrenaen bislang ungeklart.Hier wird ein Modell vorgestellt, das zusatzlich zur normalen Kruste–Mantel–Struktureinen Storkorper in der Kruste unterhalb der westlichen Pyrenaen enthalt. DieserKorper gleicht in seinen makroskopischen Eigenschaften wie Dichte und mittlererGeschwindigkeit der seismischen Wellen der umgebenen Kruste. Unterschiede be-stehen in der intrinsischen Dampfung und der Starke der Streuung. In einer erstenInversionsphase werden die Parameter der ungestorten Kruste aus Aufzeichnungenvon Wellenwegen durch die ostlichen Pyrenaen bestimmt. Auf diesen Parameternaufbauend, werden die Aufzeichnungen von Wellenwegen durch die westlichen Py-renaen nach den Parametern des Storkorpers invertiert.Dabei werden drei Modelle getestet: 1. Die Lg-Ausloschung wird durch intrinsischeDampfung verursacht. 2. Die Ausloschung wird durch Streuung verursacht. Und 3.sowohl Streuung als auch intrinsische Dampfung tragen zur Ausloschung der krus-talen Phasen bei. Alle drei Modelle sind in der Lage, die durch die Kruste laufendeEnergie effizient zu dampfen, aber es gibt deutliche Unterschiede in der Form der


A

-6˚

-6˚

-5˚

-5˚

-4˚

-4˚

-3˚

-3˚

-2˚

-2˚

-1˚

-1˚

0˚

0˚

1˚

1˚

2˚

2˚

3˚

3˚

4˚

4˚

5˚

5˚

6˚

6˚

37˚ 37˚

38˚ 38˚

39˚ 39˚

40˚ 40˚

41˚ 41˚

42˚ 42˚

43˚ 43˚

44˚ 44˚

45˚ 45˚

46˚ 46˚

47˚ 47˚

48˚ 48˚

49˚ 49˚

200 km

low high

amplitudeB

-6˚

-6˚

-5˚

-5˚

-4˚

-4˚

-3˚

-3˚

-2˚

-2˚

-1˚

-1˚

0˚

0˚

1˚

1˚

2˚

2˚

3˚

3˚

4˚

4˚

5˚

5˚

6˚

6˚

37˚ 37˚

38˚ 38˚

39˚ 39˚

40˚ 40˚

41˚ 41˚

42˚ 42˚

43˚ 43˚

44˚ 44˚

45˚ 45˚

46˚ 46˚

47˚ 47˚

48˚ 48˚

49˚ 49˚

200 km

low high

amplitude

Abbildung 1: Momentaufnahmen der Energieverteilung 105 s (A) und 180 s (B) nacheinem Erdbeben in Spanien. Stern: Epizentrum. Rechteck: Storkorper.

resultierenden Seismogrammeinhullenden. Es wird gezeigt, dass Modell 1 eine deutli-che Diskrepanz zu den Beobachtungen hinterlasst, woraus geschlussfolgert wird, dassintrinsische Dampfung nicht die alleinige Ursache der Lg-Ausloschung ist. Modell 2dagegen ist in der Lage, die Beobachtungen vernunftig zu reproduzieren. Die Kom-bination aus verstarkter Streuung und erhohter intrinsischer Dampfung in Modell 3kann die Beobachtungen am besten erklaren.

Fur die Parameter des besten Modells zeigen Abbildungen 1A und B Momentauf-nahmen der Energieverteilung 105 s und 180 s nach einem Erdbeben in Spanien. DieBilder verdeutlichen die Dampfung der durch die Kruste laufenden Wellen, die sichin einer Lucke in der Wellenfront im Schatten des Storkorpers bemerkbar macht.Die durch den Mantel laufende Energie breitet sich dagegen ungestort aus.

Die vorliegende Untersuchung bietet damit eine Erklarung fur das Phanomen der Lg-Ausloschung in den Pyrenaen als Folge lokal verstarkter Heterogenitat und erhohterDampfung. Dieses Modell lasst sich unter Umstanden auch fur die Erklarung ahnlicherBeobachtungen in den Alpen heranziehen. In den Pyrenaen ist die betreffende Regiondurch besonders starke Deformation und durch Materialaustausch zwischen verschie-denen Tiefenbereichen aus Kruste und Mantel im Zuge der Auffaltung der Pyrenaengekennzeichnet. Sehr wahrscheinlich stehen diese Prozesse mit den veranderten Streu-und Dampfungseigenschaften in Zusammenhang.

118 SUMMARY

Phasen–Korrelationen

Im Gegensatz zur im ersten Teil verwendeten Energietransfertheorie, die die Phasen-beziehungen der gestreuten Wellenfelder vernachlassigt, stellen diese Phasenbezie-hungen die Informationsquelle der in Teil zwei verwendeten Methoden dar. Grundla-ge dieser Methoden ist die passive Gewinnung von Greenschen Funktionen zwischenzwei Punkten durch Korrelation der Aufzeichnungen eines diffusen Wellenfeldes andiesen Punkten. In der Seismologie hat dieses Konzept mit einer Veroffentlichungvon Campillo and Paul (2003) Einzug gehalten. Es grundet auf der Tatsache, dassPhasenbeziehungen von Wellen, die gleiche Wege zuruckgelegt haben, erhalten blei-ben, auch wenn sie von einer Vielzahl von Wellen uberlagert werden, die andereWege zuruckgelegt haben. Als diffuses Wellenfeld dient hier das seismische Umge-bungsrauschen.

Interferometrie rauschgenerierter Seismogramme am Merapi (Indonesi-en)

Hier wird eine Untersuchung beschrieben, die in Sens-Schonfelder and Wegler (2006a)veroffentlicht wurde. Sie enthalt einige wesentliche methodische Entwicklungen so-wie grundlegende Erkenntnisse. Die Untersuchung beginnt mit dem Nachweis, dassdie durch Korrelation von seismischem Rauschen gewonnenen Greenschen Funktio-nen nicht nur, wie bis dahin angenommen, die ballistischen Phasen enthalten, son-dern daruber hinaus auch die Codawellen, also indirekte, reflektierte und gestreuteWellen. Diese Erkenntnis eroffnet die Moglichkeit, die passiv gewonnenen Green-schen Funktionen in der Reflexionsseismik zu benutzen. Eine andere Moglichkeit,der hier nachgegangen wird, ist deren Verwendung fur hochgenaue interferometrischeMessungen von Veranderungen im Ausbreitungsmedium. Es wird gezeigt, dass sichaus eintagigen Registrierungen des seismischen Rauschens Greensche Funktionen re-konstruieren lassen, die sich mit dem an anderen Tagen aufgezeichneten Rauschenreproduzieren lassen. Darauf aufbauend wird die Interferometrie rauschgenerierterSeismogramme (Passive Image Interferometry) vorgestellt, die eine kontinuierlicheUberwachung der seismischen Geschwindigkeit im Untergrund ermoglicht. Sie nutzteine verbesserte Variante der als Codawellen-Interferometrie bekannten Technik zurMessung von Geschwindigkeitsanderungen.

Es werden Daten vom Vulkan Merapi im Zeitraum von August 1997 bis Juli 1999analysiert. Diese zeigen deutliche Schwankungen der seismischen Geschwindigkeit imBereich von wenigen Prozent. Die Variationen zeigen einen jahreszeitlichen Verlaufund hangen in ihrer Amplitude von der Laufzeit der seismischen Wellen ab. Bild 2Azeigt die gemessenen Geschwindigkeitsanderungen fur Laufzeiten zwischen 2-8 s.

Der Vorteil gegenuber bestehenden Methoden zur Uberwachung der seismischen Ge-schwindigkeit wird in dieser Abbildung deutlich. Die sechs Messpunkte der blauenund grunen Linien, die vergroßert dargestellt sind, stammen von aufwendigen ak-tiven Messungen mit wiederholbaren Quellen, die mit Codawellen-Interferometrieausgewertet wurden. Die kontinuierlichen Messungen mit der hier vorgestellten Me-thode wurden ohne zusatzlichen Feldeinsatz aus bestehenden Daten gewonnen.


−0.04

−0.03

−0.02

−0.01

0

0.01

0.02

0.03

0.04

dv/v Jun Jul Aug

−0.01

−0.005

0

Sep Nov 1998 Mar May Jul Sep Nov 1999 Mar May Jul

40

30

20

10

0

−0.02

−0.01

0

0.01

0.02

GW

L [m

]

δv

/v

0

30

60

prec

ipita

tion

[mm

/d]

A

B

Abbildung 2: Anderungen der seismischen Geschwindigkeit am Vulkan Merapi (A) undderen hydrologische Modellierung (B)

In einer nichtlinearen Inversion werden die Parameter eines hier entwickelten hy-drologischen Modells angepasst, das in der Lage ist, die gemessenen Geschwin-digkeitsanderungen durch Schwankungen im Grundwasserspiegel zu erklaren. DieVorwartsrechung dieser Inversion nutzt ortsabhangige Sensitivitatsfunktionen undbeschreibt neben dem Verlauf auch die Laufzeit–abhangigen Amplituden der Ge-schwindigkeitsvariationen. Die Anpassung des hydrologischen Modells an die Ge-schwindigkeitsvariationen wird in Abbildung 2B verdeutlicht. Im unteren Teil derAbbildung ist die tagliche Niederschlagsmenge in Blau und der daraus bestimm-te Grundwasserspiegel in Schwarz dargestellt. Im oberen Teil sind die gemessenen(Punkte) und vom Modell vorhergesagten (durchgezogene Linien) Geschwindigkeits-variationen fur zwei Laufzeitfenster dargestellt. Rot entspricht Laufzeiten zwischen2-4 s und Blau zwischen 6-8 s. Die gute Ubereinstimmung dieser Kurven belegt,dass das hydrologische Modell den korrekten physikalischen Zusammenhang zwi-schen Niederschlag und Anderungen der seismischen Geschwindigkeit beschreibt.Dies wiederum zeigt, dass die hier vorgestellte Methode tatsachlich in der Lage ist,geringfugige Anderungen der seismischen Geschwindigkeit zu uberwachen. Die Me-thode ist einfach, auf existierende Daten anwendbar und hat ein großes Potentialfur Anwendungen in anderen Bereichen (Snieder and Page, 2007).

120 SUMMARY

Tagliche Variationen und Instrumentenfehler

Hier werden zwei weitere Anwendungen der Methode zur Gewinnung von Green-schen Funktionen aus seismischem Umgebungsrauschen fur Veranderungsmessung-en vorgestellt. Die erste Untersuchung dieses Kapitels stellt eine Verbesserung derzeitlichen Auflosung der Interferometrie rauschgenerierter Signale vor. Sie beschreibtanhand der Daten vom Vulkan Merapi tageszeitliche Anderungen der seismischenGeschwindigkeit und deren spektrale Zusammensetzung. Der Vergleich mit den Am-plitudenspektren verschiedener meteorologischer Großen und deren Phasen deutetdarauf hin, dass die tageszeitlichen Schwankungen mit Anderungen der Lufttempe-ratur zusammenhangen.Eine andere Art von Veranderungen wird in der zweiten Studie mittels rauschge-nerierter Seismogramme untersucht. Hier wird der Effekt von Zeitfehlern in denSeismometern betrachtet. Es werden die Zeitdifferenzen zwischen den paarweise un-tersuchten Seismometern bestimmt. Eine Matrixinversion fuhrt anschließend zu denZeitfehlern der einzelnen Seismometer. Diese Methode bietet die Moglichkeit, dieGenauigkeit der Zeithaltung unabhangig von Zusatzinformationen direkt anhandder Daten zu uberwachen. Sie ist damit ideal geeignet, um die Zeithaltung einzelnerSeismometer auch im Nachhinein zu korrigieren.

Ausblick

Mit der Energietransfertheorie und der Nutzung von Phasenkorrelationen im seismi-schen Rauschen werden zwei Methoden verwendet, die in der Seismologie ein großesPotential fur weitere Anwendungen und neue Entwicklungen haben.Die Simulation der Energieausbreitung in einem kombinierten Modell aus statistischbeschriebenen kleinskaligen Heterogenitaten und deterministisch beschriebener ma-kroskopischer Struktur erlaubt weitergehende Untersuchungen von Problemen, indenen Diskontinuitaten eine Rolle spielen. Unter anderem ließe sich die Ausbreitungteleseismisch einfallender Energie in der Erdkruste studieren, die sehr empfindlichfur die Eigenschaften des Wellenleiters ist. Der nachste Schritt bei der Entwick-lung von Algorithmen zur Simulation der Energieausbreitung in der Erde sollte dieBehandlung von Oberflachenwellen und deren Kopplung mit Raumwellen sein.Die Interferometrie rauschgenerierter Seismogramme kann in verschiedensten, weituber die Seismologie hinausgehenden, Bereichen angewendet werden, bei denen keineaktiven Quellen eingesetzt werden konnen. Die Kombination dieser Methode mitLokalisierungstechniken fur die Quelle der Geschwindigkeitsvariationen scheint hiereine sinnvolle Entwicklung zu sein.

Literaturverzeichnis

Aki, K. and Richards, P. G. (1980). Quantitative Seismology, Theory and Methods,volume I. W. H. Freeman and Company, San Francisco.

Campillo, M. and Paul, A. (2003). Long-range correlations in the diffuse seismiccoda. Science, 299:547 – 549. doi:10.1126/science.1078551.

Ryzhik, L., Papanicolaou, G., and Keller, J. B. (1996). Transport equations forelastic and other waves in random media. Wave Motion, 24:327–370.

Sens-Schonfelder, C. and Wegler, U. (2006a). Passive image interferomety and sea-sonal variations of seismic velocities at Merapi volcano, Indonesia. Geophys. Res.Lett., 33:L21302. doi: 10.1029/2006GL027797.

Sens-Schonfelder, C. and Wegler, U. (2006b). Radiative transfer theory for estimati-on of the seismic moment. Geophys. J. Int., 167(3):1363–1372. doi: 10.1111/j.1365-246X.2006.03139.x.

Snieder, R. and Page, J. (2007). Multiple scattering in evolving media. Phys. Today,60(4):49–55.

121

Scientific developement

Christoph Sens-Schonfelder

10/1997 – 09/2003 Studies in geophysics at the Universitat Leipzig, (Leipzig,Germany)

05/2000 – 09/2003 Fellow of the German National Merit Foundation (Studien-stiftung des deutschen Volkes)

08/2000 – 06/2001 Studies in geophysics at the Haskoli Islands (Reykjavık,Iceland)

09/2003 Diploma

01/2004 – 07/2004 Scientific co-worker at the Universitat Kiel (Kiel,Germany)

since 08/2004 Scientific co-worker at the Universitı¿ 12

Leipzig (Leipzig,Germany)

11/2006 – 07/2007 Research visit at the Laboratoire de Geophysique Interneet Tectonophysique (Grenoble, France) with support of theGerman Academic Exchange Service (DAAD)

03/2007 Granting of the Gnter Bock award by the Deutsche Geo-physikalische Gesellschaft (German Geophysical Society)for an outstanding publication

Sens Schoenfelder Dissertation

Documents

Transcript of Sens Schoenfelder Dissertation