Megathrust Earthquakes: Study of Fault Slip and Stress Relaxation Using Satellite Gravity Observations

Megathrust Earthquakes: Study of Fault Slip and Stress Relaxation

Using Satellite Gravity Observations

Taco Broerse

Uitnodiging

voor de verdediging van het proefschrift

Megathrust Earthquakes:

Study of Fault Slip and

Stress Relaxation Using

Satellite

Gravity Observations

Taco Broerse

vrijdag 14 maart

12:00 lekenpraatje - senaatszaal Aula - 12:30 publieke verdediging - senaatszaal Aula - aansluitend receptie Aula

-Invitation

for the defense of the thesis

Megathrust Earthquakes:

Study of Fault Slip and

Stress Relaxation Using

Satellite

Gravity Observations

Friday March 14

12:00 introduction - senaatszaal Aula - 12:30 public defense - senaatszaal Aula - reception afterwards - Aula - 90˚E 100˚E 0˚ 10˚N
1 0 1 mm 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N
2
1 0 1 2 mm 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N
3
2
1 0 1 2 3 mm 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N
4
2 0 2 4 mm 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N
5 0 5 mm 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N
5 0 5 mm 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N
5 0 5 mm 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N
10 0 10 mm 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N 90˚E 100˚E 0˚ 10˚N

Megathrust

Earthquakes:

Study of Fault Slip

and Stress

Relaxation Using

Satellite Gravity

Observations

Taco Broerse

This thesis relates changes in the Earth’s gravity field

to mass displacements generated by very large

earthquakes. Large scale mass displacements take

place when the Earth’s crust and underlying ductile upper

mantle deform. These displacements are observed both

during earthquakes and in the years that follow, when

the Earth relaxes earthquake-induced stress.

Special attention is given to the contribution of ocean

mass redistribution to co-seismic gravity changes, and

its relation to changes in bathymetry.

Observations of earthquake-related gravity changes

are obtained from the GRACE twin-satellite mission that

produces monthly maps of the Earth's gravity field. Using

a combination of GRACE data and GPS observations this

thesis interprets ongoing gravity changes and crustal

motions after the 2004 Sumatra-Andaman earthquake

as dominantly caused by mantle creep. Contrasts

in relaxation styles from both observation types are

related to lateral variations in mantle rheology below the

subduction zone.

Megathrust Earthquakes: Study of Fault Slip and Stress Relaxation

Using Satellite Gravity Observations

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. ir. K.C.A.M. Luyben, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op vrijdag 14 maart 2014 om 12:30 uur

door

Dirk Bernard Taco BROERSE

Ingenieur luchtvaart en ruimtevaart geboren te Utrecht

Dit proefschrift is goedgekeurd door de promotor: Prof.ir. B.A.C. Ambrosius

Copromotor: Dr. L.L.A. Vermeersen

Samenstelling promotiecommissie:

Rector Magnificus, voorzitter

Prof.ir. B.A.C. Ambrosius, Technische Universiteit Delft, promotor

Dr. L.L.A. Vermeersen, Technische Universiteit Delft, copromotor

Prof.dr.-ing. habil. R. Klees, Technische Universiteit Delft

Prof.dr. W. Spakman, Universiteit Utrecht

Prof.dr. A. Hooper, University of Leeds

Dr. R.E.M. Riva, Technische Universiteit Delft

Dr. R. Govers, Universiteit Utrecht

Prof.dr. S.B. Kroonenberg, Technische Universiteit Delft, reservelid

Dr. R.E.M. Riva heeft als begeleider in belangrijke mate aan de totstandkoming van het proefschrift bijgedragen.

Dit proefschrift is financieel ondersteund door het Earth Research Centre Delft en ISES. De publicatie is mede mogelijk gemaakt door een financi¨ele bijdrage van de vakgroep Astrodynamica en Ruimtemissies, Faculteit Luchtvaart- en Ruimtevaarttechniek, Technische Universiteit Delft.

ISBN: 978-94-6186-282-2

Gedrukt door: CPI-W¨ohrmann Print Service Zutphen

Opmaak omslag: Janneke ten Kate

Contents

Summary v

Samenvatting viii

1 Introduction 1

1.1 The earthquake cycle . . . 1

1.2 Subduction zones . . . 1

1.3 Deformation of the Earth due to subduction earthquakes . . . 2

1.4 Flow and deformation within the solid Earth: rheology . . . 3

1.5 Observations of seismic deformation: GPS . . . 4

1.6 The role of the gravity field as an earthquake observation . . . 5

1.7 GRACE satellite mission . . . 7

1.8 Motivation for this study . . . 8

1.9 Dissertation outline . . . 9

1.9.1 Part 1: Co-seismic and post-seismic deformation and the use of satellite gravimetry . . . 9

1.9.2 Part 2: Co-seismic and post-seismic deformation and gravity mod-els: applications . . . 10

I

Co-seismic and post-seismic deformation: modeling

11

2 Normal modes relaxation models 12 2.1 Introduction . . . 12

2.1.1 Rheology of the mantle . . . 12

2.1.2 Normal modes . . . 15

2.2 Normal mode relaxation theory . . . 16

2.3 Using normal modes for co-seismic deformation and gravity changes . . . . 19

2.4 Burgers rheology in the normal mode context . . . 29

2.5 Relaxation modes . . . 29

2.6 Finding modes . . . 32

2.6.1 Extended root-finding strategy for compressible models . . . 33

2.6.2 Detected roots . . . 35

2.7 Post-seismic relaxation . . . 47

2.7.1 Contribution of individual modes . . . 48

2.7.2 Influence of model stratification on post-seismic deformation . . . . 54

3 A sea level equation for internal loading 62

3.1 Introduction to absolute and relative sea level changes . . . 62

3.1.1 The gravity field of the Earth and sea level . . . 62

3.2 A sea level equation for internal loading: co-seismic sea level changes . . . 65

3.3 The sea level equation deconstructed . . . 67

3.3.1 The global ocean assumption . . . 67

3.3.2 Signal loss due to truncation at low degree spherical harmonics . . 74

3.3.3 Loading, self-gravitation and ocean mass conservation . . . 75

3.4 Conclusions . . . 77

II

Co-seismic and post-seismic deformation and gravity models:

applications

78

4 Observing changes in the Earth’s gravity field due to earthquakes 79 4.1 Introduction . . . 79

4.2 Extracting earthquake signals from GRACE data . . . 83

4.3 Sensitivity of long wavelength gravity changes to fault slip . . . 84

4.4 Application to the 2004 December 26 Sumatra-Andaman earthquake . . . 92

4.4.1 Forward model results Sumatra-Andaman earthquake . . . 95

4.4.2 Co-seismic gravity signals from GRACE: Sumatra-Andaman . . . . 99

4.4.3 The effect of compressibility on co-seismic gravity changes . . . 109

4.4.4 The ocean contribution: discussion . . . 116

4.4.5 Discussion GRACE and forward model Sumatra-Andaman earthquake117 4.5 Application to the 2010 Febraury 27 Maule earthquake . . . 119

4.5.1 Forward model results Maule earthquake . . . 119

4.5.2 Co-seismic gravity signals from GRACE: Maule . . . 120

4.5.3 Discussion GRACE and forward model Maule earthquake . . . 122

4.6 Application to the 2011 March 11 Tohoku earthquake . . . 129

4.6.1 Forward model results Tohoku earthquake . . . 129

4.6.2 Co-seismic gravity signals from GRACE: Tohoku . . . 131

4.6.3 Discussion GRACE and forward models Tohoku earthquake . . . . 141

4.7 Conclusions and Summary . . . 145

4.7.1 Sensitivity of the gravity field to shallow earthquakes . . . 145

4.7.2 Ocean contribution to co-seismic gravity field changes . . . 146

4.7.3 Solid Earth compressibility and co-seismic gravity changes . . . 147

4.7.4 The contribution of GRACE to the study of earthquakes sources . 148 4.7.5 Differentiation between co-seismic and post-seismic gravity changes 149 4.8 Co-seismic deformation and satellite gravimetry: outlook . . . 150

5 Mantle relaxation after the 2004 Sumatra-Andaman earthquake 151 5.1 Introduction . . . 151

5.1.1 Macro-rheology vs. micro-rheology . . . 152

