Computational Seismology: An Introduction Li Zhao Institute of - - PowerPoint PPT Presentation

Computational Seismology: An Introduction

Li Zhao

Institute of Earth Sciences Academia Sinica, Taipei 11529, Taiwan

e-Science Application Workshop 2011.03.19 ISGC 2011

Pathways of Seismology Study

source Earth instrument

instrument response

Observation Modeling

Before ¡earthquake ¡ A.er ¡earthquake ¡

2008 Wenchuan, China, Earthquake

Rupture duration: ~100 s Rupture length: ~300 km Wang et al. (2008) 2008.05.12 ¡ ¡14:28:04-­‑14:29:23

Rupture of Wenchuan Earthquake

March 11, 2011, M9.0 Japan Earthquake

Rupture duration: ~180 s Rupture length: ~500 km Hayes, USGS (2011)

Major Factors Influencing Ground Motion

Amplification by soft near-surface material Earthquake source rupture directivity Focusing of energy by material boundaries

Komatitsch et al. (2005) P wave (fast, weak) S waves (slower, stronger) Surface waves (slowest, strongest) Snapshot of Denali, Alaska, Earthquake

Surface Waves Cause Most Damages

University ¡of ¡Tokyo ¡(1990)

Amplification by Soft Sediments

Sedimentary basin Basement rock

Soft sedimentary basin can amplify ground motion

1994 Northridge (CA) Earthquake

Pacific Ocean Los Angeles Basin

A A B B Basin geometry focusing can further amplify ground motion.

Click to show movie

Directivity Effect on Ground Motion

Reliable Ground Motion Prediction Requires

Reliable model of the source Computational tools Earth structure model

Observation

Deploy instruments to collect waveform records

Waveforms carry information about the sources of earthquake and the structure along the paths between source and station

Paradigm of Seismology Problems

Prediction

Calculate waveforms for models

• f source and structure

Issue: efficiency and accuracy

Measurement

Compare predicted and recorded waveforms

Differences (residuals, anomalies) serve as data to refine source and structure models.

Set Up Inverse Problem

Relate data and model perturbations

Issue: (1) linear relation (2) both source and structure models

Inversion

Solve the inverse problem

Model Update

Collection of Global Earthquake Data

Observations vs. Predictions

Yang, Zhao & Hung (2010)

A major task of computational seismology is to solve the wave equation.

[ : ( ) ] ρ − ∇ ⋅ ∇ = u C u f 

Wave Equation: Realistic Earth model: 3D, irregular geometry, topography finite-difference method (FDM) finite-element method -- spectral-element method (SEM) Simplified Earth model: 1D, simple geometry Semi-analytic method: normal modes High-frequency approximation: ray theory (optics)

Equation of motion (2nd Newton’s Law) for continuum:

ρ = ∇ ⋅ + u τ f 

{ }

T 1 2

: : [( ) ( ) ] : ( ) = = ∇ + ∇ = ∇ τ C ε C u u C u

Hooke’s Law:

density body force (earthquake) stress acceleration elasticity tensor

Theoretical Foundation of Seismology

Finite-Difference Method (FDM)

(1) Discretize the medium and time by a (usually uniform) spatial-temporal grid. (2) Approximate the differential equation by finite difference schemes.

ρ = ∇ ⋅ + u τ f 

{ }

T 1 2

: : [ ( ) ( ) ] : ( ) = = ∇ + ∇ = ∇ τ C ε C u u C u

Simple Finite-Difference Schemes

