Earthquake Rupture dynamics in complex geometries using coupled high-order finite difference methods and finite volume methods

Earthquake Rupture dynamics in complex geometries using coupled high-order finite difference methods and finite volume methods

Ossian O´Reilly, Eric M. Dunham, Jeremy E. Kozdon and Jan Nordström

Linköping University Post Print

N.B.: When citing this work, cite the original article.

Original Publication:

Ossian O´Reilly, Eric M. Dunham, Jeremy E. Kozdon and Jan Nordström, Earthquake Rupture dynamics in complex geometries using coupled high-order finite difference methods and finite volume methods, 2012, 2012 AGU Fall Meeting.

Postprint available at: Linköping University Electronic Press

Earthquake Rupture dynamics in complex geometries using coupled high-order finite difference methods and finite volume methods

Ossian J. O’Reilly(1,3), Eric M. Dunham(1), Jeremy E. Kozdon(2) and Jan Nordström(3)

Abstract

We present a 2-D multi-block method for earthquake rupture dynamics in complex geometries using summation-by-parts (SBP) high-order finite differences on structured grids coupled to finite volume methods on unstructured

meshes. The node-centered finite volume method is used on unstructured triangular meshes to resolve earthquake ruptures propagating along non-planar faults with complex geometrical features. The unstructured meshes discretize the fault geometry only in the vicinity of the faults and have collocated nodes on the fault boundaries. Away from

faults, where no complex geometry is present, the seismic waves emanating from the earthquake rupture are resolved using the high-order finite difference method on coarsened structured grids, improving the computational efficiency while maintaining the accuracy of the method.

In order for the SBP high-order finite difference method to communicate with the node-centered finite volume

method in a stable manner, interface conditions are imposed using the simultaneous approximation term (SAT) pen-alty method. In the SAT method the interface conditions and boundary conditions are imposed in a weak manner. Fault interface conditions (rate-and-state friction) are also imposed in a provably stable manner using the SAT

method. Another advantage of the SAT method is the ability to impose multiple boundary conditions at a single boundary node, e.g. at the triple junction of a branching fault.

The accuracy and stability of the numerical implementation are verified using the method of manufactured solutions and against other numerical implementations. To demonstrate the potential of the method, we simulate an

earth-quake rupture propagating in a nonplanar fault geometry resolved with unstructured meshes in the fault zone and structured grids in the far-field.

Introduction

Numerical approach: Hybrid method

Weak boundary conditions and interface conditions

In implementation, characteristic inteface conditions used. Characteristic interface conditions for friction lead to a nonlinear equation to solve for .

(1) Department of Geophysics, Stanford University, CA; (2) Department of Applied Mathematics, Naval Postgraduate School, Monterey, CA; (3) Department of Mathematics, Linköping University, Sweden

Sliding plug in volcanic conduit

Rupture propagation in branching fault geometry

Verification

Computational domain tessellated into multiple blocks. Complex geometry meshed with unstructured meshes. Remaining blocks meshed with structured grids.

Wave speed

Shear impedance

Boundary conditions imposed in a provably stable manner using the SAT method. In the SAT method source terms are added to governing equations.

Dynamic modeling of rupture propagation in complex geometries in earthquake seismology and volcano seismology need to account for non-planar fault geometries. For example:

faults with bends, step-overs, branches, kinks and cylindrical plugs in volcanoes.

Into the page Out of the page Hypocenter Initial conditions Fault Friction (Rate-and-state) Upper section: Lower section: (Velocity weakening) (Velocity weakening) (Velocity strengthening)

Into the page Out of the page

Fault

PP NAP

Kink external side Kink internal side

Solve nonlinear least squares problem on a small disc to estimate A for each problem

Local mesh refinement is used to resolve singularity Crack tip

Stress-state

Crack tip

Kink Kink Kink

Crack tip Crack tip

Cross-section Cross-section −0.01 −0.005 0 0.005 0.01 −0.01 −0.005 0 0.005 0.01

Asymptotic solutions near singularities for kinks and crack tips in elastos-tatics are compared with numerical solutions.

External side Internal side −0.01 −0.005 0 0.005 0.01 −0.01 −0.005 0 0.005 0.01 −0.01 −0.005 0 0.005 0.01 −0.01 −0.005 0 0.005 0.01 FVM Hybrid 4th