5.1.2 Lateral heterogeneities in subduction zones . . . 153

5.2 GPS observed displacements . . . 156

5.2.1 Processing . . . 156

5.2.2 Analysis of post-seismic displacements . . . 157

5.3 Viscoelastic relaxation model . . . 158

5.3.1 Completeness of the normal mode solutions . . . 160

5.3.2 Gravity modeling . . . 161

5.4 Comparison GPS observations and model . . . 161

5.4.1 Horizontal deformation . . . 162

5.4.2 Vertical deformation . . . 166

5.5 GRACE . . . 173

5.6 Comparison GRACE and model . . . 176

5.7 Discussion . . . 180

5.7.1 Observations: lateral contrast in viscosity . . . 180

5.7.2 Depth of the low viscosity mantle wedge . . . 181

5.8 Conclusions . . . 182

5.9 Outlook . . . 184

Appendices 185

A Redundancy of negative orders in seismic forcing 186

B Modes and residuals 188

C Ocean contribution to post-seismic gravity changes 189

References 195

About the author 208

List of publications 209

Summary

Large earthquakes in subduction zones pose a big threat to communities close to faults,

as recently encountered after the Mw ∼ 9 2004 Sumatra-Andaman, 2010 Maule (Chile)

and 2011 Tohoku (Japan) events. On the other hand these events provide scientists a natural laboratory for the study of earthquake mechanisms, tectonics and large scale Earth deformation in general.

Plate-tectonics comprises the movement of several thin and rigid plates on top of a hot, weak and ductile upper mantle. Deformation of the plates concentrates along faults that may both be plate boundaries as well as internal fractures. One of the main features of tectonics is subduction, which involves the sinking of slabs of oceanic lithosphere deep into the mantle at convergent plate boundaries. Subduction creates converging plate motions at opposing sides of the plate boundaries. When friction forces do not allow stable sliding of the plate interfaces, strain accumulates in the crust surrounding a fault. During an earthquake, strain that may have built up for centuries is (partly) released, resulting in seismic waves and so-called co-seismic deformation. An earthquake causes uplift and subsidence combined with horizontal displacement of crust in a wide area around a rupture.

In the days to decades after the main shock ongoing deformation is usually observed, which we describe as post-seismic deformation, and that is related to relaxation of earthquake induced stresses. Main contributors are thought to be ductile, viscoelastic, creep of the upper mantle and a slow continuation of the rupture (afterslip). Co-seismic as well as post-seismic deformation in subduction zones are the topics of this thesis.

Observations of earthquake-induced deformation

For studying earthquake related deformation I rely on two satellite techniques: firstly, GPS measurements from sites around subduction zones that directly measure displacement of the Earth’s crust; secondly, the monthly maps of gravity field changes from the Gravity Recovery and Climate Experiment (GRACE) satellite mission. GRACE relies on the simple but beautiful concept of a pair of proof masses (the satellites themself) in a very low Earth orbit that are sensitive to very small variations in the Earth’s gravity field. By measuring changes in their 200 km inter-satellite distance, large scale movement of mass over the Earth’s surface and in its interior can be inferred on a monthly basis (chapter 1). In this thesis I relate changes in the Earth’s gravity field with mass displacements due to deformation of the Earth’s crust and underlying ductile upper mantle. The GRACE twin-satellite mission maps the Earth’s gravity field since March 2002 and as such has been able to observe gravity changes due to the three mentioned megathrust earthquakes. While GRACE provides a coarse spatial resolution (features of 40,000 to approximately 400 km can be detected) the scale of these earthquakes and the massive amount of fault slip is sufficient to produce gravity field changes to be detected.

As mentioned before, I use earthquake-related changes in the gravity field to study two distinct, but related phenomena: the near-instantaneous change (co-seismic) in the gravity field during an earthquake that reflects slip of the fault interface in the Earth’s crust; and the gradual change in the gravity field in the months to years after a large seismic event (post-seismic), which indicates slow relaxation of earthquake-stresses in the Earth’s interior. For the latter topic I also combine GRACE data with GPS observations to study lateral changes in the properties of the ductile upper mantle around subduction zones. I interpret co-seismic and post-seismic gravity changes and deformation based on forward models, where I use a semi-analytical normal modes model to compute surface deformation and changes in the gravitational potential due to seismic forcing in a spherically layered, compressible and self-gravitating Earth.

Co-seismic gravity changes and the effect of oceans

Beginning with co-seismic gravity changes: investigating the distribution of co-seismic slip on a subduction fault is of great importance for understanding the rupture process and for determination of seismic coupling of a subduction zone (i.e. what portion of the tectonic convergence leads to seismic slip). However, the use of GRACE gravity data for estimating co-seismic slip is still in its infancy. The main reason for this lies in the complex influence of co-seismic uplift and subsidence of the ocean floor on water redistribution, and its effect on gravity change. Part of this dissertation is devoted to the influence of changes in the ocean floor topography on seismic gravity changes and how this affects the interpretation of fault slip as inferred from GRACE data.

In continental areas the patterns of co-seismic change in gravity are similar to the patterns of co-seismic uplift and subsidence as these largely coincide with locations of gravity increase and decrease (section 2.3). In this thesis I show that for sub-oceanic earthquakes this symmetry is broken due to water redistribution, especially in the long wavelength gravity field as observed by GRACE (chapter 4). A largely positive gravity change is changed into a dominantly negative gravity change due to the ocean’s presence. I show that the main driver is a broad uplift of the ocean floor due to co-seismic expansion of the Earth, which effectively replaces a large amount of water mass (section 4.4.3). Furthermore, horizontal crustal movement in the direction of the trench can reduce the water column height when the ocean floor slopes towards the trench. I demonstrate that for the 2011 Tohoku earthquake, when slip is assumed to peak close to the surface as follows from tsunami studies, its effect on gravity is very large: peak amplitude of the gravity anomalies increases by 40% at GRACE resolution (section 4.6.2).

Studies published in the last couple of years that estimated the seismic moment (which can be read as an integral of fault slip over the rupture area) of a few recent megathrust earthquakes based on GRACE data generally used simple models to account for the water redistribution. In this dissertation I show that these simple models fail to properly account for the co-seismic ocean effects when faults rupture close to continental areas. Subsequently, GRACE data is interpreted incorrectly in the case of the 2010 Maule earthquake and to

a lesser extent in case of the 2011 Tohoku earthquake (chapter 3). This dissertation

interpretation of seismic gravity changes, allowing for a proper inclusion of GRACE data in geodetic or seismic studies of seismic slip. I apply the developed models to the case of the 2004 Sumatra-Andaman, 2010 Maule and 2011 Tohoku earthquakes and compare results with co-seismic gravity changes obtained from GRACE. Using slip distributions derived from independent data, I indicate what GRACE derived gravity changes can contribute to these earthquake models in terms of updates of main rupture depth and dip angles (chapter 4). Post-seismic gravity changes and surface deformation: indications of lateral viscosity contrasts around subduction zones

Observations of ongoing post-seismic gravity changes and deformation offer an opportunity for studying mantle deformation in subduction zones. The amplitude, spatial pattern and temporal behavior of the post-seismic trend in gravity or deformation can be explained as a result of mantle flow in reaction to earthquake stresses. In the case of earthquake-induced solid Earth deformation, crustal motions and gravity changes immediately after a loading event can be observed. This means that the initial transient deformation can be examined, which may be dominant up to a few years after an earthquake and that produces a faster deformation rate than what is observed in the long-term. Following a number of studies that interpret transient deformation using a bi-viscous Burgers rheology I adopt this rheology for modeling post-seismic relaxation (chapter 2).

More than seven years of observations of post-seismic relaxation after the 2004 Sumatra-Andaman earthquake provide an improving view at the mechanisms that deform the wide vicinity of the rupture. I use both GRACE gravity field data and GPS observations and interpret post-seismic changes as dominantly coming from mantle creep (chapter 5). With increasing time GPS and GRACE show contrasts in relaxation styles that could not easily be seen using shorter time series. Namely, the GPS time series show large contrasts be-tween initial displacement rates (which are fast) versus velocities at the end of the seven year observation period, while GRACE shows a much more gradual change of post-seismic gravity changes during this period. I find the depth extent of the currently ongoing mantle creep using the spherically layered normal modes relaxation models and a Burgers rheology combined with vertical GPS deformation and horizontal displacements at far-field sites. However, models that fit GPS data well do not match the much slower increase in post-seismic gravity changes as observed by GRACE. As GRACE is observing the largest changes over oceanic area, while the GPS observations measure continental displacements, I suggest that the contrasts seen in relaxation behavior are due to lateral changes in the properties of the asthenosphere below the subduction zone. Finally, I argue that viscoelastic relax-ation models using 1-D viscosity profiles cannot reproduce GPS and GRACE observrelax-ations of post-seismic relaxation simultaneously.

