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 Observation source Earth instrument instrument response Modeling

2008 Wenchuan, China, Earthquake Before ¡earthquake ¡ A.er ¡earthquake ¡

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

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 Focusing of energy by material boundaries Earthquake source rupture directivity

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

Amplification by Soft Sediments Sedimentary basin Basement rock University ¡of ¡Tokyo ¡(1990) � Soft sedimentary basin can amplify ground motion

1994 Northridge (CA) Earthquake A B A B Los Angeles Basin Pacific Ocean Basin geometry focusing can further amplify ground motion.

Directivity Effect on Ground Motion Click to show movie

Reliable Ground Motion Prediction Requires Earth structure model Reliable model of the source Computational tools

Paradigm of Seismology Problems 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 Prediction Measurement Calculate waveforms for models Compare predicted and of source and structure recorded waveforms Issue: efficiency and accuracy Differences (residuals, anomalies) serve as data to refine source and structure models. Model Update Set Up Inverse Problem Relate data and model perturbations Inversion Issue: (1) linear relation (2) both source and Solve the inverse problem structure models

Collection of Global Earthquake Data

Observations vs. Predictions Yang, Zhao & Hung (2010)

Theoretical Foundation of Seismology Equation of motion (2 nd Newton’s Law) for continuum:  density u τ f body force (earthquake) ρ = ∇ ⋅ + elasticity tensor stress acceleration { T } Hooke’s Law: τ C : ε C : 1 [( u ) ( u ) ] C : ( u ) = = ∇ + ∇ = ∇ 2  Wave u [ C : ( u ) ] f ρ − ∇ ⋅ ∇ = Equation : Simplified Earth model: 1D, simple geometry Semi-analytic method: normal modes High-frequency approximation: ray theory (optics) Realistic Earth model: 3D, irregular geometry, topography finite-difference method (FDM) finite-element method -- spectral-element method (SEM) A major task of computational seismology is to solve the wave equation.

Finite-Difference Method (FDM) { } T  τ C : ε C : [ ( u ) ( u ) ] C : ( u ) u τ f 1 = = ∇ + ∇ = ∇ ρ = ∇ ⋅ + 2 (1) Discretize the medium and time by a (usually uniform) spatial-temporal grid. (2) Approximate the differential equation by finite difference schemes.

Simple Finite-Difference Schemes From the Taylor expansion for function Φ (x): 1 1 2 2 3 3 4 ( x h ) ( x ) h ( x ) h ( x ) h ( x ) O ( h ) Φ + = Φ + ∂ Φ + ∂ Φ + ∂ Φ + x x x 2 6 1 1 2 2 3 3 4 ( x h ) ( x ) h ( x ) h ( x ) h ( x ) O ( h ) Φ − = Φ − ∂ Φ + ∂ Φ − ∂ Φ + x x x 2 6 Forward- and backward difference formula for 1 st -order derivative: 1 2 ( x ) [ ( x h ) ( x )] O ( h ) ∂ Φ = Φ + − Φ + x h 1 st -order accuracy 1 2 ( x ) [ ( x ) ( x h )] O ( h ) ∂ Φ = Φ − Φ − + x h Central-difference formula for 1 st -order derivative: 1 2 nd -order accuracy 3 ( x ) [ ( x h ) ( x h )] O ( h ) ∂ Φ = Φ + − Φ − + x 2 h Central-difference formular for 2 nd -order derivative: 1 3 rd -order accuracy 2 4 ( x ) [ ( x h ) 2 ( x ) ( x h )] O ( h ) ∂ Φ = Φ + − Φ + Φ − + x 2 h