Points Log₁₀ Rate Points Log₁₀ Rate 7 4 . 1 -1 2 4 1 7 -1.39 1 2 , 4 16,571 -2.01 2.07 3,094 -2.10 3.11 65,728 -2.64 2.09 4,118 -2.68 2.60 261,776 -3.24 2.02 32,807 -3.33 2.87 1,044,790 -3.83 1.94 105,204 -3.82 1.91 4,174,645 -4.37 1.81 353,169 -4.36 2.08

Particle velocity field (m/s) . Rupture takes the SWFZ branch.

Finite volume method (FVM) in the whole domain vs the hybrid method. Points removed in the grids until the errors and closely match. Analytic tion constructed using method of manufactured solu-tions (MMS). MMS adds extra terms to boundary con-ditions to satisfy a priori known solution.

Contact:

ooreilly@stanford.edu

Log10 radius (m) Log10 radius (m) Log10 error −2.5 −2 −1.5 −1 −0.5 0 −0.5 0 0.5 1 1.5 2

Log10 radius (m) Log10 radius (m) Log10

error

Convergence

Log10 Number of nodes

Log10 Number of nodes

1 m 1 m East Plug Fault interface Hybrid interface

We demonstrate our numerical method motivated by the 2004 M6.0 Parkfield earthquake. Figure 1 shows a cross-section of the San Andreas fault (SAF) look-ing from the south to the north. The Pacific plate (PP) is movlook-ing into the page

and the North-American plate (NAP) is moving out of the page. Observations in-dicate aftershock events on the Southwest fracture zone (SWFZ) branch . In our model, we consider the upper section of the SAF as velocity strengthening, the lower section of the SAF and the SWFZ as velocity weakening.

Plug Fault Hybrid interface Nucleation zone Friction Velocity Cross-section Cross-section Pressure Plug

Normal not defined at corner. Two source terms are imposed; one for each neighboring edge segment. For example: friction at node “a”

−1 −0.5 0 0.5 1 −2 0 2 4 6 8 10 −1 −0.5 0 0.5 1 −2 0 2 4 6 8 10 −1 −0.5 0 0.5 1 −2 0 2 4 6 8 10 −1 −0.5 0 0.5 1 −2 0 2 4 6 8 10 10 km 10 km 100 m 100 m 100 m 100 m 100 m 100 m 100 m 100 m 10 km 10 km 10 km 10 km 10 km 10 km Numerical solution Asymptotic solution Numerical solution Asymptotic solution Numerical solution

Asymptotic solution Numerical solutionAsymptotic solution

−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 1 1.5 2 2.5 3 3.5 4 −4.5 −4 −3.5 −3 −2.5 1 1.5 2 2.5 3 3.5 4 −3.4 −3.2 −3 −2.8 −2.6

Plug sealing of the vent of a volcano . Constant driving force pushes the

plug upward . Sliding resisted by friction. No restoring force accounted for. Horizontal cross-section modeled.

Density

Shear modulus Antiplane shear

Particle velocity field. Color scale is saturated to emphasize wave field in the body. Rupture propagates downward along the smooth margin of the plug.

Slip velocity as a function of angular coordinate. Wave focusing after ~60 ms amplifies particle velocity in plug. Friction parameters

(m/s)

t = 42 ms t = 62 ms t = 76 ms t = 112 ms

One application problem for our numerical method is motivated by repeated drumbeat earthquakes as observed during the Mt. St. Helens 2004 - 2008

eruption. 40 0 (MPa) 10 0 (MPa)

Numerical solution, kink,

Numerical solution, crack tip

1 m 1 m 1 m 1 km a b (m/s) -5 0 5 (m/s) -1 0 1 (m/s) -1 0 1 (m/s) -1 0 1 (m/s) -1 0 1 (m/s) -1 0 1 (m/s) -1 0 1 (m/s) -1 0 1 (m/s) -1 0 1 (m/s) -5 0 5 (m/s) -5 0 5 (m/s) -5 0 5 (m/s) -5 0 5 (m/s) -5 0 5 (m/s) -5 0 5 (m/s) -5 0 5 t = 0.5 s t = 6 ms t = 16 ms t = 28 ms t = 40 ms t = 50 ms t = 62 ms t = 78 ms t =90 ms t = 1.5 s t = 2.0 s t = 2.5 s t = 3.0 s t = 3.5 s t = 4.5 s t = 6.0 s