Samenvatting

Megathrust aardbevingen: Een studie naar breukverschuiving en postseismische ontspanning met gebruik van satellietzwaartekrachtwaarnemingen

Zware aardbevingen vormen een levensgrote bedreiging voor bewoners van gebieden nabij

breukzones, zoals opnieuw duidelijk is geworden na de Mw ∼ 9 Sumatra-Andaman (2004),

Maule (2010, Chili) en Tohoku (2011, Japan) aardbevingen. Aan de andere kant bieden deze allesvernietigende natuurverschijnselen wetenschappers juist een natuurlijk laboratorium voor onderzoek naar aardbevingsmechanismen, plaattektoniek of meer in het algemeen het de-formatiegedrag van de aarde.

Plaattektoniek beschrijft de beweging van een dunne en rigide lithosfeer, gefragmenteerd in verschillende platen, die drijft bovenop een hete, zwakke en plastisch vervormbare asthe-nosfeer. De lithosfeer bevat zowel de korst als het topje van de mantel en de asthenosfeer vormt een op langere tijdschalen zwakke laag van de bovenmantel. Vervorming van de lithosfeer vindt met name plaats langs bestaande breuken die zowel aan de plaatranden, als binnenin platen te vinden zijn. Een van de meest in het oog springende kenmerken van plaat-tektoniek is subductie. Hierbij zinkt het uiteinde van een oceanische plaat diep weg in de onderliggende mantel, daar waar twee platen elkaar ontmoeten. Belangrijk hierbij is dat sub-ductie zorgt voor convergente bewegingen aan beide zijden van de plaatgrenzen. Wanneer geleidelijk over elkaar schuiven van de betrokken platen door wrijving op de breukvlakken niet mogelijk is, komt de korst rond een dergelijke breuk onder toenemende spanning te staan. Deze spanning kan zich gedurende eeuwen opbouwen voordat een aardbeving ontstaat en de spanning weer (gedeeltelijk) ontladen wordt, daarbij resulterend in wat we coseismische ver-vorming noemen. Een groot gebied rondom een breuk kan ofwel omhoog geduwd worden, dan wel plotseling dalen waarbij de aardkorst eveneens in horizontale beweging verschuift. Grootschalige postseismische vervorming wordt meestal waargenomen in de dagen tot de-cennia na de hoofdschok en dit gedrag hangt samen met het ontspannen van spanningen die nieuw opgewekt zijn door de aardbeving. Hierbij spelen visco-elastische kruip van de asthenosfeer en een geleidelijke en aseismische (stille) voortzetting van de aardbeving (in het Engels afterslip genoemd) vermoedelijk de hoofdrol. Coseismische en postseismische vervorming vormen de hoofdthema’s van dit proefschrift.

Observaties van aardbevinggerelateerde korstvervorming

Bij het bestuderen van door aardbevingen veroorzaakte vervorming van de aarde maak ik intensief gebruik van twee satelliettechnieken: allereerst, GPS-metingen uitgevoerd in de omgeving van subductiezones, waarbij rechtstreeks verplaatsing van de aardkorst gemeten wordt; ten tweede, observaties van de Gravity Recovery and Climate Experiment (GRACE) satellietmissie die op maandelijkse basis veranderingen in het aardse zwaartekrachtveld in

kaart brengt. GRACE berust op het simpele doch elegante principe van een tweetal proef-massa’s (in de vorm van twee satellieten) in een zeer lage aardbaan, die letterlijk de speel-bal zijn van zeer kleine variaties in de aantrekkingskracht van de aarde. Door variaties op de 200 km afstand tussen de tweelingsatellieten te meten door middel van een micro-golfverbinding worden grootschalige massaverplaatsingen op en onder het aardoppervlak in kaart gebracht. In dit proefschrift breng ik zwaartekrachtsveranderingen in verband met massaverplaatsing door vervorming van de aardkorst en onderliggende laagvisceuze

asthe-nosfeer. De GRACE-tweelingsatellietmissie is in bedrijf sinds maart 2002 en was zo in

staat zwaartekrachtveranderingen veroorzaakt door de drie eerder genoemde aardbevingen waar te nemen. Ondanks GRACE’ grove ruimtelijke resolutie (patronen met afmetingen tussen de 40.000 en 400 km kunnen worden waargenomen) waren de schaal en aanzienlijke verschuiving tijdens deze aardbevingen genoeg om voor GRACE waarneembare zwaarte-krachtsveranderingen te veroorzaken.

Zwaartekrachtsveranderingen gebruik ik om de twee eerder genoemde verschijnselen te be-studeren: de vrijwel instantane veranderingen (coseismisch) in zwaartekracht die duiden op verschuiving van de platen op het breukvlak, binnenin de aardkorst; en de langzame ver-anderingen in zwaartekracht in de maanden tot jaren na een grote schok (postseismisch), welke wijzen op een ontspanning van spanningen die zijn veroorzaakt door de aardbeving. Voor de postseismische veranderingen maak ik naast de GRACE-metingen gebruik van GPS-tijdseries, om met deze combinatie laterale veranderingen (parallel aan het aardoppervlak) in vervormingsgedrag van de ondiepe mantel in subductiezones te herleiden.

Voor het verklaren van coseismische en postseismische veranderingen in zwaartekracht en vervorming van de aarde baseer ik mij op voorwaartse modellen, gebruikmakend van het semi-analytische normal mode (staande golf-) model om vervorming van het aardoppervlak en veranderingen in de zwaartekrachtspotentiaal te berekenen, ten gevolge van een seis-mische belasting binnen een radiaal symmetrische, compressiebele en zichzelfaantrekkende aarde.

Coseismische zwaartekrachtsveranderingen en het effect van oceanen

Om met coseismische zwaartekrachtsveranderingen te beginnen: het is van groot belang voor het begrijpen van aardbevingsmechanismen en het in kaart brengen van seismische koppeling (d.w.z. welk gedeelte van de tektonische beweging aan de plaatranden daad-werkelijk door middel van aardbevingen ontladen wordt) om de reikwijdte en de ruimtelijke

verdeling van de verschuiving op breukvlakken te bestuderen. Echter, het gebruik van

GRACE-gegevens voor dit doel staat nog in de kinderschoenen. De belangrijkste reden die hiervoor kan worden aangewezen is de bijdrage van zeebodemdaling en -stijging ten gevolge van aardbevingen aan waterverplaatsingen, en de wijze waarop dit doorwerkt in zwaarte-krachtsveranderingen. De studie hiervan beslaat een groot gedeelte van dit proefschrift waarbij ik bekijk hoe de aanwezigheid van oceanen rondom subductiezones de interpretatie van GRACE-gegevens be¨ınvloedt.

In de afwezigheid van zee¨en zijn de ruimtelijke patronen in coseismische zwaartekrachtsver-anderingen vergelijkbaar met die van de coseismische bodemdaling en -stijging, aangezien

deze grotendeels samenvallen met gebieden met negatieve dan wel positieve zwaartekrachts-veranderingen (paragraaf 2.3). Later laat ik echter zien dat waterverplaatsing ten gevolge van een aardbeving deze symmetrie grotendeels teniet doet, wanneer, zoals het geval is voor alle bestudeerde aardbevingen, deze plaatsvindt in de nabijheid van de zee, of onder het zeeoppervlak. Dit geldt nog het meest voor de lange golflengtes van het zwaartekrachtveld, het deel waar GRACE het meest gevoelig voor is (hoofdstuk 4). Voor subductiegerelateerde aardbevingen verandert de aanwezigheid van een oceaan de coseismische zwaartekrachts-verandering van een hoofdzakelijke toename tot een overwegende afname. Zoals ik in para-graaf 4.4.3 aantoon, hangt deze afname samen met een uitgestrekte expansie van de aarde, resulterend in het wegstuwen van een grote hoeveelheid oceaanwater. Daarnaast gaat de horizontale verplaatsing van de zeebodem in de buurt van de trog een grote rol spelen indien de zeebodem een sterke helling kent in de richting van de trog, waarbij de bovenliggende waterkolom wederom verder verkleint. In het geval van de Tokoku-aardbeving is deze invloed op de zwaartekrachtsverandering aanzienlijk indien wordt aangenomen dat de beweging op de breuk vlakbij het oppervlak piekte zoals uit tsunami-studies blijkt: uit het model volgt