Higher-Order Finite-Difference Schemes From the Taylor expansions: 1 1 1 2 2 3 3 4 4 5 ( x h ) ( x ) h ( x ) h ( x ) h ( x ) h ( x ) O ( h ) Φ + = Φ + ∂ Φ + ∂ Φ + ∂ Φ + ∂ Φ + x x x x 2 6 24 1 1 1 2 2 3 3 4 4 5 ( x h ) ( x ) h ( x ) h ( x ) h ( x ) h ( x ) O ( h ) Φ − = Φ − ∂ Φ + ∂ Φ − ∂ Φ + ∂ Φ − x x x x 2 6 24 9 27 64 2 2 3 3 4 4 5 ( x 3 h ) ( x ) 3 h ( x ) h ( x ) h ( x ) h ( x ) O ( h ) Φ + = Φ + ∂ Φ + ∂ Φ + ∂ Φ + ∂ Φ + x x x x 2 6 24 9 27 64 2 2 3 3 4 4 5 ( x 3 h ) ( x ) 3 h ( x ) h ( x ) h ( x ) h ( x ) O ( h ) Φ − = Φ − ∂ Φ + ∂ Φ − ∂ Φ + ∂ Φ − x x x x 2 6 24 We obtain these relations: 1 3 3 5 ( x h ) ( x h ) 2 h ( x ) h ( x ) O ( h ) Φ + − Φ − = ∂ Φ + ∂ Φ + x x 3 3 3 5 ( x 3 h ) ( x 3 h ) 6 h ( x ) 9 h ( x ) O ( h ) Φ + − Φ − = ∂ Φ + ∂ Φ + x x 5 27 [ ( x h ) ( x h )] [ ( x 3 h ) ( x 3 h )] 48 h ( x ) O ( h ) Φ + − Φ − − Φ + − Φ − = ∂ Φ + x and the fourth-order finite-difference scheme: 9 1 5 ( x ) [ ( x h ) ( x h )] [ ( x 3 h ) ( x 3 h )] O ( h ) ∂ Φ = Φ + − Φ − − Φ + − Φ − + x 16 h 48 h

Finite-Difference Implementation: 1D Case 1D wave equations (2 nd order): 2 k u ρ t u ( ), τ = ∂ ∂ = ∂ τ x x Introducing velocity: v u = ∂ t and we have two 1 st -order equations: ρ t v , k v ∂ = ∂ τ ∂ τ = ∂ x t x 1 Applying the central-difference scheme: 3 ( x ) [ ( x h ) ( x h )] O ( h ) ∂ Φ = Φ + − Φ − + x 2 h n 1 n 1 n n n 1 n 1 n n v v v v + − + − − τ − τ τ − τ − ρ i i i 1 i 1 , i i k i 1 i 1 + − + − = = i i 2 t 2 x 2 t 2 x Δ Δ Δ Δ We get the 2 nd -order finite-difference equations: 1 t Δ n 1 n 1 n n v v ( ) + − = + τ − τ i i i 1 i 1 + − ρ x Δ i t Δ n 1 n 1 n n + − k ( v v ) τ = τ + − i i i i 1 i 1 x + − Δ

Parallelization (128 Processes) 1 t Δ n 1 n 1 n n v + v − ( ) = + τ − τ i i i 1 i 1 ρ x + − Δ i Proc 127 t Δ n 1 n 1 n n k ( v v ) + − τ = τ + − i i i i 1 i 1 + − x Δ Up (Z) � North (Y) � N 4 n = East z z N 4 n N 8 n = (X) � = y y x x 3 n n n x 2 n z y Proc 007 Proc 008 x n x 1 Proc 001 Proc 000

Finite-Difference Method (FDM) { } T  τ C : ε C : [ ( u ) ( u ) ] C : ( u ) u τ f 1 = = ∇ + ∇ = ∇ ρ = ∇ ⋅ + 2 (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.

Spectral-Element Method (SEM)   u τ f ρ = ∇ ⋅ + { T } τ C : ε C : [ ( u ) ( u ) ] C : ( u ) 1 = = ∇ + ∇ = ∇ 2 3 3 3  u d r τ d r f d r ∫ ρ = ∫ ∇ ⋅ + ∫ ⊕ ⊕ ⊕ (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. (1,1,1) (4) Expand all functions (displacement, stress, model, etc.) in spectral basis functions in each cell. N = ∑ i i M cells (elements) with N u ( x ) c B ( x x ), k 1 , 2 ,..., M − = k k k control points in each cell i 1 = (-1,-1,-1) G c 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 Click to show movie Tromp et al. (2011)

Shake Movie: Chi-Chi (Taiwan) Earthquake Near-surface velocity Click to show movie

Thank you!