From the Taylor expansion for function Φ(x): Forward- and backward difference formula for 1st-order derivative: Central-difference formula for 1st-order derivative: Central-difference formular for 2nd-order derivative:

) ( )] ( ) ( [ 1 ) (

2

h O x h x h x

x

+ − + = ∂ Φ Φ Φ ) ( )] ( ) ( [ 1 ) (

2

h O h x x h x

x

+ − − = ∂ Φ Φ Φ ) ( ) ( 6 1 ) ( 2 1 ) ( ) ( ) (

4 3 3 2 2

h O x h x h x h x h x

x x x

+ ∂ + ∂ + ∂ + = + Φ Φ Φ Φ Φ ) ( ) ( 6 1 ) ( 2 1 ) ( ) ( ) (

4 3 3 2 2

h O x h x h x h x h x

x x x

+ ∂ − ∂ + ∂ − = − Φ Φ Φ Φ Φ ) ( )] ( ) ( [ 2 1 ) (

3

h O h x h x h x

x

+ − − + = ∂ Φ Φ Φ ) ( )] ( ) ( 2 ) ( [ 1 ) (

4 2 2

h O h x x h x h x

x

+ − + − + = ∂ Φ Φ Φ Φ

1st-order accuracy 2nd-order accuracy 3rd-order accuracy

Higher-Order Finite-Difference Schemes

) ( ) ( 3 1 ) ( 2 ) ( ) (

5 3 3

h O x h x h h x h x

x x

+ ∂ + ∂ = − − + Φ Φ Φ Φ ) ( ) ( 9 ) ( 6 ) 3 ( ) 3 (

5 3 3

h O x h x h h x h x

x x

+ ∂ + ∂ = − − + Φ Φ Φ Φ ) ( ) ( 48 )] 3 ( ) 3 ( [ )] ( ) ( [ 27

5

h O x h h x h x h x h x

x

+ ∂ = − − + − − − + Φ Φ Φ Φ Φ ) ( )] 3 ( ) 3 ( [ 48 1 )] ( ) ( [ 16 9 ) (

5

h O h x h x h h x h x h x

x

+ − − + − − − + = ∂ Φ Φ Φ Φ Φ

From the Taylor expansions: and the fourth-order finite-difference scheme:

) ( ) ( 24 1 ) ( 6 1 ) ( 2 1 ) ( ) ( ) (

5 4 4 3 3 2 2

h O x h x h x h x h x h x

x x x x

+ ∂ + ∂ + ∂ + ∂ + = + Φ Φ Φ Φ Φ Φ ) ( ) ( 24 1 ) ( 6 1 ) ( 2 1 ) ( ) ( ) (

5 4 4 3 3 2 2

h O x h x h x h x h x h x

x x x x

− ∂ + ∂ − ∂ + ∂ − = − Φ Φ Φ Φ Φ Φ ) ( ) ( 24 64 ) ( 6 27 ) ( 2 9 ) ( 3 ) ( ) 3 (

5 4 4 3 3 2 2

h O x h x h x h x h x h x

x x x x

+ ∂ + ∂ + ∂ + ∂ + = + Φ Φ Φ Φ Φ Φ ) ( ) ( 24 64 ) ( 6 27 ) ( 2 9 ) ( 3 ) ( ) 3 (

5 4 4 3 3 2 2

h O x h x h x h x h x h x

x x x x

− ∂ + ∂ − ∂ + ∂ − = − Φ Φ Φ Φ Φ Φ

We obtain these relations:

), ( ρ

2

τ ∂ = ∂

x t u

u k

x

∂ = τ v k

x t

∂ = τ ∂ , ρ τ ∂ = ∂

x tv

, 2 2 ρ

1 1 1 1

x t v v

n i n i n i n i i

Δ τ − τ = Δ −

− + − +

x v v k t

n i n i i n i n i

Δ − = Δ τ − τ

− + − +

2 2

1 1 1 1

1D wave equations (2nd order):

Finite-Difference Implementation: 1D Case

Introducing velocity: Applying the central-difference scheme:

u v

t

∂ =

) ( ρ 1

1 1 1 1 n i n i i n i n i

x t v v

− + − +

τ − τ Δ Δ + = ) (

1 1 1 1 n i n i i n i n i

v v x t k

− + − +

− Δ Δ + τ = τ

and we have two 1st-order equations: We get the 2nd-order finite-difference equations:

) ( )] ( ) ( [ 2 1 ) (

3

h O h x h x h x

x

+ − − + = ∂ Φ Φ Φ

Proc 000 Proc 007 Proc 001 Proc 008 Proc 127

East (X) North (Y)

Parallelization (128 Processes)

Up (Z)

) ( ρ 1

1 1 1 1 n i n i i n i n i

x t v v

− + − +

τ − τ Δ Δ + = ) (

1 1 1 1 n i n i i n i n i

v v x t k

− + − +

− Δ Δ + τ = τ

x

n

x

n 2 1

x

n 3

x x

n N 8 =

y y

n N 4 =

z z

n N 4 =

y

n

z

n

Finite-Difference Method (FDM)

(1) Discretize the medium and time by a (usually uniform) spatial-temporal grid. (2) Approximate the differential equation by finite difference schemes. Only neighboring points are

• linked. Easy implementation of domain-

decomposition for distributed parallel programs. (3) FD grid is usually uniform for easy implementation and bookkeeping (advantage), but can not handle irregular geometry (shortcoming). Recent advances enable irregular grid.

ρ = ∇ ⋅ + u τ f 

{ }

T 1 2

: : [ ( ) ( ) ] : ( ) = = ∇ + ∇ = ∇ τ C ε C u u C u

(1) Discretization of the medium by a volumetric mesh: a collection of (irregular) hexahedral cells (elements). (2) Transformation from hexahedral cells to unit cubes. (3) Integrate the wave equation in space to obtain its weak form. (4) Expand all functions (displacement, stress, model, etc.) in spectral basis functions in each cell.

M k B c u

i k N i i k k

,..., 2 , 1 ), ( ) (

1

= − = ∑

=

x x x

F c G = ⋅

M cells (elements) with N control points in each cell

(-1,-1,-1) (1,1,1)

Spectral-Element Method (SEM)

3 3 3

d d d

⊕ ⊕ ⊕

ρ = ∇ ⋅ +

∫ ∫ ∫

u r τ r f r 

{ }

T 1 2

: : [ ( ) ( ) ] : ( ) = = ∇ + ∇ = ∇ τ C ε C u u C u

ρ = ∇ ⋅ + u τ f  

(5) SEM uses higher-order basis function; and the coefficient matrix G is

• diagonal. Easy for inversion and parallelization.

Examples of Wave Propagation Simulations

Caltech’s Southern California ShakeMovie Portal

E-mail Dissemination of SoCal ShakeMovie

Click to show movie

PGV Map and Waveform Fitting

Princeton’s Global ShakeMovie Portal

Global ShakeMovie

Tromp et al. (2011)

Click to show movie

Shake Movie: Chi-Chi (Taiwan) Earthquake

Click to show movie

Near-surface velocity

Thank you!