een toename van 40% van de piek in zwaartekrachtanomali¨en op de GRACE-meetresolutie

(paragraaf 4.6.2).

In de afgelopen paar jaar zijn enkele studies verschenen die op basis van GRACE-gegevens het seismisch moment schatten van enkele recente subductiezone-aardbevingen, waarbij de auteurs merendeels eenvoudige modellen gebruikten om de rol van de oceaan te verwerken. Het seismisch moment kan hier gelezen worden als een sommatie van verschuiving op de breuk, genomen over het oppervlak waar het verzet plaatsvond. In dit proefschrift laat ik zien dat deze modellen niet in staat zijn de coseismische oceaanbijdrage aan zwaartekracht te benaderen wanneer een aardbeving nabij een continent plaatsvindt. Daardoor wordt de GRACE-data onjuist ge¨ınterpreteerd, met name in het geval van de Maule-aardbeving en

in zekere mate in het geval van de Tohoku-aardbeving (hoofdstuk 3). Dit proefschrift

biedt de gereedschappen (in de vorm van een seismische zeespiegelvergelijking) voor een verbeterde uitleg van seismische zwaartekrachtsveranderingen, waardoor het gebruik van GRACE-gegevens binnen geodetische of seismische studies naar aardbevingsmechanismen mogelijk wordt. Ik pas de ontwikkelde modellen toe op de Sumatra-Andaman-aardbeving (2004), de Maule-aardbeving (2010) en de Tohoku-aardbeving (2011) en vergelijk resul-taten met de coseismische zwaartekrachtsverandering zoals ik die uit de GRACE-data win. Gebruikmakend van bestaande modellen van breukverschuiving op basis van onafhanke-lijke metingen, geef ik aan hoe coseismische zwaartekrachtsveranderingen verkregen met GRACE een andere optiek levert op de diepte van de verschuivingen en op de diphoek, de hoek waaronder de breuk met het oppervlak staat (hoofdstuk 4).

Postseismische zwaartekrachtsveranderingen en vervorming van het aardoppervlak: aanwijzingen voor laterale viscositeitsverschillen rondom subductiezones

Waarnemingen van postseismische korstvervorming en zwaartekrachtsveranderingen bieden een mogelijkheid om de reologie van de mantel (het kruipgedrag) onder subductiezones te bestuderen. De amplitude, het ruimtelijke patroon en de tijdsafhankelijkheid van de post-seismische trend in zwaartekracht of korstvervorming (in GPS-waarnemingen laat deze trend

een snel beginverloop zien waarna deze in de loop van enkele jaren convergeert naar een lagere snelheid) kunnen herleid worden tot mantelstromingen, ontlokt door spanningsveran-deringen ten gevolge van de breukverschuivingen tijdens de laatste aardbeving. In het geval van vervorming van de aarde na aardbevingen, kunnen vervormingen en zwaartekrachtsver-anderingen kort na een belastingsverandering waargenomen worden. Dit houdt in dat de aanvankelijke, kortstondige (in het Engels transient) vervorming kan worden waargenomen, het type vervorming dat de overhand heeft tot enkele jaren na een aardbeving en die tot snellere ontspanning leidt dan de vervorming die wordt vernomen over de langere termijn. In de duiding van deze kortstondige vervorming volg ik recente studies die een verklaring zoeken in een Burgers reologie, die twee verschillende viscositeiten in zich verenigt die elk op verschillende tijdschalen werken (hoofdstuk 2).

Een beter inzicht in de mechanismen achter postseismische ontspanning, zichtbaar in de wijde omtrek van een breuk, is mogelijk door de inmiddels zeven jaar aan postseismische waarnemingen sinds de Sumatra-Andaman aardbeving van 2004. Hierbij maak ik gebruik van zowel GRACE-data als GPS-waarnemingen van korstverplaatsing in Thailand en Maleisi¨e. Hierbij beschouw ik mantelontspanning als belangrijkste aandrijver van post-seismische ver-anderingen (hoofdstuk 5). Met het verstrijken van de tijd blijken GPS- en GRACE-tijdseries contrasterend gedrag in ontspanningssnelheid te vertonen. Hierbij is in de GPS-tijdseries een groot verschil te zien tussen de aanvankelijke verplaatsingssnelheid en de langzamere verplaatsingen aan het eind van de zeven jaar van de meetserie. Aan de andere kant laten de tijdseries verkregen met de GRACE-gegevens een veel constanter verloop zien in de postseismische zwaartekrachtsveranderingen gedurende dezelfde tijdspanne. Op basis van de radie¨el symmetrische normal mode modellen vind ik met de vertikale GPS-verplaatsingen en de horizontale GPS-verplaatsingen van de meest ver weg gelegen meetstations een grote gevoeligheid voor de diepte van de postseismische mantelontspanning. Aan de andere kant blijken de modellen die de GPS-data in het back-arc-gebied goed beschrijven geen verklaring te bieden voor het veel langzamere verloop van de postseismische zwaartekrachtsverander-ingen. Aangezien ik uit GRACE-data de grootste postseismische veranderingen vind boven oceanisch gebied en GPS continentale verplaatsingen meet, stel ik dat de verschillen in ontspanningsgedrag te herleiden zijn tot laterale viscositeitsverschillen in de asthenosfeer onder de subductiezone. Hierbij betoog ik ook dat de veelgebruikte visco-elastische mo-dellen die aannemen dat viscositeit voornamelijk in radie¨ele richting verandert niet geschikt zijn om gelijktijdig GPS- en GRACE-waarnemingen uit subductiezones te verklaren.

1

Introduction

1.1

The earthquake cycle

Earthquakes can be extremely destructive natural phenomena causing many deaths, destroy-ing infrastructure and deprivdestroy-ing large numbers of people from their homes. On the other hand earthquakes give scientists the opportunity to study the response of the solid Earth to tectonic loading. Ironically, often the largest earthquakes provide the best natural labora-tories for these studies, due to their far-reaching spatial influence. In this thesis I study the large scale crustal deformation that is caused by earthquakes, both the near-instantaneous deformation during a seismic event as well as the gradually increasing deformation that fol-lows in the days to centuries afterwards. Earthquakes should be understood in the framework of an earthquake cycle: when adjacent plates move relative to one another, friction often prohibits movement at the fault interfaces. Consequently stress is building up at the plate boundaries, and once a frictional threshold has been overcome the fault slips during an earthquake. This releases (part of) the strain that has accumulated for decades to cen-turies. Afterwards the fault is locked again and the earthquake cycle repeats in a similar - but not equal - fashion. Three stages of deformation are distinguished: inter-seismic deformation (the long-term accumulation of strain), co-seimic deformation (the release of strain during an earthquake) and post-seismic deformation (the response of the Earth due to co-seismic stress changes); in this thesis I examine the last two types. Note that while earthquakes show cyclic behavior, a rupture is never a repetition of a previous earthquake nor have fault zone fixed earthquake recurrence intervals (e.g. Heimpel (1997); Satake and Atwater (2007); Scholz and Campos (2012)).

1.2

Subduction zones

I focus on so-called megathrust earthquakes in subduction zones, which can be found along continental margins and produce the largest earthquakes; bear in mind the recent Mw ∼ 9

2004 Sumatra-Andaman, 2010 Maule (Chile) and 2011 Tohoku (Japan) earthquakes. The theory of plate-tectonics prescribes a thin and rigid outer layer of the Earth called litho-sphere (which includes both the crust and the uppermost rigid part of the mantle) that is fragmented into several plates with different motions, floating on top of a weak and ductile upper mantle, called asthenosphere. Deformation of the lithosphere is accommo-dated by weak fractured zones, especially at plate boundaries. Subduction, the process of lithospheric slabs slowly dipping into the hot and ductile underlying mantle, occurs at the boundaries of tectonic plates with converging motions (figure 1.1). Subduction always applies to oceanic lithosphere, which due to its high density is able to slide under another continental (or oceanic) plate and sink into the mantle due to gravity. This not only recy-cles the oceanic crust while it plunges into the mantle and slowly dissolves, but also drives global tectonic motions, even though how exactly subducting slabs control plate motions

is still heavily debated (e.g. Conrad and Lithgow-Bertelloni (2002); Becker and Faccenna (2009)). While tectonic plate movement over long time scales is a continuous process, at the surface friction between the two plate interfaces often inhibits continuous slip between two plates. Earthquakes usually nucleate above the ductile-brittle transition depth, a depth that appears to be controlled by either temperature or slab intersection with the fore-arc mantle (Hyndman et al., 1997) and that determines the dominant mode of deformation. The seismogenic zone where stick-slip earthquake behavior determines the intraplate mo-tion is located above this transimo-tion, while below this transimo-tion the increased temperature

(above 350◦ C, e.g. Hyndman and Wang (1993)) allows episodical slow slip events and

further down stable sliding at plate tectonic rates (e.g. Scholz and Campos (2012)).

1   viscous   asthenosphere   sea  level   locked  fault   volcanic  islands   magma  

Figure 1.1: Subduction zone at a convergent ocean-continent plate boundary. The top layer con-sists of an elastic crust and elastic upper mantle (lithosphere) above the ductile (visco-elastic) part of the upper mantle (asthenosphere). The oceanic elastic slab subducts below the overriding continental lithosphere and sinks into the mantle. Fault locking and subsequent seismic slip occurs in the elastic part of the continental crust and lithosphere, while stable sliding takes place below the ductile-brittle transition due to increased temperature with depth. Around the transition it-self episodical slow slip events take place (e.g. Scholz and Campos (2012)). While the oceanic slab descends into the hot and low viscous asthenosphere, sediments and water are transported along, changing the composition and strength of the mantle wedge while material upwelling causes vulcanism (e.g. van Keken (2003)).

1.3

Deformation of the Earth due to subduction

earthquakes

While the plate motions are typically less than 10 cm/y, on the timescales of decades and centuries the integrated slip deficit at subduction zones can lead to very large earthquakes and on far longer timescales these motions can lead to complete separation of continents.

Earthquake related deformation may result in uplift or subsidence of areas, a clearly visible example being the emergence of coral reefs above sea level off the coast of Sumatra in 2004 (Meltzner et al., 2006). Horizontal displacement may be less visible from a layman’s point of view but is well observable with techniques such as GPS; after the 2011 Tohoku earthquake, using sea floor reference points and a combination of GPS and acoustic techniques up to 24 meter of horizontal displacement has been measured offshore Japan (Sato et al., 2011). The horizontal motions of the crust, for tectonic plates as a whole or within plates, are indicators of ongoing tectonic processes and indicate the accumulation of strain in the period between earthquakes. The release of energy and deformation of the crust during an earthquake, the so-called co-seismic deformations, are a direct result of the inter-seismic stress accumulation (and deformation) around fault zones. Using deformation models and measured co-seismic deformation next to related observations of seismic waveforms, the distribution of slip during a rupture can be estimated using inversion techniques. By comparing this slip distribution to observations of inter-seismic crustal motion, the physics of fault zones and seismic risks can be better understood.

While earthquakes co-seismically release decade or century long stress accumulation, at the same time the rupture introduces new stresses on the parts of the fault that did not slip co-seismically, as well as in the underlying upper mantle. Aseismic slip (slip too slow to be recorded seismically) can occur for more than a year after an earthquake (Pritchard and Simons, 2006) and is accompanied by numerous smaller aftershocks (e.g. Perfettini and Avouac (2007)). Afterslip is thought to take place on parts of the fault that due to their frictional behavior did not allow rupture during the actual earthquake but that experienced large co-seismic stress perturbations (Marone et al., 1991). Secondly, the ductile mantle that is found below the brittle lithosphere experiences stress changes as well and slowly starts to flow as to relax the newly introduced co-seismic stresses (e.g. Thatcher and Rundle (1984)). Depending on the viscosity of the mantle this post-seismic relaxation process may take years to centuries (Wang et al., 2012a). While it has been known for several decades that both afterslip and visco-elastic mantle relaxation likely occur simultaneously, with the advance of dense GPS networks around subduction fault zones and with the launch of satellite gravity measurements next to the occurrence of several large Mw ∼ 9 earthquakes,

post-seismic deformation can now be studied in much more detail.

1.4

Flow and deformation within the solid Earth:

rheology

When discussing rheology, the study of flow and deformation of matter, it is almost unavoid-able not to put forward Heraclitus’ notion of everything being in constant flux, repeated by several classical philosophers such as Plato (in his work Cratylus) and Simplicius as ”παντ α ρι- everything flows”. The phrase rheology was coined by Eugene C. Bingham in the early 20th century after this Greek aphorism. This statement of constant change has even a broader meaning in the field of geodynamics as subduction establishes a huge conveyor belt constantly creating new oceanic lithosphere at the mid-oceanic ridges and recycling sedi-ments, water and oceanic lithosphere as a slab descends into the hot mantle. At the same

time this process slowly, but insistently, changes the outlook of the world by creating topo-graphy, changing the arrangement of continents and changing climate, causing vulcanism and redirecting ocean currents.

Apart from the elastic (immediate and recoverable) deformation, for example occurring during earthquakes, rheology focuses on the ductile flow and viscous (gradual and irrecov-erable) deformation of the solid Earth under loading conditions. While there are several approaches to rheology, for example the microphysical laboratory study of mantle materials under high pressure and temperature to mimic mantle conditions, in this work I focus on the macro scale of geodetic observations in the natural laboratory of post-seismic viscoelastic relaxation. As the forcing (the earthquake) is relatively well known, from observations of post-seismic deformation we may obtain information of the viscosity of the upper mantle and time dependency of the rheology in the region of the fault. The time after which a Newtonian1_{material changes from dominantly elastic (solid) to dominantly viscous}

deforma-tion (fluid) is characterized by τ = η_µ (e.g. Karato (2010)), where τ is the relaxation time of the material, η is the viscosity and µ the shear modulus, which quantifies the material stiffness. Geodetic estimates for the upper mantle viscosity differ many orders of magni-tude: 5· 1017_{− 10}21_{Pa s where the lower values have been found for plate boundaries and}

the high end values apply to continental interiors, see for a review B¨urgmann and Dresen (2008). Depending on location, the timescales at which viscoelastic deformation acts thus differs by a few orders of magnitude as well. In contrast to the steady-state viscosity values derived from geodetic observations, from laboratory rock mechanical studies it follows that viscosity under some conditions is a stress-dependent function (viscosity decreases when stress increases), which implies that the time scale of viscoelastic relaxation of a load be-comes a function of the load itself (e.g. Ranalli (1995)). Knowledge on mantle rheology is both essential in understanding long term dynamics of the Earth as well as understanding the several phases of deformation of the earthquake cycle.

1.5

Observations of seismic deformation: GPS

A substantial part of this type of research is aided by observations of earthquake related surface deformation that can be made using the Global Positioning System commonly known as GPS. The use of geodetic quality GPS receivers at either permanent GPS stations or routinely visited GPS markers, combined with advanced data processing allows to precisely (in the order of a few mm, e.g. Simons et al. (2007)) record co-seismic displacements as well as post-seismic movements at GPS station locations in a wide region around a rupture. While post-seismic deformation is a slow process, the continuous mode of measurements

allows tracking of deformation that would otherwise carry on unnoticed. The Mw 7.3

1992 Landers earthquake, California, was one of the first seismic events after which GPS observations were used to infer post-seismic crustal motions (Ivins, 1996; Savage and Svarc, 1997), and since then GPS has become the most important method for observing post-seismic deformation.

1.6

The role of the gravity field as an earthquake

observation

Until now I have described co-seismic and post-seismic relaxation in terms of deformation, which can be directly observed as surface motions. A second important source of infor-mation that I use in this work is the gravity field. The Earth’s gravity field is to a high degree homogeneous, while relatively small latitudinal differences are caused by the Earth’s equatorial flattening and the centrifugal force. However next to these variations, slight local deviations in gravity disclose local mass variations inside the Earth and differences in topography or bathymetry.

As large earthquakes result in massive movements of crust and mantle, the resulting uplifted or subsided areas imply redistribution of mass and thus changes in the gravity field. The gravity field can be described by various observables, which all are functions of the gravity potential W . First, the gravitational potential V is defined as the density ρ of an attracting mass divided by the distance to the observer l , integrated over the volume v of the attracting mass and times the gravitational constant G:

V = G

Z Z Z _ρ

ld v (1.1)

To include the centrifugal force, the gravity potential W is defined as: W (x , y , z ) = V (x , y , z ) +1

2ω

2_d2

z (1.2)

where ω is the angular velocity and dzis the distance to the rotation axis z . The gravitational

attraction is the gradient of the potential

g =|∇W | (1.3)

In this thesis I use two functionals derived from the gravity potential: geoid height and gravity anomalies. The geoid is that equipotential surface of the gravity potential W (there is an infinite set of equipotential surfaces) that best fits an idealized global ocean surface at rest. Being an equipotential surface, at the geoid holds for the gravitational potential: W = constant = W0. As the Earth’s shape resembles, due to its equatorial flattening,

an ellipsoid of revolution, the height of the geoid surface is referenced to an equipotential surface U of an ellipsoid. This ellipsoid is defined such that it optimally fits the geoid and that it has the same mass as the Earth. At this equipotential surface holds U = U0= W0.

Subsequently, the geoid height N is then the difference between the geoid surface and the reference ellipsoid. As the shape of the geoid is quite smooth and the reference ellipsoid closely approximates the geoid, maximal geoid undulations - geoid height differences - are in the order of 100 meters referenced with respect to the reference ellipsoid. Figure 1.2 shows a global field of the static geoid height.

As equipotential surface the local gravitational attraction is everywhere perpendicular to the geoid, and as such the geoid is the zero elevation surface for the oceans if they were at rest

Figure 1.2: Global static geoid heights (in meters) derived from the GRACE Gravity Model 03 (Tapley et al., 2007). Geoid height has a range of approximately -85 to 105 meter, with respect to the GRS 1980 reference ellipsoid. Peaks in geoid height indicate areas with increased mass concentrations, while geoid lows show mass deficits with respect to the global average. As the geoid closely follows the ocean surface when dynamic effects are removed, also the sea level shows the same peaks and valleys as shown in this figure. The static geoid height variations are largely due to density variations inside the Earth, while time variations are dominated by mass redistribution over the Earth’s surface, i.e. variations in ice sheet dynamics, hydrology of large river systems and ocean circulation. However, part of the time variable signal can be explained by solid Earth deformation phenomena such as ongoing uplift of regions that hosted large ice sheets during glacial maxima or earthquake related deformation around subduction zones.

and undisturbed by tides, currents or atmosphere. The geoid thus coincides with absolute sea level (as measured from the ellipsoid). For an ocean at rest the geoid is the lowest possible level given the amount of water in the oceans, even while the height of the geoid globally (smoothly) varies. As a representation of gravity, the geoid is sensitive to lateral inhomogeneities of density inside the Earth and changes in topography and bathymetry. Next to absolute sea level, I consider seismic changes in relative sea level. Relative sea level S, the difference between absolute sea level and sea floor topography, is the sea level as experienced from the Earth’s surface and is very much affected by subsidence and uplift in coastal areas. Relative sea level change will be further discussed in chapter 3.

The other expression of gravity consists of gravity anomalies ∆g that denote a difference between the observed gravity acceleration (on the geoid) and the normal gravity reference value on the ellipsoid potential γ. Normal gravity is the gradient of the normal potential U, compare equation 1.3 for the gravitational acceleration. Here Q is a point on the ellipsoid and P is a projection from the ellipsoid onto geoid (Hofmann-Wellenhof and Moritz, 2006):

∆g = gP− γQ (1.4)

To be able to describe the small deviations in the Earth’s gravity field, gravity anomalies are usually expressed in the unit milligal, which equals to 1· 10−5_m/s2_{. Throughout this}

thesis gravity anomalies are often used as these give larger weights to smaller wavelength signals, compared to geoid height, which proves to be useful for smaller scale features in GRACE gravity fields.

For studying earthquake related deformation I focus on small gravitational changes in time (in the order of a few cm geoid height), for which I rely on satellite gravimetry. The gravity field can change due to seismic deformation in a number of ways:

• Changes in topography: uplift means a local increase in mass, hence an increase in local strength of the gravity field.

• Expansion (dilatation) inside the Earth: a local decrease in density causes a local decrease in mass and thus gravity field.

• Uplift of the ocean floor: water mass is pushed away, implying a local decrease of the gravity field.

1.7

GRACE satellite mission

Changes in observed gravity can thus serve as an indirect measure of seismic deformation. Since March 2002 the Gravity Recovery and Climate Experiment (GRACE) satellite mission maps the Earth’s gravity field in a monthly fashion, which allows an unprecedented view of global mass changes over the Earth. The GRACE mission consists of two satellites that follow the same orbit, but separated approximately 220 km along-track, and can observe monthly gravity changes with spatial resolutions of 400 to 40,000 km wavelength (Tapley et al., 2004). Since satellite orbits are largely determined by the gravitational attraction of the Earth (and a number of other smaller and well-known forces), any irregularity in the Earth’s gravity field results in slight changes of the satellite’s motion. Because of their low orbit (∼500 km) the GRACE satellites are sensitive to small variations in the gravity field. Due to the satellites’ along-track separation, the same variations in the gravity field are experienced with a small time delay (figure 1.3). Precise K-band microwave inter-satellite ranging observes the satellite separation and is combined with other instrument information and complex modeling to compute monthly gravity field solutions by a few data centers. From small changes in satellite motions we can thus indirectly observe large scale mass variations in the Earth and over the Earth’s surface caused by earthquakes; provided that the earthquakes are very large: currently five events with a moment magnitude Mw 8.5 or

GRACE  A   GRACE  B   GPS  satellite   GPS  antenna   upli3   incoming  GPS  signal  

Figure 1.3: Schematic overview of the observation types used in this thesis for the study of earthquake related deformation. A broad area that shows uplift implies a local increase in mass, hence an increase in gravity. When one of the low altitude GRACE satellites passes over the area with increased mass, gravity has slightly increased at satellite altitude, introducing minute changes in satellite acceleration. As GRACE A is 200 km ahead of GRACE B, both satellites visit the same locations with a small time delay. Consequently, any irregularities in the gravity field and the resulting accelerations, will be experienced by GRACE A prior to GRACE B. When a mass change at the Earth surface, or inside the Earth, is sufficiently large, accelerations lead to changes in inter-satellite distance that can be picked up by the precise microwave ranging of the twin-inter-satellites. The figure indicates that when the satellites pass over an area with a positive mass anomaly (e.g. due to uplift) GRACE B accelerates when it approaches the anomaly, while GRACE A is decelerated again as the attraction of the mass anomaly has become opposite to its orbital motion. Next to observing deformation indirectly by gravity measurements, displacements can be directly observed by GPS measurements. Conducting continuous or episodical GPS measurements at fixed sites, changes in position in time of that point can be observed.

1.8

Motivation for this study

As I have written previously, in this dissertation I study two earthquake-related phenom-ena: co-seismic deformation and post-seismic deformation. Both deformation types are part of the earthquake cycle, but are distinct in the timeframes in which they act: episodic (co-seismic) versus continuous (post-seismic). For the co-seismic case I study how

satel-lite measurements (especially by the GRACE mission) can be used to observe the seismic displacement at fault interfaces. I focus on earthquakes occurring along subduction zones where often oceans cover the ruptured area. Several studies have focused on using GRACE data for observing co-seismic changes in the gravity field, or even on inferring seismic source parameters from GRACE gravity fields. However, especially at the start of this research, the role of ocean water to changes in the ocean floor topography, while being a first order effect, was insufficiently understood and is still commonly modeled improperly in literature (see chapters 3 and 4). This has severe implications for the usefulness of GRACE observa-tions for the study of large earthquakes. Therefore, a large part of the study of co-seismic gravity changes is devoted to understanding the, sometimes counterintuitive, changes in the gravity field due to ocean water movements after large earthquakes. Only when the contri-bution of the ocean to changes in the gravity field is well understood, satellite gravimetry can finally contribute to studies of the slip in subduction zones.

The other main topic of this study is the origin of post-seismic deformation observed af-ter large subduction earthquakes. While this dissertation is not the first to study this phenomenon, the use of gravity observations and the occurrence of several very large earthquakes since 2004 provides new possibilities of studying the mechanism behind the post-seismic relaxation. Due to the large magnitude of the 2004 Sumatra-Andaman earth-quake seven years of continuous and episodical GPS measurements in the wide area around the rupture allow to constrain the origin of the mechanism that causes post-seismic sur-face deformation and gravity field changes. To the best of my knowledge, up to date the publications that investigated the post-seismic relaxation after the 2004 event never used more than four years of GPS data, or authors used longer periods but did not combine GPS with GRACE gravity data (section 5.1.2). This means that a study combining the longer time span of data of both GPS and GRACE can better constrain the upper mantle rheology than previously. As I will show chapter 5, GPS observed motions and post-seismic GRACE gravity field changes show spatial variations in relaxation behavior that may provide an added value of space gravimetry to study of subduction zone rheology.

1.9

Dissertation outline

This dissertation is divided in two parts. The first part encompasses the theoretical mod-els of solid Earth deformation and gravity changes and ocean water redistribution due to

earthquakes. The second part consists of applications of the models to recent Mw ∼9

earthquakes.

1.9.1

Part 1: Co-seismic and post-seismic deformation and the use

of satellite gravimetry

In the following chapter 2 I introduce the normal mode viscoelastic deformation models. In the existing normal modes model I have introduced a Burgers rheology that allows for a transient deformation phase prior to steady-state relaxation of stress. The chapter fur-thermore discusses and tests a number of approximations on which I later build in the part

2 of this thesis. Next, chapter 3 introduces the sea level equation that I have adapted to the seismic case. Using the case of the 2010 Maule earthquake I show the essential im-provements of the sea level equation for modeling the co-seismic gravity changes compared to commonly used approximate methods of modeling co-seismic gravity field changes over oceanic areas.

1.9.2

Part 2: Co-seismic and post-seismic deformation and gravity

models: applications

Chapter 4 starts with a discussion of sensitivity of the gravity field changes to fault slip. It continues in applying the developed models of part 1 to three recent megathrust earth-quakes: the 2004 Sumatra-Andaman, 2010 Maule and 2011 Tohoku megathrust events. Here I compare models (with input of land based earthquake observations) of co-seismic gravity change with actual GRACE measured gravity changes. I explain the contribution of both ocean water redistribution and solid Earth dilatation to the observed long wave-length gravity changes. As there are significant differences between model and GRACE derived changes, I indicate how GRACE complements slip models using land-based (GPS and seismological) observations. The final chapter 5 descends into post-seismic territory and investigates the transient deformation in southeast Asia measured with GPS and ob-served by GRACE. I study upper mantle viscosities and depth extent of mantle relaxation using the widespread deformation patterns as measured with GPS. Even more so, here I compare differences in relaxation as measured by GPS and GRACE and make a link with lateral variations in subduction zone rheology, thereby concluding this research.

Part I

Co-seismic and post-seismic

deformation: modeling

2

Normal modes relaxation models

2.1

Introduction

The viscoelastic nature of the Earth’s interior can be seen in various geodetic or geophysi-cal observations of crustal deformations, changes in sea-level or in the gravity field due to loading of the Earth. Loading or unloading can be due to past events such as the retreat of the large ice sheets of the last glacial cycle or the loading of the interior of the Earth due to tectonic processes such as earthquakes. The response of the Earth to these loadings has an immediate, elastic, nature which is followed by a slow, often exponentially decaying, viscoelastic response (see e.g. B¨urgmann and Dresen (2008)). Viscoelastic creep shows that the mantle is essentially a fluid, be it with a high viscosity causing its response to be very slow. Tectonically induced loading, such as by earthquakes that often occur inside the elastic upper part of the Earth, causes differential stresses inside the underlying viscous mantle. The resulting mantle creep gradually transfers stresses back into the crust caus-ing displacements that can be observed at the Earth’s surface (e.g. Freed et al. (2010)). Thousands of years may be needed for the mantle to fully relax after a loading event; for example in Fennoscandia vertical deformation up to 10 mm/yr is still observed (Haskell, 1935; Johansson et al., 2002), caused by the melt of glaciers at the end of the last glacial period roughly 10,000 years ago, a process which is referred to as glacial isostatic adjust-ment. Indicators of changes in coastlines can show us the amount and rate of past crustal movements (Farrell and Clark, 1976). Geodetic observation techniques such as GPS help to measure either ongoing deformation of the crust due to recent processes, such as post-seismic deformation, or the current crustal deformation that is still occurring due to the unloading of formerly glaciated areas since the late Pleistocene.

2.1.1

Rheology of the mantle

Modeling deformation of the Earth due to loading combined with actual observations of crustal deformation is instrumental in understanding the physics of flow properties of the Earth’s mantle. Rheology, the study of flow, or creep, of materials, deals with the relation between rate of deformation, strain rate, and stress. The simplest linear relations that serve well to describe stress and strain rate relations as observed in materials are Hookean, New-ton, Maxwell, Kelvin and Burgers bodies. These mathematical formalisms have analogies in the form of combinations of springs and dashpots. Elastic, instantaneous deformation can be described using a linear relation between stress and strain (Hooke solid) and can be visualized by a spring. A Hooke solid has the property that strain is recoverable; as soon as a load is removed the solid instantaneously returns to its original shape. This relation explains well the rheology of the crust and upper mantle (lithosphere) as long as strains are small enough and no brittle failure occurs (faults) (e.g. B¨urgmann and Dresen (2008)). At higher temperatures rocks are viscous on very long timescales, meaning they can flow

!"

_#"

$

#

!!

a

!

_"

##

$

_"

!!

!#

_%#

$

%

!!

b

Figure 2.1: Schematic representation of a) a Maxwell rheology, consisting of a spring and dashpot in series, b) Burgers rheology, consisting of a Maxwell element in series with a Kelvin element where a dashpot and spring are connected in parallel. Viscosity is denoted by η and rigidity (or shear modulus) by µ.

(strain rate) when the Earth is loaded. Material having a linear relation between stress and strain rate is called a Newtonean fluid, which can be schematically represented by a dashpot. However, in intermediate times rocks show both elastic and viscous behavior and can be described as a combination springs and dashpots.

A Maxwell rheology is the most common viscoelastic rheology used in literature, and can be schematically represented by a spring with rigidity µm that represents the immediate

elastic response in series with a the spring with the dashpot with viscosity ηm that stands

for the viscoelastic response that slowly relaxes stresses in the Earth’s mantle (figure 2.1). Especially in the case of interpreting crustal deformation due to glacial isostatic adjustment where we measure a response long after the load changes (withdrawal of ice sheets) a Maxwell rheology seems to fit the observations of steady deformation rates rather well (e.g. Milne et al. (2001); Paulson et al. (2007); van der Wal et al. (2011)).

In the case of a Kelvin rheology the spring with rigidity µk and dashpot with viscosity ηk are

coupled in parallel, which results in delayed elastic strain that is recoverable after the load is removed. Both Maxwell and Kelvin models have characteristic times defined as τ = η/µ, denoting the timescales for which the mechanical behavior changes. A Maxwell body will behave like an elastic solid for timescales t << τm and like a viscous material for t >> τm;

a Kelvin body on the other hand will behave like an elastic material for t > τk (Karato,

2010).

The simplest model that incorporates both elastic, transient creep and steady state creep is the bi-viscous Burgers rheology (Peltier et al., 1981). The Burgers model essentially extends the Maxwell model by taking a Maxwell element and coupling this in series with a Kelvin element (Burgers, 1939). Here the Kelvin element schematically consists of a spring coupled in parallel with a dashpot (see figure 2.1). Here the Kelvin viscosity ηk is usually

between Kelvin and Maxwell rigidity γ is typically smaller than 0.5 (Hetland and Hager, 2006; Hoechner et al., 2011). As the newly introduced spring is coupled in parallel by a dashpot it does not introduce any elastic deformation, but it allows for an anelastic transient viscoelastic deformation that causes a fast relaxation quickly after the application of a load (such as the occurrence of an earthquake). The influence of the transient deformation then diminishes, since it allows for only limited strain as the effective viscosity increases with time (Rumpker and Wolff, 1996), after which the Maxwell relaxation becomes dominant at longer timescales.

In contrast to the long term crustal movement seen from glacial isostatic adjustment studies, post-seismic deformation from recent earthquakes can be studied as a short term Earth response. Observations of post-seismic deformation can be done by geodetic techniques such as GPS or gravity satellite missions such as GRACE, which can detect large scale mass movements related to solid Earth deformation. GPS observations of areas where large earthquakes occurred such as in south-east Asia after the Sumatra-Andaman earthquake of 2004 indicate a transient crustal movement that slows down quickly in the time span of a few years (Gahalaut et al., 2008; Panet et al., 2010). This fast initial crustal motion is difficult to reproduce using viscoelastic models incorporating Maxwell rheologies (Paul et al., 2012) but can be well fitted using a Burgers rheology (Pollitz et al., 2006; Hoechner et al., 2011). Transient deformation after the 2004 Sumatra-Andaman earthquake is discussed in detail in chapter 5.

Next to empirical, geodetic studies, laboratory rock mechanical experiments enable a dif-ferent perspective on mantle rheology. Difdif-ferent modes of rock creep are studied from a microphysical perspective, where these experiments are conducted under high pressure and high temperatures. These studies show two main creep processes: diffusion creep and dislo-cation creep which involve diffusion and movement of vacancies in crystal lattices and along grain boundaries. Experimental data has been fitted through constitutive equations where strain is shown to be dependent on a number of variables such as stress, grain size, water content and temperature (e.g. Hirth and Kohlstedt (2003)). Dislocation creep exhibits a non-linear strain behavior due to a power-law dependence on stress, which results in a decrease of effective viscosity (which is stress divided by strain rate) when stress levels rise (e.g. B¨urgmann and Dresen (2008)).

Next to the previously mentioned empirical, geodetic studies, also laboratory experiments suggest a transient phase in mechanical behavior of several mantle materials. In microcreep tests of olivine polycrystals three phases of deformation are observed: elastic (represented by the spring of the Maxwell element), recoverable transient creep (represented by the Kelvin element) and permanent viscoelastic deformation (represented by the Maxwell element) (Faul and Jackson, 2005). Because of these strong indications of importance of transient as well as steady state viscosity, both from the lab as well as space geodetic observations of post-seismic deformation, I will use the Burgers rheology in my models of solid Earth deformation.

2.1.2

Normal modes

Modeling solid Earth deformation can consist of numerical simulations but also of analyti-cal models of viscoelastic deformation. In this chapter I rely on normal mode models that analytically allow to determine the Earth’s responses due to loading (e.g. Sabadini and Ver-meersen (2004)). Normal modes represent standing waves (or free oscillations) in terms of Earth deformation that arise due to instantaneous loading and are dependent on the internal structure of the Earth. The concept of free oscillation of a sphere is an important concept in both seismology (as the Earth’s normal modes can be observed after large earth-quakes) as well as the static deformation of a spherical Earth. Here static deformation can be modeled as an oscillation with an infinite period (Ben-Menahem and Singh, 1981). For example, one of the first realistic estimates for the seismic moment of the 2004 Sumatra-Andaman earthquake was based on a seismological analysis of normal modes that could be observed for weeks after the earthquake (Stein and Okal, 2005; Shearer and B¨urgmann, 2010). On the other hand, normal mode relaxation models have been used for a many decades to model viscoelastic deformation of the solid Earth, with important applications such as glacial isostatic adjustment (e.g. Peltier (1974)) or post-seismic deformation (e.g. Piersanti et al. (1995)).

Normal mode models assume a spherical symmetry of the Earth and a linear viscosity, which only varies with depth. These models can capture the large scale features of the inner Earth based on a phenomenological approach to mantle rheology. That is, rheology is seen from a perspective obtaining consistency with observations of deformation while using the simplest possible viscoelastic model. Here is not accounted for the physical mechanisms which explain time-dependent deformation at microscopic scale. A feature of normal modes models is that the total viscoelastic response to loading is a summation of a number of individual contributions (modes) that are clearly linked to the way the Earth model is stratified in spherical layers. Each layer boundary leads to one or more so-called relaxation modes that represent the physical interaction between the layers. This means that surface displacements can be seen as a summation of discrete contributions which in turn can be traced back to the structure of the chosen Earth model. This characteristic will be used later in this chapter.

While a number of authors use a Burgers rheology to explain the transient deformation as seen in observations of crustal deformation (Yuen and Peltier, 1982; Sabadini et al., 1985; Rumpker and Wolff, 1996; Pollitz, 2003; Hetland and Hager, 2006; Pollitz et al., 2006; Melini et al., 2008; Cannelli et al., 2010; Spada et al., 2011; Hoechner et al., 2011; Hu and Wang, 2012), none showed the complete set of resulting modes in detail nor described the influence of the individual modes resulting from a Burgers rheology. In this theoretical chapter will be shown how to implement a Burgers rheology in the normal mode approach. Furthermore a discussion of the relaxation modes is presented which illustrates the behavior of a Burgers body. Because for a compressible Earth model successfully finding modes for each subsequent degree is nontrivial due to the presence of an infinite amount of dilatation modes (Han and Wahr, 1995) I also describe a method that performs this task.

The goal of this chapter is to implement a compressible Burgers rheology in the normal mode frame work and study its application in modeling solid Earth deformation due to slip on finite faults. Furthermore, a main topic of the chapter is to investigate the sensitivity for co-seismic and post-seismic deformation and gravity changes for several model parameters, especially the level of radial stratification. Results of this analysis will prove to be important in the second part of this thesis where models will be applied to several cases of co-seismic and post-seismic deformation (chapter 5).

2.2

Normal mode relaxation theory

In the normal modes approach the elastic and viscoelastic response of the Earth due to a load are analytically solved by the momentum equation (which guarantees the conser-vation of momentum) and the Poisson equation (which describes the perturbation of the gravitational potential due to deformation), combined with conservation of mass and the constitutive equations that provide a stress-strain relation. Those equations are solved for a global Earth in the Laplace-transformed domain, such as first used by Peltier (1974). Here the Earth is represented by a spherically layered, self-gravitating model. A transfor-mation of the equations to the Laplace domain and assumption of linearity allows us to apply the Correspondence Principle to determine the viscoelastic response using the same set of equations that is used to solve the purely elastic equations (Peltier, 1974). The dis-placement vector u = (ur, uθ, uφ) in spherical coordinates (with r radius, θ colatitude and φ

longitude) is decomposed in two parts: a spheroidal (or poloidal) and a toroidal part (e.g. Smylie (1965); Ben-Menahem and Singh (1981); Piersanti et al. (1995)). Here spheroidal deformation is associated with vertical mass transport and toroidal deformation corresponds to vortex-like motions (Schubert et al., 2001).

u = us+ ut =∇ × ∇ × [ ˜S(r )· er] +∇ × [ ˜T (r )· er] (2.1)

Here ˜S(r ) denotes a spheroidal scalar function, ˜T (r ) a toroidal scalar function,∇× the curl of a vector and er the radial unit vector. Furthermore, all variables are expanded using the

spherical harmonic function Ym

l (θ, φ), e.g. ˜ T (r ) = ∞ X l =0 l X m=−l ˜ T_lm(r )Y_lm(θ, φ) (2.2) Y_lm(θ, φ) = P_lmc os(θ)ei mφ (2.3)

with associated Legendre polynomials P_lm of degree l and order m. When expressing the spheroidal part of equation 2.1 in spherical harmonics and expanding the expression, the surface deformation (Piersanti et al., 1995) and the gravitational potential perturbation φ1

(Sabadini and Vermeersen, 2004) read:     ur uθ uφ φ1     = ∞ X l =0 l X m=−l     U_lm Vm l ∇θ Vm l ∇φ −φm l     Y_lm(θ, φ) (2.4)