Off-fault plasticity in three-dimensional dynamic rupture simulations using a modal Discontinuous Galerkin method on unstructured meshes: implementation, verification and application

SUMMARY

The dynamics and potential size of earthquakes depend crucially on rupture transfers between adjacent fault segments. To accurately describe earthquake source dynamics, numerical models can account for realistic fault geometries and rheologies such as nonlinear inelastic processes off the slip interface. We present implementation, verification and application of off-fault Drucker–Prager plasticity in the open source software SeisSol (www.seissol.org). SeisSol is based on an arbitrary high-order derivative modal Discontinuous Galerkin method using unstructured, tetrahedral meshes specifically suited for complex geometries. Two implementation approaches are detailed, modelling plastic failure either employing subelemental quadrature points or switching to nodal basis coefficients. At fine fault discretizations, the nodal basis approach is up to six times more efficient in terms of computational costs while yielding comparable accuracy. Both methods are verified in community benchmark problems and by 3-D numerical h- and p-refinement studies with heterogeneous initial stresses. We observe no spectral convergence for on-fault quantities with respect to a given reference solution, but rather discuss a limitation to low-order convergence for heterogeneous 3-D dynamic rupture problems. For simulations including plasticity, a high fault resolution may be less crucial than commonly assumed, due to the regularization of peak slip rate and an increase of the minimum cohesive zone width. In large-scale dynamic rupture simulations based on the 1992 Landers earthquake, we observe high rupture complexity including reverse slip, direct branching and dynamic triggering. The spatiotemporal distribution of rupture transfers are altered distinctively by plastic energy absorption, correlated with locations of geometrical fault complexity. Computational cost increases by 7 per cent when accounting for off-fault plasticity in the demonstrating application. Our results imply that the combination of fully 3-D dynamic modelling, complex fault geometries and off-fault plastic yielding is important to realistically capture dynamic rupture transfers in natural fault systems.

Numerical approximations and analysis, Numerical modelling, Earthquake dynamics, Earthquake hazards, Fractures, faults, and high strain deformation zones

1 INTRODUCTION

Understanding the physics of earthquake source processes enhances seismic hazard assessment for natural fault systems. Specifically, directivity effects as well as potential rupture transfers to adjacent fault segments are determined by earthquake source dynamics. However, the fundamentals of rupture dynamics are difficult to infer from observations. Numerical simulations pose a powerful tool to further our understanding of earthquake rupturing (complex) faults.

In modelling earthquake rupture dynamics, the fault geometry in conjunction with fault stress and strength constitute essential initial conditions, determining frictional failure, rupture propagation and seismic wave emanation off the fault. Additionally, such models can describe the interaction of fault slip with the surrounding host-rock material, for example, by considering off-fault plastic deformation.

Accounting for plastic rock failure causes inelastic energy dissipation which, in turn, influences rupture dynamics (e.g. Andrews 1976, 2005; Kaneko & Fialko 2011). Moreover, unreasonably high slip velocities on the fault are limited (Andrews 2005; Dunham et al.2011a). Plasticity also affects the source mechanical characteristics such as rupture speed and rupture style (e.g. Duan & Day 2008; Templeton & Rice 2008; Dunham et al.2011a; Gabriel et al.2013). For example, the transition from pulse-like to crack-like, as well as from sub- to supershear rupture may be delayed or even prevented by plastic material response.

Dynamic rupture simulations reveal that off-fault plastic yielding is enhanced by the interaction of rupture with a free surface boundary condition (Ma 2008; Ma & Andrews 2010). As a consequence, the formation of depth-dependent flower-like distribution of plastic strain is observed (Ma & Andrews 2010), which is consistent with observations of damage zones in natural fault systems (e.g. Chester et al.1993; Ben-Zion & Sammis 2003; Mitchell & Faulkner 2009). In addition, off-fault plasticity restricts shallow fault slip (Kaneko & Fialko 2011; Erickson et al.2017; Roten et al.2017), contributing to the well-observed shallow slip deficit for strike-slip events (e.g. Fialko et al.2005). Purely kinematic earthquake scenarios of the Southern San Andreas Fault (Roten et al.2014), as well as spontaneous dynamic rupture on simplified planar faults (Roten et al.2015) confirm that plasticity impacts peak slip rates, dominantly close to the free surface. As a result, peak ground velocities (PGVs) are reduced dramatically compared to purely elastic modelling.

Previous studies including plastic deformation are mostly based on planar fault geometries. However, geological observations reveal that natural faults are complex geometrical systems, which may include bends, branches and distinct fault segments (e.g. King & Nabelek 1985; Wesnousky 1988, 2006). Accounting for full geometric fault complexity jointly with a homogeneous (regional tectonic) background stress state in the modelling domain leads to highly heterogeneous initial fault stresses, which influences rupture propagation (e.g. Aochi & Fukuyama 2002; Gabriel & Pelties 2014) and enhances plastic yielding (e.g. Dunham et al.2011b).

Combining inelastic processes and complex fault geometries may crucially affect rupture transfers in terms of branching and dynamic triggering (‘jumping’). For example, plasticity inhibits the activation of branching segments on compressional sides, but promotes branching on extensional sides of strike-slip fault events (DeDontney et al.2012). Considering 2-D discontinuous, yet planar fault segments (step-over faults), off-fault plastic failure has been shown to effectively lengthen faults by enhancing coseismic slip and increasing the slip gradient at fault tips (Nevitt & Pollard 2017).

In spite of off-fault plasticity being widely studied, its numerical implementation remains challenging. Plastic material behaviour introduces a nonlinear constitutive equation relating stress and strain. Incorporating this relation directly into the underlying wave equation requires distinct (nonlinear) numerical methods and additional stabilization techniques for solving this problem based on a damage rheology description (e.g. Lyakhovsky & Ben-Zion 2014). Alternatively, plasticity can be readily implemented in already existing numerical solvers of the linear wave equation via a predictor–corrector approach: First, an elastic trial stress state needs to be calculated and checked against the corresponding elastic yield surface. In a second step, the trial stress state is adjusted whenever this state exceeds the elastic domain, that is, plastic deformation occurs. As a consequence, additional calculations are required for every time step and for the complete discretization space independent of plastic yielding actually occurring. Thus, incorporating plasticity can be computationally demanding. In particular, additional interpolation may become necessary to define stresses on all spatial discretization points, for example, in staggered-grid methods. For instance, Roten et al. (2016) report an increase of computational costs by 65 per cent in comparison to a purely elastic simulation.

In addition, numerical solutions of plasticity models tend to be mesh dependent regarding strain localization (e.g. Templeton et al.2009; Dunham et al.2011a), and thus require regularization. To ensure reliable simulation results with mesh refinement, numerical solvers need to be verified by convergence tests. However, detailed convergence studies of 3-D dynamic rupture problems are generally scarce, and are specifically missing with respect to heterogeneous on-fault initial conditions and off-fault plasticity.

Simulating dynamic rupture requires a discretized model of the area of interest, including a prescribed fault surface. Numerical methods based on hexahedral element discretization of the domain such as Finite Differences (FD; e.g. Day 1982; Dalguer & Day 2007; Cui et al.2010), Spectral Element methods (Kaneko et al.2008; Galvez et al.2014) or Spectral Boundary Integral Equation (BIE) methods (e.g. Lapusta et al.2000) are often restricted to planar fault models. The incorporation of curvilinear elements (Cruz-Atienza & Virieux 2004; Kozdon et al.2012; Zhang et al.2014; Duru & Dunham 2016) enables the discretization of geometrical heterogeneities such as fault roughness across fault segments, whereas fault branching remains challenging to incorporate. Finite Element (FE) methods based on tetrahedral elements (e.g. Duan & Oglesby 2006; Ma 2008; Barall 2009; Pelties et al.2012) as well as BIE methods (e.g. Rice 1993; Aochi et al.2000; Ando 2016) are well suited for complex fault geometries. Discontinuous Galerkin (DG) and BIE methods are the only methods able to accurately handle intersecting faults (e.g. Pelties et al.2012; Tago et al.2012; Ando et al.2017). However, the BIE method is restricted to purely elastic material properties. The DG method is, therefore, increasingly becoming attractive as a method of choice for representing Earth’s complex structure in a high-order accurate manner (e.g. Wilcox et al.2010; Pelties et al.2012; Tago et al.2012; Duru et al.2017; Tavelli & Dumbser 2018).

Here, we analyse numerical and physical characteristics of dynamic rupture simulations accounting for off-fault plastic yielding in direct conjunction with complex 3-D fault geometry as well as heterogeneous initial stress and strength conditions. Such simulations pose high demands on numerical methods in terms of computational costs and mesh generation. 3-D models featuring highly resolved fault zones, necessary to accurately capture the dynamic rupture process, exhibit quickly millions of degrees of freedom.

To this end, we present the implementation of plastic yielding in the open-source software package SeisSol via a Return Mapping Algorithm (e.g. Ortiz & Simo 1986). SeisSol is based on an Arbitrary high-order DERivative DG (ADER-DG) method (e.g. Dumbser & Käser 2006; Pelties et al.2012). The software package is highly optimized for the efficient use on modern high-performance computing infrastructure (Breuer et al.2014; Heinecke et al.2014; Breuer et al.2015, 2016; Heinecke et al.2016; Uphoff & Bader 2016) enabling the simulation of realistic large-scale dynamic rupture models.

The remainder of this paper is structured as follows: In Section 2 we detail two distinct implementations of plastic material response based on a non-associated Drucker–Prager plastic yield criterion. Both implementations are suitable for any high-order modal DG approach and not restricted to the ADER-DG method. We verify both implementations in Section 3 in well-established community benchmarks for two naturally arising faulting mechanisms (Harris et al.2009, 2011). In Section 4, we investigate the numerical convergence of key dynamic rupture parameters with mesh refinement and increasing polynomial order with respect to a reference solution. The setup includes depth-dependent initial stresses and plastic yielding, to reflect the current state-of-the-art of realistic dynamic earthquake simulations. The 3-D h- and p-refinement study extends previous studies considering only elastic material response and homogeneous initial stress conditions (e.g. Day et al.2005; Pelties et al.2012). Both plasticity implementations are then analysed in terms of efficiency (computational cost) versus accuracy.

We demonstrate the specific advantages of the modal DG approach on unstructured meshes by presenting a large-scale 3-D dynamic rupture earthquake scenario based on the 1992 Landers earthquake. We focus on the influence of inelastic material behaviour on rupture transfers across a geometrical complex fault system. Furthermore, we compare macroscopic earthquake source characteristics, such as peak slip rate and the distribution of slip in purely elastic versus plastically yielding simulations. Plasticity highly impacts 3-D source dynamics, specifically at locations of geometrical complexities such as fault branches and fault bends. We critically discuss our numerical analysis as well as application example and conclude that it may be essential to combine fully 3-D dynamic modelling, complex fault geometries and off-fault plastic yielding to realistically capture dynamic rupture transfers in natural fault systems.

2 METHODOLOGY

In the first part of this section (Section 2.1), we summarize the governing equations, underlying numerical method and the recent optimization of the software package SeisSol (www.seissol.org) simulating spontaneous earthquake rupture coupled to elastic seismic wave propagation. We then detail the theory of non-associated Drucker–Prager plastic yielding and its numerical regularization (Section 2.2). In Section 2.3 we present two approaches for the numerical implementation of plastic yielding in modal DG methods, and discuss their computational efficiency.

2.1 Spontaneous earthquake rupture within the ADER-DG framework

2.1.1 Numerical method solving the elastic wave equation

SeisSol numerically solves the weak form of the elastic wave equation in velocity–stress formulation. The underlying system of equations can be written in a compact matrix–vector form as

\begin{equation*} \frac{\partial }{\partial t} \boldsymbol {Q} + \boldsymbol {A} \frac{\partial }{\partial x} \boldsymbol {Q} + \boldsymbol {B} \frac{\partial }{\partial y} \boldsymbol {Q} + \boldsymbol {C} \frac{\partial }{\partial z} \boldsymbol {Q} = \boldsymbol {0}, \end{equation*}

(1)

for the solution |$\boldsymbol {Q}=(\sigma _{xx}, \sigma _{yy},\sigma _{zz}, \sigma _{xy}, \sigma _{yz},\sigma _{xz}, u, v,w)^T$| including the stress tensor components σ_ij and the velocities u, v, w in the x, y, z directions, respectively. The space-dependent Jacobian matrices |$\mathbf {A},\mathbf {B},\mathbf {C}$| contain the Lamé parameters λ and μ as well as the density ρ, encapsulating the elastic material properties. A detailed definition can be found in Dumbser & Käser (2006).

For spatial discretization, a modal DG approach is employed. The computational domain is subdivided into first-order tetrahedral elements τ_m and all material properties are constant within an element. To approximate the solution |$\boldsymbol {Q}$| at a point |$\boldsymbol {x}= (x,y,z)^T$| and at time t we use a linear combination of orthogonal polynomial basis functions Φ_i, namely the Dubiner’s basis functions (Cockburn et al.2000), and time-dependent coefficients |$\boldsymbol {\hat{Q}}_i(t)$|⁠:

\begin{equation*} \boldsymbol {Q}(\boldsymbol {x},t) = \sum _{i=1}^{L} \boldsymbol {\hat{Q}}_i(t) \ \Phi _i(\boldsymbol {x}), \end{equation*}

(2)

with L = (p + 1)(p + 2)(p + 3)/6 for polynomial degree p or order of accuracy |$\mathcal {O}=p+1$|⁠. Multiplication of eq. (1) by a test function Φ_k and integration over one element τ_m leads to the semi-discrete weak formulation. To evaluate the resulting mass, stiffness and flux matrices (see Dumbser & Käser 2006) independently of the element shape of τ_m, all elements are transformed to a tetrahedral reference element τ_E. These matrices can then be pre-calculated analytically leading to a quadrature-free scheme (e.g. Atkins & Shu 1996).

The exchange of information between elements is purely local, based on the concept of numerical fluxes established in Finite Volume methods: SeisSol uses the exact solution of the Riemann problem, namely the Godunov flux. This upwind flux poses a differential equation with discontinuous initial conditions (Toro 1999; LeVeque 2002). The numerical properties of the ADER-DG algorithm are extremely sensitive to the choice of flux function. The numerical dissipation and dispersion properties of the Godunov upwind flux have been widely studied (e.g. Hu et al.1999; Käser et al.2008; Hesthaven & Warburton 2010). While numerical dissipation leads to a decay in waveform amplitudes of the marginally resolved wave numbers, numerical dispersion leads to a gradual separation of waves of different wave periods. Our ADER-DG scheme offers favourably low numerical dispersion properties, allowing to accurately recover phase velocities propagating over a large number of wavelengths (Käser et al.2008). Numerical dissipation decreases with increasing polynomial degree and is stronger beyond a cut-off frequency that depends on the mesh size h and on the order of accuracy |$\mathcal {O}$|⁠. The cut-off frequency is expected to be inversely proportional to the traveltime of s-waves over a typical grid spacing |$\approx h/\mathcal {O}/V_S$| (Pelties et al.2014).

For the integration in time, the DG method is combined with an ADER scheme (Titarev & Toro 2002; Dumbser & Käser 2006; Käser & Dumbser 2006), which provides equivalent high-order accuracy as in space using a single explicit time integration step.

2.1.2 Dynamic rupture as internal boundary condition

De la Puente et al. (2009) and Pelties et al. (2012) incorporate nonlinear frictional failure as internal boundary condition into SeisSol in 2-D and 3-D. The Coulomb failure criterion and consecutive constitutive laws are enforced by solving a modified inverse Riemann problem on prescribed element interfaces (fault surfaces) in contrast to the typically applied traction at split-node approach (e.g. Andrews 1999). This requires, in distinction to the quadrature-free approach used for solving the wave equation, the introduction of space-time quadrature points (QP) across dynamic rupture interfaces. The flux functions across those interfaces are integrated by quadrature based on (p + 2)² Gaussian points (Stroud 1971). These points are located within the triangular faces of the tetrahedral elements connected to the fault surface. Their distribution is visualized for p = 3 in Fig. C1 in Appendix  C. At these element internal points, we also enforce the Coulomb failure criterion of the frictional boundary condition. Numerical convergence towards a reference solution was demonstrated for subelemental treatment of the dynamic rupture boundary condition in Pelties et al. (2012). Fault initial stress and friction parametrizations are assigned individually to each QP, resulting in an excellent agreement with established numerical methods even under highly heterogeneous conditions (Pelties et al.2014).

SeisSol is verified for a wide range of advanced dynamic rupture problems such as branched faults, dipping fault geometries and laboratory-derived constitutive laws such as the rate-and-state friction law (Pelties et al.2014). It allows for unstructured, tetrahedral mesh discretization, facilitating automatized mesh generation for naturally arising fault geometries, topography and subsurface structures. Additionally, local mesh refinement and coarsening allow for adaptive resolution, for example, to resolve the stress drop at the rupture tip in the cohesive zone (process zone; Ida 1972).

SeisSol’s on-fault slip rates remain notably free of high-frequency oscillations (De la Puente et al.2009; Pelties et al.2012). The generation of such non-physical high-frequency modes may contaminate the solution over all space-time scales (Duan & Day 2008). Due to the numerical properties of our flux, higher frequency modes are subdued while the physically meaningful lower frequencies are minimally affected (see Section 2.1.1). This is advantageous for dynamic rupture simulations: spurious high-frequency oscillations do not affect the vicinity of the fault and typical damping procedures as used by other methods (e.g. Rojas et al.2008) do not need to be applied. The presence of spurious oscillations in DG schemes using a different flux, such as the non-dissipative central flux (Tago et al.2012), suggests that on-fault solution behaviour is directly related to the properties of the numerical flux. Still, it remains unclear whether this behaviour results from the dissipative properties of the (modified) flux across the fault boundary elements (on-fault) or the upwind flux across non-frictional elements (off-fault) or both.

2.1.3 Computational optimization

SeisSol is optimized for the latest CPU architectures. It delivers high performance on a single compute node as well as on thousands of compute nodes in parallel. For example, SeisSol was one of the first software packages that were optimized for Intel’s Knights Landing architecture (Heinecke et al.2016). It exceeded 1 PFLOP s⁻¹ performance on the SuperMUC supercomputer (Breuer et al.2014), and it scaled up to 1.6 million compute cores of the Tianhe-2 supercomputer (Heinecke et al.2014).

Recent optimization include a hybrid OpenMP/MPI (Open Multi-Processing/Message Passing Interface) parallelization (Heinecke et al.2014), high-performance compute kernels (Breuer et al.2014), asynchronous I/O (Rettenberger & Bader 2015; Rettenberger et al.2016) and clustered local time-stepping (Breuer et al.2016).

SeisSol is suitable for large-scale (spatial extend and modelling duration) earthquake simulations (Pelties et al.2012; Gabriel & Pelties 2014; Heinecke et al.2014; Rettenberger et al.2016; Weingärtner et al.2016; Madden et al.2017) including modelling challenges due to the geometrical complexity of the Earth, for example, shallowly dipping megathrust faults, topography, 3-D subsurface structure and fault roughness (Ulrich & Gabriel 2017; van Zelst et al.2017).

All previous SeisSol dynamic rupture simulations were based on the assumption of purely elastic material properties. Motivated by the recent gain in efficiency, we present in this work the incorporation of the physics of plastic deformation into SeisSol.

2.2 Off-fault plastic yielding

This section provides an overview of the physics and the derivation of a numerical update scheme for off-fault plasticity. We first summarize a non-associated plasticity model widely used to describe plastic material behaviour of soils and rocks. We then discuss the required numerical regularization using viscoplasticity with a mesh-independent regularization factor. Finally, we describe the implementation of viscoplasticity via a Return Mapping algorithm into an existing dynamic rupture solver.

2.2.1 Non-associated Drucker–Prager plasticity

Plastic material response can be represented by a yield function F, which defines the onset of plastic yielding and thus limits the elastic domain of a material, and by a plastic potential function g defining the direction of the plastic strain rate (plastic flow) in case of plastic yielding. The corresponding plastic strain rate can be defined as the derivative of the flow rule with respect to the stress. In so-called associated plasticity models, for example, used for metals and alloys (e.g. Vermeer & de Borst 1984), the plastic strain increment is collinear to the normal of the yield surface, hence g = F. However, to model the plastic behaviour of soil and rocks, non-associated plasticity formulations are used, in which F and g are defined individually.

The total strain |$\boldsymbol {\epsilon }$| is the sum of an elastic strain component and a plastic strain component, |$\boldsymbol {\epsilon } = \boldsymbol {\epsilon }^e + \boldsymbol {\epsilon }^p$|⁠, whereas |$\boldsymbol {\epsilon } = \boldsymbol {\epsilon }^e$| in case of pure elasticity. The stress increment can then be expressed by the time derivative of Hooke’s law:

\begin{equation*} \dot{\sigma }_{ij} = \sum _{k,l}C_{ijkl} (\dot{\epsilon }_{kl} - \dot{\epsilon }^p_{kl}), \end{equation*}

(3)

for the isotropic, elastic, fourth-order stiffness tensor C. In the following, we derive the plastic strain increment

\begin{equation*} \sum _{k,l}C_{ijkl} \dot{\epsilon }^p_{kl} = \lambda \delta _{ij}\left(\sum _k\dot{\epsilon }_{kk}^p\right) + 2 \mu \dot{\epsilon }_{ij}^p, \end{equation*}

(4)

where δ_ij is the Kronecker delta. We consider a plasticity model without dilatancy, that is, no volumetric changes due to plastic yielding, such that |$\sum _i \dot{\epsilon }_{ii}^p=0$| holds true.

The formulation employs a Drucker–Prager yield criterion which is defined as

\begin{equation*} \tau _{c}= c\ \cos (\phi ) - \sigma _m \ \sin (\phi ), \end{equation*}

(5)

for cohesion c, an internal angle of friction ϕ = tan ⁻¹(⁠|$v$|⁠) with bulk friction |$v$| and mean stress |$\sigma _m = (\sum _{i=1}^3\sigma _{ii})/3$|⁠, assuming compressional stresses to be negative. Defining the second invariant of deviatoric stresses s_ij as

\begin{equation*} I_2= \frac{1}{2}\sum _{i,j} s_{ij} s_{ji} = \frac{1}{2}\sum _{i,j} (\sigma _{ij}- \delta _{ij}\ {\sigma _m})(\sigma _{ji}- \delta _{ji}\ {\sigma _m}), \end{equation*}

(6)

the yield function can be expressed as

\begin{equation*} F(\boldsymbol {\sigma }) = \sqrt{I_2} - \tau _c. \end{equation*}

(7)

If the current stress state |$\boldsymbol {\sigma }$| reaches the elastic limit (⁠|$F(\boldsymbol {\sigma })= 0$|⁠), the plastic strain rate |$\boldsymbol {\dot{\epsilon }}^p$| is determined with the plastic potential function |$g(\boldsymbol {\sigma }) = \sqrt{I_2}$| as follows:

\begin{equation*} \dot{\epsilon }^p_{ij} = \frac{1}{2\mu } \frac{\partial g(\mathbf {\sigma })}{\partial \sigma _{ij}} = \frac{1}{2\mu } \frac{s_{ij}}{2\sqrt{I_2}}. \end{equation*}

(8)

With this, we define a so-called perfect plasticity model: The condition F = 0 is enforced strongly by instant multiplication of the deviatoric stresses with |$\tau _c / \sqrt{I_2}$| when the yield surface is reached.

In case of plastic yielding, plastic strain at time t can be mapped into the scalar quantity η(t) following Ma (2008):

\begin{equation*} \eta (t) = \int _0^t \mathrm{ d}\eta = \int _0^t \sqrt{\frac{1}{2} \dot{\epsilon }_{ij}^p \dot{\epsilon }^p_{ij}}. \end{equation*}

(9)

Fig. 1 depicts the exemplary accumulated plastic strain η surrounding a strike-slip fault in SCEC benchmark TPV27 (defined in Section 3).

Figure 1.

SeisSol’s simulated accumulated plastic strain, as defined in eq. (9), surrounding a strike-slip fault (TPV27, after 8 s of simulated time, see Section 3). At depth, plastic strain occurs mainly on the compressional side of the fault, but shifts to both fault sides close to the surface. SeisSol uses an unstructured tetrahedral mesh featuring coarsening away from the fault.

2.2.2 Viscoplastic regularization

Perfect plasticity models, as discussed in the last section, are known to be mathematically ill-posed; their numerical solution tends to show mesh-dependent behaviour in case of strain localization (e.g. De Borst et al.1996; Templeton & Rice 2008; Dunham et al.2011a; Xu et al.2012). Dias da Silva (2004) summarizes different approaches enabling FE methods to produce mesh-independent results or, at least, to control mesh dependencies. These include employing Cosserat media, gradient plasticity theories or the introduction of a rate-dependent material behaviour emulated by viscoplasticity. In the latter case, stresses are allowed to exceed the yield criterion and are subsequently relaxed to the yield surface over a specified amount of time (relaxation time).

In dynamic rupture simulations, viscoplastic relaxation is widely used to regularize the numerical implementation of off-fault plastic yielding (Andrews 2005; Duan 2008; Templeton & Rice 2008; Dunham et al.2011a; Xu et al.2012; Gabriel et al.2013). We employ hereafter a rate-dependent Duvaut–Lions (Duvaut & Lions 1976) formulation of viscoplasticity, as specified for SCEC benchmark test problems (Harris et al.2009; Barall & Harris 2014). The viscoplastic strain rate is given as the differential equation

\begin{equation*} \dot{\epsilon }^{vp}_{ij} = \frac{1}{2\mu T_v}(\sigma _{ij}-P_{ij}(\boldsymbol {\sigma })), \end{equation*}

(10)

for the viscoplastic relaxation time T_|$v$| > 0 and projection P with

\begin{equation*} P_{ij}(\boldsymbol {\sigma }) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(\tau _c / \sqrt{I_2}) \ s_{ij} + \sigma _m \ \delta _{ij} &\text{for } F(\boldsymbol {\sigma }) \ge 0\\ \sigma _{ij} &\text{for } F(\boldsymbol {\sigma }) < 0. \end{array}\right. \end{equation*}

(11)

Note that |$P_{ij}(\boldsymbol {\sigma })$| is the adjusted stress state in the rate-independent plasticity case with |$\dot{\epsilon }^{p}_{ij}$| defined in eq. (8). In the absence of additional plastic yielding, viscoplastic regularization permits the current stress state σ_ij to reach the inviscid stress state |$P_{ij}(\boldsymbol {\sigma })$| after a specified time T_|$v$|. In distinction to previous studies, we choose a constant relaxation time T_|$v$| independent of spatial discretization as detailed in the next section.

2.2.3 Implementation via a Return Mapping algorithm

Incorporating plasticity is commonly achieved via Return Mapping algorithms (e.g. Andrews 2005, for dynamic rupture solvers), requiring the following steps:

First, a trial stress state |$\boldsymbol {\sigma }^{\text{trial}}$| at time step t_{n + 1} is calculated, which is based on the numerical solution of the wave equation coupled to frictional boundary conditions assuming purely elastic material properties (i.e. |$\boldsymbol {\dot{\epsilon }}^p=0$| in eq. 3). Second, if |$F(\boldsymbol {\sigma }^{\text{trial}}) < 0$|⁠, no plastic yielding occurs and the trial stress becomes the new stress, |$\boldsymbol {\sigma }^{n+1}=\boldsymbol {\sigma }^{\text{trial}}$|⁠. In case of plastic yielding (⁠|$F(\boldsymbol {\sigma }^{\text{trial}}) \ge 0$|⁠) the stress state |$\boldsymbol {\sigma }^\text{trial}$| is updated in a plasticity corrector-step by assuming non-zero plastic strain rates (i.e. |$\boldsymbol {\dot{\epsilon }}^p \ne 0$|⁠)

\begin{equation*} {\sigma }_{ij}^{n+1} = \sigma _{ij}^{\text{trial}} - 2 \mu ({\epsilon }_{ij}^{vp})^{n+1}. \end{equation*}

(12)

To determine |${\epsilon }_{ij}^{vp}$| at t_{n + 1} we need to integrate the viscoplastic strain rate in eq. (10) over one time step. In classic Return Mapping algorithms an implicit backward Euler scheme is used for numerical integration (Simo et al.1988). This approach has been presented for a viscoplastic rheology in dynamic rupture simulations by Dunham et al. (2011a). For the plasticity formulation presented here, an updated scheme based on a closed form integration of eq. (12) is used (Andrews 2005; Duan & Day 2008). Therein, the problem is solved explicitly by assuming the projection P_ij to be constant over one time step.

The full derivation of the updated stress state |$\sigma _{ij}^{n+1}$| is detailed in Appendix  A and summarized here as

\begin{equation*} \sigma _{ij}^{n+1} = f^* s_{ij}^{\text{trial}} + \sigma _{m}^{\text{trial}}\delta _{ij} , \end{equation*}

(13)

with the adjustment factor

\begin{equation*} f^* = (1-\exp (-\Delta t/T_v))\frac{\tau _c}{\sqrt{I_2}} + \exp (-\Delta t/T_v), \end{equation*}

(14)

for time step width Δt. Combining eqs (13) and (12), the viscoplastic strain is then defined as

\begin{equation*} (\epsilon _{ij}^{vp})^{n+1} = \frac{1}{2\mu }(1-f^*) \ s_{ij}^{\text{trial}}. \end{equation*}

(15)

Due to the assumption of no change of volumetric plastic strain, the same mean stress before and after the plastic adjustment is preserved.

In previous studies, T_|$v$| in eq. (14) is chosen to be of the order of the P- or S-wave traveltime across the discretization length dx of the modelling domain (e.g. Duan 2008; Ma & Andrews 2010; Xu et al.2012; Gabriel et al.2013). This approach results in a discretization-dependent T_|$v$|, that is, in a shorter relaxation time for meshes of smaller discretization lengths. However, the time step Δt, discretizing an explicit set of hyperbolic partial differential equations in time, needs to satisfy the Courant–Friedrich–Lévy criterion to guarantee a stable numerical solution. In that case the time step is determined according to the fastest wave velocity in the medium and the smallest discretization length dx. Therefore, Δt itself introduces a mesh-dependent parameter into eq. (14). Thus, throughout this study, we use a constant relaxation time T_|$v$| = 0.03 independent of the spatial discretization. This value has been chosen to regularize plasticity for typical dynamic rupture setups (Harris et al.2011). Note that this choice is not connected to a physically motivated description of viscoplasticity but is a purely numerical factor. We prove convergence towards a reference solution of this parametrization choice with increasing polynomial degree and decreasing mesh size in Section 4.

2.3 Novel approaches for incorporating plasticity in modal DG schemes

In most numerical methods used for solving dynamic rupture problems, such as FD (e.g. Andrews 2005; Kozdon et al.2012), FE (e.g. Duan 2008) or nodal high-order FE (e.g. Kaneko et al.2008; Tago et al.2012), the unknowns of the system of equations define the solution at mesh points (FD, FE) or at polynomial interpolation points within elements (nodal high-order FE). Such points can be readily used to evaluate a plastic yield criterion.

In SeisSol’s modal approach, the unknowns are given by the coefficients of a local polynomial expansion which are not associated with spatial points within an element. Hence, it is impossible to directly apply a nonlinear function at mesh or interpolation points.

We thus present in the following two implementation approaches for incorporating plastic yielding:

1. evaluation of the stress state at subelemental quadrature points and

2. switching to a system of nodal basis functions.

Both methods are fundamentally similar. The yield criterion is first evaluated at a set of subelemental points, followed by a projection of a non-polynomial function (the adjusted stresses) back to the polynomial space of the underlying high-order numerical discretization. To this end, our first approach approximates an L²-projection using a quadrature rule, whereas the second approach uses interpolation. Both schemes can be generalized to any DG approach and are not restricted to the modal ADER-DG method of SeisSol. Even nodal methods need to decide whether the nonlinearity is enforced at the present interpolation points or if a different set of points might be used.

We analyse both approaches in terms of efficiency (computational cost) in Section 2.3.3 and in terms of accuracy in Section 4.4.

2.3.1 Quadrature points approach

For both approaches, we need to retrieve a trial stress tensor |$\boldsymbol {\sigma }^{\text{trial}}$| to be evaluated against the yield function F. We here detail the approximation of the required stress state at quadrature points |$\boldsymbol {\xi }_1, \ldots ,\boldsymbol {\xi }_n$|⁠. This approach is referred to in the following as the QP approach.

The QPs are based on a multidimensional quadrature formulae of degree (2n − 1) for n³ points for tetrahedral elements (Stroud 1971) using Gauss–Jacobi polynomials of degree p. The number of QPs increases with increasing p as n_QP = (p + 2)³ corresponding to a quadrature of degree 2p + 3. Their distribution for p = 2 across a tetrahedral reference element is exemplarily visualized in Fig. 2. We achieve resolution of the solution distinctively higher than given by the vertices of the computational mesh. The trial stress is evaluated at every QP via the following matrix–matrix product:

\begin{equation*} \sigma _{iq}^{\text{QP}} = \sum _{l=1}^L \Phi _l(\boldsymbol {\xi }_i) \hat{Q}_{lq} =: \sum _{l=1}^L G_{il} \hat{Q}_{lq}, \end{equation*}

(16)

where |$\hat{Q}_{lq} := (\boldsymbol {\hat{Q}}_l)_q$| is the coefficient of the lth basis function in the polynomial expansion of the qth quantity and |$\mathbf {G}$| is an n_QP × L matrix. We obtain the n_QP × 6 matrix |$\boldsymbol {\sigma }^{\text{QP}}$|⁠, which contains the stress tensor evaluated at every QP.

Figure 2.

Reference tetrahedron with 64 QPs of the QP approach (blue triangles) and 10 NB points of the NB approach (red circles) for polynomial degree p = 2. Both approaches achieve subelement resolution. QPs cluster towards element edges.

After checking the plastic yield criterion for each QP separately, we adjust the stresses in case of plastic yielding using the update formula defined in eq. (13). To recover the polynomial coefficients from the non-polynomial adjusted stresses |$\boldsymbol {\sigma }^{\text{new}}$| we can make use of an L² projection for the qth quantity requiring

\begin{equation*} \int _{\tau _E} \sigma _q^{\text{new}} \ \Phi _l \ \mathrm{ d}V = \int _{\tau _E} \sum _{k=1}^{L}(\hat{Q}_{kq}^{\text{new}} \Phi _k) \Phi _l \ \mathrm{ d}V \end{equation*}

(17)

for all basis functions Φ_l (i.e. the best polynomial representation of |${\sigma _q^{\text{new}}}$| in the L² sense). Due to the orthogonality of the basis functions, the integral on the right side of the equation is only non-zero for k = l. We can recover the coefficients |$\hat{Q}_{lq}^{\text{new}}$| of the polynomial representation by solving eq. (17)

\begin{equation*} \hat{Q}_{lq}^{\text{new}} = \frac{\int _{\tau _{E}}\sigma _{q}^{\text{new}} \ \Phi _l \,\mathrm{d}V}{\int _{\tau _{E}} \Phi _l \Phi _l \,\mathrm{d}V} \approx \frac{\sum _{i=1}^{n_{\text{QP}}} \sigma _{iq}^{\text{new}} \ w_i \Phi _l(\boldsymbol {\xi }_i) \,\mathrm{d}V}{\int _{\tau _{E}} \Phi _l \Phi _l \,\mathrm{d}V}, \end{equation*}

(18)

where the integral in the numerator is approximated by quadrature using the same QPs |$\boldsymbol {\xi _i}$| in conjunction with corresponding weights |$w$|_i. Hence, the choice of approximating the solution at the specific points |$\boldsymbol {\xi }$| is motivated by the (numerical) evaluation of the integral in eq. (18). The denominator is simply the (analytically precomputed) lth entry M_ll of the diagonal mass matrix |$\mathbf {M}$| as defined in Dumbser & Käser (2006). The L² projection can also be written as matrix–matrix product,

\begin{equation*} \hat{Q}_{lq}^{\text{new}} = \sum _{i=1}^{n_{\text{QP}}} \tilde{G}_{li} \sigma _{iq}^{\text{new}}, \end{equation*}

(19)

with the L × n_QP matrix |$\mathbf {\tilde{G}}$| with entries |$\tilde{G}_{li} = \frac{1}{M_{ll}}w_i\Phi _l(\mathbf {\xi }_i)$|⁠. Note that both matrices |$\mathbf {G}$| and |$\mathbf {\tilde{G}}$| depend only on the reference element and the basis functions and can be precomputed.

2.3.2 Nodal basis approach

As an alternative formulation, we consider switching from SeisSol’s modal to an nodal basis for incorporating plastic yielding. This approach provides polynomial coefficients representing the solution directly at a set of interpolation points that we may readily use to check the plastic yield criterion. To transform modal coefficients |$\hat{Q}_{lq}$| to nodal coefficients |$\sigma _{iq}^{\text{NB}}$| defined at L points |$\boldsymbol {\zeta }_1, \ldots , \boldsymbol {\zeta }_L$| we make use of a generalized Vandermonde matrix

\begin{equation*} \sigma _{iq}^{\text{NB}} = \sum _{l=1}^L \Phi _l(\boldsymbol {\zeta _i}) \hat{Q}_{lq} =: \sum _{l=1}^L\mathcal {V}_{il} \hat{Q}_{lq}, \end{equation*}

(20)

where |$\boldsymbol {\mathcal {V}}$| is an L × L matrix. In contrast to matrix |$\mathbf {G}$| in eq. (16), the Vandermonde matrix |$\boldsymbol {\mathcal {V}}$| is invertible, which allows to easily switch between different polynomial coefficients.

However, in order to ensure that computing the inverse of |$\boldsymbol {\mathcal {V}}$| is well-conditioned for a given set of basis functions, we use a special set of 3-D nodal points |$\mathcal {\boldsymbol \zeta }_i$| [after Warburton (2006) and Hesthaven & Warburton (2010), script available at https://github.com/tcew/nodal-dg/tree/master/Codes1.1/Codes3D]. The approach presented here is not restricted to this particular set of points but other choices might yield similar Lebesgue constants (Warburton 2006). The main advantage of these points here is that they can easily be constructed on runtime avoiding the storage in additional files. The distribution of nodal points |$\mathcal {\boldsymbol \zeta }_i$| is exemplarily visualized for p = 2 inside the reference tetrahedron in Fig. 2. The number of NB points n_NB = (p + 1)(p + 2)(p + 3)/6 conforms with the number of basis functions L for a given polynomial degree p.

In case of plastic yielding, we adjust the stresses at the NB points using eq. (13) and obtain |$\boldsymbol {\sigma }^\text{new}$|⁠. An L² projection to the polynomial space as in Section 2.3.1 is not required since |$\sigma _{iq}^{\text{new}}$| is already a polynomial. To reverse the adjusted nodal coefficients into the modal formulation, we use the inverse Vandermonde matrix |$\boldsymbol {\mathcal {V}}^{-1}$|⁠:

\begin{equation*} \hat{Q}_{lq}^{\text{new}} = \sum _{i=1}^{L} {\mathcal {V}}^{-1}_{li} \sigma _{iq}^{\text{new}}. \end{equation*}

(21)

The Vandermonde matrix |$\boldsymbol {\mathcal {V}}$| and its inverse can be precomputed.

In comparison to the previously discussed QP approach, the NB algorithm requires to check the plastic yielding criterion at considerably less subelement points, specifically for increasing polynomial degree p. Additionally, the points building the NB are spatially equally distributed, in distinction to the QPs that cluster towards element edges (see Fig. 2). Since these points are located at the element interfaces, their distance to the fault does not change with polynomial order. In contrast, Gaussian integration points are located closer to the fault surface for higher polynomial orders. Larger pointwise adjustments due to higher stresses experienced may be applicable. We demonstrate in Section 4 that neither the choice of implementation approach nor the location of the integration points impair the physical solution quality.

2.3.3 Computational efficiency

Both approaches presented above make use of a code generator for highly efficient matrix multiplications to compute eqs (16) and (19) or (20) and (21), respectively, in analogy to the optimization of viscoelastic rheologies in SeisSol (Uphoff & Bader 2016).

It is challenging to unambiguously define an overhead in terms of time-to-solution caused by incorporating plastic yielding, as it depends on the number of elements that yield, in both space and time. In addition, the performance characteristics of the target hardware architecture may shift the relative computational efficiency of the plasticity kernel compared to the elastic kernel. As a proxy metric, we provide lower and upper bounds for the computational overhead per element and time step of both approaches with respect to the number of floating point operations (flops). We expect this to be meaningful, because both kernels are generated using the same code generator.

As a lower bound estimate, we assume that the considered element is not yielding plastically: the yield condition is checked, but neither adjustment of the stress state nor mapping is required. The upper bound is based on the assumption of the element yielding plastically, therefore, requiring yield criterion evaluation, adjustment of the stress state and mapping into the original modal polynomial coefficients at each point inside this element.

In practice, only a small percentage of all elements, mostly situated in the vicinity of the fault, will experience plastic deformation in realistic dynamic rupture simulations (see e.g. Fig. 1). Furthermore, plastic yielding is restricted to the few time steps in which the rupture front and emitted waves are passing by. Therefore, the upper bound definition overestimates the computational cost considerably for application purposes.

We determine the computational overhead by counting the number of flops per element and time step in the elastic kernel as well as in the plasticity kernel. For matrix–matrix multiplications, flops are automatically determined by the code generator. For simple loops (e.g. scalar times vector), we manually count flops in the loop body. Note that the number of flops is hardware-specific (Uphoff & Bader 2016). Here, we analyse the Haswell architecture, which is used, for example, in the supercomputer SuperMUC, Phase 2 (https://www.lrz.de/services/compute/supermuc/). Furthermore, the number of flops in the elastic kernels has a slight dependence on the mesh generator, so we base the overhead on the average flops in the elastic kernel.

The estimated minimum overhead of the QP approach (Section 2.3.1) for varying polynomial degrees ranges between 23.9 per cent (p = 2) and 39.9 per cent, (p = 5), whereas the maximum increase of computational cost is estimated between 50.0 and 79.1 per cent. For the NB approach (Section 2.3.2), the lower bound is considerably lower: between 4.5 per cent (p = 2) and 6.6 per cent (p = 4); as well as the upper bound estimate that ranges between 8.6 and 13.1 per cent (Note that for NB, p = 4 has a larger relative overhead than p = 5). As a consequence, the NB approach is, in general, cheaper by up to a factor of 6.

Both approaches to incorporate plastic material response cause considerable computational overhead, as the number of QPs and the number of nodal coefficients grow with |$\mathcal {O}(p^3)$|⁠. However, the NB approach requires considerably less points per element (see Fig. 2); hence it is computationally more efficient at the cost of lower subelemental resolution. In Section 4 we analyse the accuracy of both approaches in a h- and p-refinement study of on-fault dynamics, shedding light onto the trade-off of computational cost versus accuracy.

For a typical dynamic rupture scenario (Southern California Eearthquake Center-U.S. Geological Survey (SCEC-USGS) benchmark problem TPV26/27 with plasticity, as presented in Section3.1), we measure a computational overhead of 12.5  per cent with the NB approach versus 46  per cent with the QP approach for polynomial degree p = 3. Computational performance (FLOP s⁻¹) is not affected by off-fault plasticity. In comparison, computational costs in the AWP-ODC code of Roten et al. (2016) using an FD method increases by 65  per cent for the same benchmark setup. Note that AWP-ODC is less computationally demanding with respect to time-to-solution for these type of benchmark problems.

3 VERIFICATION UNDER REALISTIC MODELLING ASSUMPTIONS

Two community benchmark problems established by the SCEC-USGS Spontaneous Rupture Code Verification Project [http://scecdata.usc.edu/cvws/, Harris et al. (2009, 2011, 2018)] combine off-fault plastic yielding with realistic fault geometries and heterogeneous initial conditions. In the following, we compare SeisSol to two state-of-the-art numerical methods in terms of synthetic ground motions and key dynamic rupture parameters for these two test problems.

In Section 3.1, we verify the implementation using a strike-slip fault as used in several previous studies with plasticity (e.g. Andrews 2005; Templeton & Rice 2008; Ma & Andrews 2010; Dunham et al.2011a). Additionally, asymmetric near-source ground motions induced by dipping faults (Oglesby et al.2000) are of particular interest and challenging in terms of mesh generation and numerical stability: High stresses are induced due to rupture-free-surface interaction that have to be resolved in relatively small elements and enhance plastic yielding in the hanging wall (Ma 2008, 2009).

The initial fault stress and strength conditions vary with depth in both benchmarks. Capturing these heterogeneities is crucial, since they may fundamentally affect rupture nucleation and propagation (Day 1982; Boatwright & Quin 1986; Oglesby & Day 2002; Ripperger et al.2007; Pelties et al.2014). To accurately capture heterogeneous fault initial conditions, we assign initial stress values to each of DG’s fault QPs, allowing subelement sampling (Pelties et al.2014).

We verify, in the following, the application of a high polynomial degree in conjunction with larger elements for resolving smaller-scale heterogeneous initial conditions for strike-slip and for dipping fault geometries. We additionally demonstrate the sensitivity of on-fault dynamic rupture measurements to the resolution of the nucleation zone (cf. Galis et al.2014). All results shown in this section are based on the QP plasticity implementation (Section 2.3.1). The QP and NB approaches are shown to yield near-identical results in Appendix  B.

3.1 Strike-slip fault benchmarks TPV26/27

Test problems TPV26 (elastic) and TPV27 (plasticity, with viscoplastic regularization) incorporate rupture on a vertical strike-slip fault geometry differing with respect to elastic versus plastic rheology assumptions. The fault frictional properties are governed by linear slip-weakening friction. For the viscoplasticity benchmark TPV27, initial stress loading is prescribed additionally off the fault throughout the domain.

Off-fault initial stresses need to be defined consistently with the stresses acting on-fault. Therefore, the commonly applied nucleation method of using an overstressed patch is not applicable for dynamic rupture simulations including plastic yielding. To initiate rupture, the fault is forced to break within a circular patch of radius r_crit around the hypocentre by gradually reducing the friction coefficient from its static to its dynamic value. In this manner, forced rupture is smoothly overtaken by spontaneous rupture which reduces numerical artefacts. The modelling parameters of both benchmarks are listed in Table 1.

Table 1.

Simulation parameters for SCEC benchmarks TPV26 (purely elastic) and TPV27 (with plasticity). The additional parametrization of viscoplasticity in TPV27 is denoted at the bottom.

Symbol	Parameter	Value
Vp	P-wave velocity	6000 m s−1
Vs	S-wave velocity	3464 m s−1
ρ	Density	2670 kg m−3
μs	Static friction coefficient	0.18
μd	Dynamic friction coefficient	0.12
Dc	Slip-weakening critical distance	0.3 m
d	Fault depth	0–20 000 m
|$w$|	Fault width	0–40 000 m
c0	Frictional cohesion
	d ≤ 5000 m	−4.0 MPa
		+0.00072 × d MPa m−1
	d > 5000 m	0.4 MPa
t0	Forced rupture decay time	0.5 s
rcrit	Nucleation patch radius	4000 m
Pf	Fluid pressure	9.8 × d MPa m−1
Ω(d)	Tapering coefficient
	d < 15 000 m	1
	15 000 m ≤d ≤ 20 000 m	(20 000 m − d)/5000 m
	d > 20 000 m	0
|$\sigma _n^0$|	Background normal stress	(−2.632.9267 × Ω(d)
		−26.166) × d MPa m−1
|$\tau _n^0$|	Background shear stress	−6.07897×
		Ω(d) × d MPa m−1
h	Smallest element edge length	250 m
p	Polynomial degree	4
	(space and time)
cplast	Plastic cohesion	1.36 MPa
ν	Bulk friction	0.1934
T|$v$|	Viscoplastic relaxation time	0.03 s

We here verify SeisSol’s results by comparing fault dynamics and ground motions to FaultMod, a well-established low-order FE software (Barall 2009). FaultMod can handle hexahedral or tetrahedral elements to represent 3-D geometries and material properties. Fault friction is implemented in FaultMod based on the traction-at-split-nodes method, featuring Newmark damping (Hughes 2000) and an optional thin viscous layer surrounding the fault zone (Day et al.2005; Dalguer & Day 2007) to suppress spurious high-frequency oscillations. Thus, comparison to the undamped ADER-DG approach is specifically of interest.

We compare SeisSol’s results of polynomial degree p = 4 and 250 m fault discretization to FaultMod’s results with 50 m element edge length and polynomial degree p = 2. All meshes are coarsened away from the fault by a gradation rate of 1.06 (Simmetrix Inc. 2017) in order to reduce computational cost. The unstructured mesh used for this benchmark is depicted in Fig. 1.

Fig. 3 illustrates the overall agreement of SeisSol and FaultMod in terms of on-fault rupture dynamics for both, the elastic problem and the benchmark accounting for off-fault plasticity. Comparing on-fault slip rates (Fig. 3a), we observe that plastic yielding significantly reduces peak slip rate (here, by approximately 50 per cent) and delays rupture arrival time (decreases rupture speed) as confirmed in various dynamic rupture models accounting for plastic energy dissipation (i.e. Andrews 2005; Dunham et al.2011a; Gabriel et al.2013; Roten et al.2015). In the elastic simulation, peak slip rate of SeisSol is slightly increased compared to FaultMod.

Figure 3.

(a) Along-strike slip rate, (b) along-strike shear and (c) normal stress at a fault location 10 km downdip and 10 km along-strike for the benchmark problems TPV26 [purely elastic (e), grey] and TPV27 [with plasticity (p), red]. Lighter colours denote the ADER-DG SeisSol solution (using the QP approach), darker colours the FEM FaultMod solution.

Remarkably, agreement improves distinctively when incorporating plastic material response. As we discuss in detail in Appendix  C, plasticity increases the cohesive zone width, the dynamic rupture problem inherent minimum length scale across which shear stress and slip rate vary abruptly (Day et al.2005). For dynamic rupture models that include plastic yielding a high resolution at the fault may be less crucial than commonly assumed: Analytical and numerical estimates based on purely elastic frameworks (in 2-D) may underestimate the minimum length scale that needs to be resolved in a plastically yielding simulation (Appendix  C). This also impacts numerical convergence behaviour as we discuss in the next section.

The excellent agreement between SeisSol and FaultMod in terms of synthetic ground motions is shown in Fig. 4 for the purely elastic and the plastically deforming problem. SeisSol captures all features of horizontal, vertical and normal ground velocity time-series even 20 km away from the fault where the mesh is already strongly coarsened to a maximum of 1300 m edge length. Within FaultMod’s grid-doubling approach, the resolution is doubled from 50 to 100 m in a distance of 10 element layers (t500 m) from the fault. SeisSol resembles off-fault velocities obtained by the second-order method FaultMod, by using larger elements with polynomial degree p = 4, which illustrates the potential strengths of a high-order method also for plastically yielding materials. For the simulations including plasticity, we overall observe smaller PGVs as a consequence of plastic deformation in the vicinity of the fault, consistent with recent studies (e.g. Roten et al.2014). Even 20 km away from the fault, the effect of the delayed rupture arrival in the model with plasticity (Fig. 3) is visible as time-shift of ground velocities (Fig. 4).

Figure 4.

(a) Horizontal, (b) vertical and (c) normal ground velocity time-series at a location 20 km off-fault (near side) perpendicular to the epicentre for the purely elastic (e) and plasticity (p) strike-slip fault benchmark TPV26/27. Lighter colours denote the ADER-DG SeisSol solution (using the QP approach), darker colours the FEM FaultMod solution.

3.2 Dipping fault benchmarks TPV12/13

Test problems TPV12 (elastic) and TPV13 (plasticity with viscoplastic regularization) incorporate rupture on a 60° dipping normal fault based on the Solitario Canyon Fault at Yucca Mountain (Harris et al.2009). All modelling parameters are summarized in Table 2. Rupture is initiated by prescribing a lower static coefficient of friction inside a predefined rectangular patch surrounding the hypocentre. In that way, the initial shear stress is greater than the yield stress and rupture initiates immediately after the simulations starts. Linear-slip weakening friction is assumed. The asymmetric unstructured tetrahedral mesh of the here presented solution of SeisSol exhibits a minimal element edge length of 250 m across the fault surface; its volume discretization is coarsened by a gradation rate of 1.06 with increasing fault distance (Simmetrix Inc. 2017).

Table 2.

Simulation parameters for SCEC benchmarks TPV12 (purely elastic) and TPV13 (with plasticity). The additional parametrization of viscoplasticity in TPV13 is denoted at the bottom.

Symbol	Parameter	Value
Vp	P-wave velocity	5716 m s−1
Vs	S-wave velocity	3300 m s−1
ρ	Density	2700 kg m−3
μs	Static friction coefficient	0.7
μs, nuc	Static friction coefficient	0.54
	(inside nucleation patch)
μd	Dynamic friction coefficient	0.1
Dc	Slip-weakening critical distance	0.5 m
d	Fault depth	0–15000 m
|$w$|	Fault width	0–30000 m
c0	Frictional cohesion	−0.2 MPa
Pf	Fluid pressure	9.8 × d MPa m−1
|$\sigma _n^0$|	Background normal stress
	(d <13800 m)	−7.39001 × d MPa m−1
	(d ≥13800 m)	−14.42798 × d MPa m−1
τ0	Background shear stress
	(d <13800 m)	−0.549847|$\cdot \sigma _n^{0}$|
	(d ≥13800 m)	0.0
Anuc	Nucleation size	3 km × 3 km
h	Smallest element edge length	250 m
p	Polynomial degree in space and time	4
cplast	Plastic cohesion	5.0 MPa
ν	Bulk friction	0.85
T|$v$|	Viscoplastic relaxation time	0.03 s

In Fig. 5 we show the overall very good agreement of SeisSol with the solution of FaultMod. We additionally compare our solution with the high-order spectral element (SE) method SPECFEM3D (Kaneko et al.2008), which is based on unstructured hexahedral meshes (Peter et al.2011) and uses a traction-at-split-nodes approach to incorporate fault dynamics. As in FaultMod, SPECFEM3D damps spurious oscillations by a thin layer of 1 or 2 Kelvin-Voigt elements surrounding the fault (Galvez et al.2014). The fault discretization is 100 m for the FaultMod and 100 m for the SPECFEM3D solution, the latter with additional four subelemental integration points in each dimension.

Figure 5.

(a) Along-strike slip rate, (b) along-strike shear and (c) normal stress at a fault location 1.5 km downdip aligned with the hypocentre for benchmark TPV13. The ADER-DG SeisSol solution (using the QP approach) is shown in black, the FaultMod solution in grey, the SPECFEM3D solution in orange.

Slip rate, shear stress, and normal stress (Fig. 5) are in equally good agreement. However, both, FaultMod and SPECFEM3D solutions, exhibit high-frequency oscillations most visible in the (relatively small) normal stress amplitudes, despite the implemented damping. Furthermore, we note a small time-shift in the arrival of the peak amplitudes for all three values of both high-order methods compared to FaultMod. Rupture arrival time is a sensitive indicator of numerical precision (Day et al.2005). Both, SeisSol and SPECFEM3D reach subelement resolution of the nucleation patch causing a slightly faster rupture initiation in comparison to the low-order method FaultMod.

4 MESH AND POLYNOMIAL DEGREE REFINEMENT STUDY

The comparison of numerical methods in well-defined benchmark problems is a commonly used approach to verify dynamic (spontaneous) earthquake rupture implementations (e.g. Pelties et al.2014; Harris et al.2018;). However, the convergence of the numerical solution with increasing mesh resolution (h-refinement) and increasing polynomial degree (p-refinement) cannot be guaranteed in such a manner.

A formal convergence analysis requires the analytical solution of the underlying problem. Analytical solutions to dynamic rupture problems have been proposed for a 2-D self-similar crack (Kostrov 1964) or self-similar pulse (Nielsen & Madariaga 2003), but are not available for typical application-based dynamic rupture simulations. Hence, a high-resolution reference solution that is invariant under additional discretization and polynomial refinement is commonly defined as stand in for the exact solution (e.g. Day et al.2005; Rojas et al.2008; Pelties et al.2012). In this manner, we provide helpful guidelines for the required fault resolution guaranteeing accurate on-fault physical results in the framework of dynamic rupture problems. In the following, the terms convergence and convergence rate refer to the given reference solution. We critically discuss this approach in Section 6.2.

Few refinement studies of this kind for dynamic rupture simulations have been presented in 2-D (Kaneko et al.2008; Rojas et al.2008; Rojas et al.2009; De la Puente et al.2009; Huang & Ampuero 2011) and 3-D (Day et al.2005; Pelties et al.2012; Tago et al.2012). All of them were based on idealized setups assuming homogeneous initial stresses, purely elastic material behaviour and abrupt nucleation procedures. However, state-of-the-art dynamic rupture simulations emulate real fault properties using depth-dependent initial stresses (e.g. Ma & Andrews 2010; Kozdon & Dunham 2013; Heinecke et al.2014), off-fault plasticity (e.g. Roten et al.2015, 2016) and smoother nucleation procedures (e.g. Bizzarri 2010; Roten et al.2017).

Such models challenge convergence studies: For example, the assignment of initial fault stress and strength, as well as of the initial bulk stress, depends on the underlying discretization. Mesh refinement results in a higher sampling of the initial fields and smaller discontinuities in between neighbouring elements. Furthermore, the onset of rupture may be affected even by slight variations in the sampling of the nucleation conditions (Galis et al.2014).

Here, we present an in-depth h- and p-refinement study of a 3-D dynamic rupture setup including now common modelling complexities: We combine elastic as well as plastic material response, following both plasticity approaches presented in Section 2.3, with depth-dependent initial stress conditions. To ensure consistency of the initial stress tensor in the bulk and on the fault, we employ a nucleation procedure smooth in space and time in contrast to overstressed nucleation patches employed for previous convergence studies.

This section aims to provide hands-on guidelines in terms of the required resolution to ensure the accurate resolution of on-fault rupture dynamics of realistic, 3-D dynamic rupture applications. To this end, we numerically investigate the convergence of key dynamic rupture parameters with mesh refinement and increasing polynomial order towards a reference solution. Furthermore, we analyse the size of the cohesive zone as the inherent length scale determining accurate resolution.

4.1 Models and procedure

We measure the root-mean square (RMS) error of dynamic rupture characteristics in comparison to a reference solution. The reference solution is generated using a small fault discretization h and high polynomial degree p (e.g. Day et al.2005; Pelties et al.2012). We ensure that the reference solution does not change with further mesh or polynomial refinement, that is, accurately resolves the cohesive zone width as described in the next section.

In particular, we analyse rupture arrival time, peak slip rate amplitude, peak slip rate time and the total amount of slip (final slip) on the fault surface. The numerical accuracy in resolving time-dependent parameters is particularly important for complex source dynamics involving rupture jumps and branching. Final slip represents an integrated value determined by rupture arrival time, peak slip rate amplitude and time, and rise time across the fault.

We base the following refinement study on SCEC benchmark problems TPV26 (elastic) and TPV27 (with plasticity; see Section3.1). To reduce the computational cost, we shorten the simulation time to 6.3 s and shrink the fault surface to 21 km × 30 km. The hypocentre is shifted to the centre of the new fault domain (at depth of 10.5 km) and we slightly decrease the size of the circular nucleation patch (r_crit =3.5 km) to account for the smaller fault area. All remaining modelling parameters, specifically the depth-dependence of initial stresses, are equivalent to Table 1.

We perform simulations with increasing polynomial degree of basis functions from p = 2 to p = 5 and decreasing fault mesh discretization from h = 1061 m to h = 106 m. The high-resolution reference solution is employing p = 5 and h = 71 m.

To quantify the accuracy of our results with respect to the reference solution, we calculate the RMS errors at a total number of 362 on-fault locations (‘receivers’). The RMS error is then normalized by the mean over all receivers of the reference solution, ensuring comparability with previous studies (Day et al.2005; Pelties et al.2012). The receivers are located at all fault gridpoints of the coarsest mesh (h = 1061 m). For our analysis, we excluded the nucleation patch, receivers within 1 km from the fault edges and areas that rupture does not reach within the simulated time. For each simulation, the same points are evaluated independent of the mesh discretization and polynomial degree.

We emphasize that pointwise measures are prone to large variations for any Galerkin-type method since the equations are solved in an integral (weak) sense, with the largest variations being expected at the mesh vertices. Thus, our method is analysing the ‘worst-case’ solution evaluation. However, pointwise evaluation of rupture dynamics characteristics is commonly used for the interpretation of simulations.

We generate a set of fully regular triangles discretizing the fault surface (except the circular nucleation patch). This is achieved by first dividing the fault surface into quadrilateral elements of edge lengths |$h/\sqrt{2} =$| 750, 375, 250, 125, 75 and 50 m. Then, each quadrilateral is subdivided into two triangles, resulting in fault discretizations of edge lengths h = 1061, 530, 354, 177, 106 and 71 m. Thus, the mesh vertices of the coarsest discretization are a subset of the gridpoints of the finest discretization. In contrast, the volume mesh is fully unstructured, coarsened by a factor of 10 per cent away from the fault up to a maximum edge length of 10h. The regular fault discretization does not affect the generality of the results but facilitates evaluating exactly the same receiver locations (Pelties et al.2012).

All meshes for the refinement tests are generated with the open-source software gmsh (Geuzaine & Remacle 2009). All simulations of this section were conducted using SeisSol (https://github.com/SeisSol) with git hash 72596ff.

4.2 Resolving the cohesive zone in 3-D simulations with heterogeneous fault initial conditions

Dynamic rupture simulations have to accurately resolve a problem inherent length scale, the cohesive zone width Λ that spans the part of the fault across which shear stress decreases from its static to its dynamic value. The width of the cohesive zone is decreasing with increasing rupture velocity (e.g. Day et al.2005; Duan & Day 2008).

We present a detailed analysis of the dynamic evolution of the cohesive zone width, along-strike and up-dip, under heterogeneous initial stress, and with and without off-fault plasticity in Appendix  C. We find that the cohesive zone width varies considerably across the fault in dependence of the rupture speed that is determined by the depth-dependent initial stress, frictional properties and propagation distance. Absolute rupture speed reaches a maximum of 3250 m s⁻¹ in the purely elastic reference simulation, whereas rupture speed is limited to 3050 m s⁻¹ in the reference simulation with plasticity. Both values are well below the theoretical terminal speed (p-wave speed) for 3-D mixed mode rupture propagation (Bizzarri & Das 2012) and below the s-wave speed of 3464 m s⁻¹. Specifically close to the free surface, the increase of cohesion results in rising frictional resistance and the decrease of initial stresses results in very low rupture speeds leading to large cohesive zone widths. In general, off-fault plasticity leads to a distinctively larger cohesive zone width than under comparable purely elastic conditions since rupture propagates at lower speeds.

The minimum cohesive zone width for both setups, purely elastic and plastic material behaviour, is located at the largest distance from the hypocentre along-strike at a depth of 14.5 km where the absolute value of the depth-dependent initial stresses are highest. For the purely elastic case, the minimum cohesive zone with |$\Lambda _{\mathrm{ min}}^\mathrm{ e}$| is 162 m, which is considerably smaller than in previous convergence setups [Λ_min = 325 m in Day et al. (2005)]. In the simulation with off-fault plasticity, this doubles to a width of |$\Lambda _{\mathrm{ min}}^\mathrm{ p}$| = 325 m. Consistent with previous studies (Day et al.2005; Pelties et al.2012; Tago et al.2012), we also evaluate the median value of the cohesive zone width over all fault receivers: |$\bar{\Lambda }^\mathrm{ e}$| is 583 m and |$\bar{\Lambda }^\mathrm{ p}$| is 722 m.

We follow Day et al. (2005) in defining a solution to be sufficiently close to the reference solution once the RMS errors reached the following thresholds: lower than 0.2  per cent for rupture arrival time, lower than 7  per cent for peak slip rate and lower than 1  per cent for final slip. In order to accurately resolve dynamic rupture characteristics (i.e. fulfilling the conditions described above), we determine the minimum required mesh resolution at the fault for a given polynomial degree p. The such-defined mesh resolution may differ for different numerical methods. Assessing this resolution requires careful analysis in the form of self-consistent refinement tests (see also Section 6).

The number of elements per median cohesive zone width are denoted by N_c (e.g. Day et al.2005). Since the cohesive zone width changes dynamically with rupture speed, also N_c varies across the fault plane for a given discretization.

In the purely elastic convergence setup, we find N_c varying between 0.55 and 5.5 elements per median cohesive zone width |$\bar{\Lambda }$| for discretizations between 1061 and 106 m. In the simulations including plasticity, N_c falls in between 0.69 and 6.91 elements.

SeisSol resolves shear and normal stress at the fault at (p + 2)² Gaussian QPs inside each fault triangle. Thus, the longest edge length h is itself discretized by p + 1 additional points (see Fig. C1 in Appendix  C). As a result, the effective fault discretization is smaller than denoted by the grid spacing h and a higher resolution than denoted by N_c is achieved. We additionally introduce the number of subelement points per median cohesive zone width |$N_\mathrm{ c}^{\mathrm{ sub}} = N_\mathrm{ c}(p+1)$|⁠. For example, in a p = 5 simulation with N_c = 0.55 elements, the median cohesive zone width is effectively resolved by |$N_\mathrm{ c}^{\mathrm{ sub}}$| = 3.3 points. However, we will still present N_c as guidance to choose the corresponding required minimum mesh discretization. This approach also allows direct comparison to low-order methods (e.g. Day et al.2005).

The minimum required fault mesh resolution h for each polynomial degree p based on purely elastic convergence tests and including plastic yielding is reported in Sections 4.3 and 4.4, respectively.

4.3 Convergence of 3-D elastic dynamic rupture simulations with heterogeneous initial conditions

4.3.1 Absolute errors of on-fault measurements

The absolute RMS errors measured in the elastic refinement study are summarized in Table 3 and visualized in Fig. 6. The errors of the rupture arrival time, the peak slip rate time, the final slip and the peak slip rate with respect to the reference solution are globally decreasing with mesh refinement and increasing polynomial degree. The smooth convergence of on-fault quantities towards the reference solution illustrates the robustness of the numerical method in handling heterogeneous initial conditions.

Figure 6.

Refinement study results for purely elastic material behaviour and heterogeneous initial conditions. RMS errors with respect to the reference solution are shown for (a) rupture arrival time, (b) final slip, (c) peak slip rate time and (d) peak slip rate. Markers denote the RMS errors in  per cent, normalized by the mean over all receivers of the reference solution (h = 71 m, p = 5). Different colours and markers represent varying polynomial degrees. N_c denotes the number of elements per median cohesive zone for a given fault mesh spacing. The solid lines represent least-squares fits.

Table 3.

RMS errors of elastic refinement tests with varying discretization h and polynomial degree p. The bottomline denotes the mean over all receivers of the reference solution.

h (m)	p	Rupture arrival time ( per cent)	Final slip (per cent)	Peak slip rate time ( per cent)	Peak slip rate (per cent)
1061	2	1.42	3.39	1.79	47.03
	3	0.60	2.87	0.89	35.10
	4	0.34	2.12	0.52	27.83
	5	0.24	2.03	0.45	22.58
530	2	0.44	1.56	0.64	30.63
	3	0.25	1.32	0.28	21.32
	4	0.18	0.99	0.24	16.04
	5	0.12	0.94	0.23	12.23
354	2	0.26	1.13	0.35	23.44
	3	0.15	0.95	0.18	15.54
	4	0.09	0.74	0.10	10.22
	5	0.08	0.69	0.11	6.57
212	2	0.12	0.63	0.16	15.61
	3	0.06	0.54	0.07	8.50
	4	0.04	0.41	0.07	4.48
	5	0.04	0.40	0.05	2.57
106	2	0.04	0.32	0.05	6.64
	3	0.03	0.28	0.03	2.60
	4	0.02	0.22	0.03	1.38
	5	0.02	0.21	0.03	1.04
71	5	3.90 s	1.67 m	4.03 s	2.98 m s−1

We note that the minimum measured RMS error of rupture arrival time is only 0.02  per cent (h = 106 m and p = 4, 5). Even for the coarsest discretization and lowest polynomial degree (h = 1061 m and p = 2), the solution for rupture arrival time exhibits an error below 1.5 per cent with respect to the reference solution. Equally low RMS errors can be found for peak slip rate time. In contrast, the RMS errors measured for the peak slip rate reach up to 47  per cent for the same polynomial degree and fault resolution. The RMS errors for final slip ranges between 3.39 per cent (p = 2, h = 1061 m) and 0.21 per cent (p = 5, h = 106 m). As an integrated value, the RMS errors for final slip reflects the low RMS errors of the time-dependent quantities in conjunction with the higher RMS errors of the peak slip rate.

Coarse (h- and p-) resolution highly impacts peak slip rate, which directly affects ground motion assessment. This suggests that the measured peak value of slip rate is very sensitive to the sampling of the depth-dependent initial conditions and hence to the resolution of the fault. In distinction, the overall lower RMS errors of rupture arrival time, peak slip rate time and final slip in our setup indicate that temporal and integrated fault characteristics are more robust: lower resolution of the cohesive zone still results in acceptable small RMS errors.

Using polynomial degrees p = 2, 3, solutions are sufficiently close to the reference solution for a minimum resolution of h = 106 m which corresponds to N_c = 5.5 (⁠|$N_\mathrm{ c}^{\mathrm{ min}}$| = 1.53, |$N_\mathrm{ c}^{\mathrm{ sub}}$| = 16.5–22.0). For p = 4 a resolution higher than 212 m (N_c = 2.75, |$N_\mathrm{ c}^{\mathrm{ min}}$| = 0.76, |$N_\mathrm{ c}^{\mathrm{ sub}}$| = 13.75) and for p = 5 a resolution higher than 354 m (N_c = 1.65, |$N_\mathrm{ c}^{\mathrm{ min}}$| = 0.46, |$N_\mathrm{ c}^{\mathrm{ sub}}$| = 9.9) result in sufficiently small RMS errors. We note that the convergence criteria for rupture arrival time and final slip are met already at much coarser discretizations. However, the RMS error of the peak slip rate is still too high for fulfilling all conditions for convergence towards the reference solution. In general, we state that the median cohesive zone width |$\bar{\Lambda }$| needs to be resolved by approximately 5–6 (p = 2, 3), 2–3 (p = 4) and 1–2 (p = 5) elements, or 22 (p = 2), 16.5 (p = 3), 13.75 (p = 4) and 9.9 (p = 5) subelement points for purely elastic simulations featuring depth-dependent initial stress conditions.

4.3.2 Low-order convergence rates of on-fault measurements

The corresponding on-fault convergence rates are defined as the slope of a least-squares fit to the RMS values and are given in Table 4. The reported rates are defined with respect to the pre-defined reference solution and are therefore only representative of the formal convergence of the numerical scheme, if the reference solution is sufficiently close to the exact solution (which is unknown). For the following discussion, supplemented by the cohesive zone size analysis in Appendix  C, we assume that our high-resolution reference solution fulfils this requirement. Even though we observe clear h- and p-convergence towards the reference solution, the following specific characteristics of the convergence rates require special consideration:

Table 4.

RMS error convergence rates with respect to the reference solution for varying polynomial degree p for the refinement studies with elasticity and plasticity.

	p	Rupture arrival time	Final slip	Peak slip rate time	Peak slip rate
Elastic	2	1.53	1.02	1.52	0.84
	3	1.38	1.01	1.43	1.12
	4	1.26	0.98	1.26	1.32
	5	1.17	0.98	1.27	1.39
Plastic (QP)	2	1.67	1.13	1.74	0.68
	3	1.40	1.05	1.37	0.83
	4	1.33	1.01	1.06	0.91
	5	1.23	0.99	1.01	0.92
Plastic (NB)	2	1.57	1.31	1.70	0.73
	3	1.16	1.07	1.26	0.85
	4	1.31	1.03	1.04	0.92
	5	1.28	1.01	1.04	0.93

Overall, no spectral convergence is achieved for on-fault dynamic rupture quantities. In particular, the reported convergence rates, ranging between 1.53 and 0.84, indicate that there is no high-order (≥2) convergence for any of the dynamic rupture observables. Additionally, only peak slip rate exhibits increasing rates with higher polynomial order. Convergence rates for rupture arrival time, peak slip rate time and slip do not clearly increase with increasing p.

Low-order (at most quadratic) convergence rates of dynamic rupture on-fault parameters have been previously reported (Kaneko et al.2008; Rojas et al.2008; Rojas et al.2009; Pelties et al.2012). Rojas et al. (2008), Kaneko et al. (2008) and Pelties et al. (2012) show that on-fault convergence rates achieved by high-order (‘spectral’) numerical methods do not exceed the reported rates of second-order methods (Day et al.2005) for this type of refinement study.

A notable exception is posed by the convergence studies presented in Kozdon & Dunham (2013), O’Reilly et al. (2015), Duru & Dunham (2016) and Erickson et al. (2017): High-order convergence is achieved with respect to an a priori formulated analytical solution of a related but simplified problem [method of manufactured solutions (MMS)]. We discuss a potential general limitation of convergence rates for dynamic rupture problems with respect to mathematical and numerical theory in Section 6.2.

4.4 Convergence of 3-D dynamic rupture simulations with off-fault plasticity and heterogeneous initial conditions

We extend the elastic refinement study to account for plastic energy dissipation. Both of the implementation schemes presented in Section 2.3 are analysed. The RMS errors of the QP-based implementation approach are summarized in Table 5 and visualized in Fig. 7. The RMS errors of the NB-based implementation approach is given in Table 6 and Fig. 8. On-fault convergence rates (slopes) with respect to the reference solution are listed in Table 4.

Figure 7.

Refinement study results for plasticity implemented via the QP approach. RMS errors with respect to the reference solution are shown for (a) rupture arrival time, (b) final slip, (c) peak slip rate time and (d) peak slip rate. Markers denote the RMS error in  per cent, normalized by the mean over all receivers of the reference solution (h= 71 m, p = 5). Different colours and markers represent varying polynomial degrees. N_c denotes the number of elements per median cohesive zone for a given mesh spacing. The solid lines represent least-squares fits.

Figure 8.

Refinement study results for plasticity implemented via the NB approach. RMS errors with respect to the reference solution are shown for (a) rupture arrival time, (b) final slip, (c) peak slip rate time and (d) peak slip rate. Markers denote the RMS error in  per cent, normalized by the mean over all receivers of the reference solution (h= 71 m, p = 5). Different colours and markers represent varying polynomial degrees. N_c denotes the number of elements per median cohesive zone for a given mesh spacing. The solid lines represent least-squares fits.

Table 5.

RMS errors of QP plasticity refinement tests with varying discretization h and polynomial degree p. The bottomline denotes the mean over all receivers of the reference solution.

h (m)	p	Rupture arrival time ( per cent)	Final slip ( per cent)	Peak slip rate time ( per cent)	Peak slip rate (per cent)
1061	2	2.01	4.94	2.48	10.34
	3	0.65	3.44	0.85	9.56
	4	0.38	2.45	0.46	9.12
	5	0.29	2.26	0.40	8.46
530	2	0.40	1.88	0.57	6.29
	3	0.20	1.52	0.27	6.89
	4	0.15	1.10	0.16	5.29
	5	0.13	1.03	0.14	4.06
354	2	0.22	1.33	0.28	5.63
	3	0.12	1.07	0.14	4.96
	4	0.09	0.81	0.11	3.66
	5	0.07	0.75	0.09	2.91
212	2	0.10	0.71	0.12	3.75
	3	0.06	0.59	0.08	2.80
	4	0.04	0.44	0.07	2.15
	5	0.03	0.42	0.06	1.77
106	2	0.04	0.36	0.04	2.06
	3	0.03	0.31	0.04	1.49
	4	0.02	0.23	0.04	1.16
	5	0.02	0.23	0.04	1.00
71	5	4.18 s	1.50 m	4.36 s	1.98 m s−1

Table 6.

RMS errors of NB plasticity refinement tests with varying discretization h and polynomial degree p. The bottomline denotes the mean over all receivers of the reference solution.

h (m)	p	Rupture arrival time ( per cent)	Final slip ( per cent)	Peak slip rate time ( per cent)	Peak slip rate ( per cent)
1061	2	3.54	8.03	3.92	11.26
	3	0.64	3.70	0.83	10.01
	4	0.36	2.46	0.40	9.24
	5	0.31	2.32	0.42	8.58
530	2	0.45	2.09	0.61	6.94
	3	0.28	1.57	0.33	6.89
	4	0.27	1.17	0.22	5.51
	5	0.22	1.08	0.17	4.16
354	2	0.32	1.42	0.35	5.83
	3	0.25	1.14	0.21	5.13
	4	0.16	0.86	0.14	3.71
	5	0.12	0.77	0.11	2.91
212	2	0.19	0.75	0.18	3.84
	3	0.13	0.61	0.12	2.81
	4	0.07	0.45	0.08	2.15
	5	0.05	0.42	0.06	1.76
106	2	0.07	0.36	0.07	2.06
	3	0.04	0.31	0.04	1.49
	4	0.02	0.24	0.04	1.17
	5	0.02	0.23	0.04	1.00
71	5	4.18 s	1.51 m	4.35 s	1.98 m s−1

For both approaches and all measured quantities, we observe h and pconvergence. The observed smooth convergence of all measured on-fault quantities towards the reference solution illustrates the robustness of the numerical method while stresses are adjusted due to plastic yielding.

The RMS errors of rupture arrival time and peak slip rate time become very low (0.02 and 0.04 per cent, respectively) at high resolutions for both plasticity implementations. Overall, RMS errors do not exceed 10 per cent even at the coarsest discretization and smallest polynomial degree under consideration. It is notable that polynomial degrees 4 and 5 yield in general very similar absolute errors for both approaches, specifically at higher resolutions.

The required discretization to accurately resolve the cohesive zone width with plasticity differs to the purely elastic case. For the QP approach, our results are sufficiently close to the reference solution at a minimum mesh resolution of h = 212 m (N_c = 3.41, |$N_\mathrm{ c}^{\mathrm{ min}}$| = 1.53, |$N_\mathrm{ c}^{\mathrm{ sub}}$| = 10.23) for polynomial degree p = 2, h = 354 m (N_c = 2.04, |$N_\mathrm{ c}^{\mathrm{ min}}$| = 0.92, |$N_\mathrm{ c}^{\mathrm{ sub}}$| = 8.16–10.2 ) for p = 3, 4 and approximately h = 530 m (N_c = 1.36, |$N_\mathrm{ c}^{\mathrm{ min}}$| = 0.61, |$N_\mathrm{ c}^{\mathrm{ sub}}$| = 8.16) for p = 5. The NB approach requires similar minimum on-fault resolutions, besides for p = 3 which requires a higher resolution of h = 212 m (N_c = 3.41, |$N_\mathrm{ c}^{\mathrm{ min}}$| = 1.53, |$N_c^{sub}$| = 13.64) than the corresponding polynomial degree in the QP approach.

Both plasticity implementations show non-spectral convergence for all on-fault quantities with respect to the reference solution, similar to the elastic refinement study in Section 4.3 (see also Section 6.2).

4.4.1 Accuracy versus efficiency for off-fault plasticity implementation schemes

We shed light on the trade-off between accuracy and efficiency based on the refinement study of the QP and NB plasticity implementations. We aim to enable conclusions on a preferable scheme for a given dynamic rupture problem at hand. The reference solutions of both schemes (last row in Tables 5 and 6), give near-identical results in terms of key dynamic rupture characteristics. As a consequence, the normalized RMS differences are well comparable.

The RMS errors of all on-fault quantities are in general smaller using the QP approach than the NB implementation. The largest difference can be found for p = 2 and h = 1061 m. The RMS error for rupture arrival is here 23 per cent larger in the NB approach than in the QP approach (3.54 versus 2.01 per cent) since the QP approach approximates the nonlinearity more accurately with respect to the L² norm which is, in particular, pronounced for low fault resolutions and low polynomial degrees. RMS errors for p > 3 and for a median cohesive zone width resolution (N_c) of more than 2.04 elements are near equivalent for both plasticity implementations.

The convergence rates with respect to the reference solution reported in Table 4 are similar for all quantities in both approaches, in particular for polynomial degrees p = 4 and 5. The convergence rates of the time-dependent quantities decrease with increasing polynomial degree p. The slope of the least-squares fit to the RMS values of the peak slip rate is increasing with increasing p, while it saturates for both implementations at 1.0 for the final slip rate.

In terms of computational costs, the QP plasticity refinement study simulations are 24–79  per cent more expensive than the corresponding purely elastic simulations, depending on polynomial degree p and fault discretization length h. The NB plasticity approach is computationally distinctively less expensive, resulting in an increase of only 4.5–13 per cent compared to the corresponding elastic models.

We conclude that it is beneficial to use the QP plasticity approach for low mesh resolutions (N_c < 1.36 ) and low polynomial degree (p = 2). In this case, it yields higher accuracy of on-fault dynamics. However, in all other cases, the NB plasticity approach pays off, that is, it is computationally less expensive yielding comparably small RMS errors.

4.4.2 Effect of off-fault plasticity on dynamic rupture convergence

Fig. 9 compares slip rate over time for varying mesh discretizations and polynomial degrees for the refinement setups with elasticity and the two plasticity implementations. The figure highlights how recorded slip rates vary distinctively in terms of timing, peak and encapsulated slip.

Figure 9.

Exemplary slip rate measurements in the elastic (red) and plastic (orange QP approach, black NB approach) refinement tests. The time-series are recorded at 11 250 m along-strike and at 7500 m depth (location visualized by a triangle in Fig. C1 in Appendix  C). Three different mesh discretizations and polynomial degrees are marked by different line styles. Numbers indicate the peak slip rate of the corresponding curve. The dashed (530 m, p = 3) and solid lines (71 m, p = 5) align almost perfectly for both plasticity implementations.

Time-dependent RMS errors, such as rupture arrival time and peak slip rate time, are significantly larger for coarse mesh discretizations (h = 1061 m) and low polynomial degrees (p = 2, 3) for the simulations with plasticity in comparison to the elastic results (Tables 3, 5, 6 and visualized in Fig. 9). In the models with pure elasticity rupture, arrival time is already very consistent with the reference solution even for large on-fault discretization.

Remarkable are the small absolute RMS errors with respect to the reference solution for peak slip rate in the plastically yielding simulations: For all simulations with plasticity, the RMS errors fall in between 1 and 11  per cent (Tables 5 and 6) whereas the RMS error in the elastic case reach up to 47  per cent (Table 3). Fig. 9 exemplarily visualizes the high variation of peak slip rate in the elastic case (up to 20  per cent) while differences with plasticity are below 7 per cent even for the coarsest resolution.

With plasticity, rupture speed is decreased by off-fault energy absorption compared to elastic settings. The cohesive zone width does not shrink as pronouncedly as in elastic simulations with propagation along-strike (see Appendix  C and Andrews 2005). As a result, the median cohesive zone width |$\bar{\Lambda }$| is larger by up to a factor of 2 in the models with plastic yielding . Thus, |$\bar{\Lambda }$| is naturally better resolved in the simulations with plasticity than in elastic models at comparable mesh discretization h and polynomial degree p.

Analysing the median cohesive zone width resolution N_c, which is independent of Λ, we find that less elements per cohesive zone width are required with plasticity to gain comparably small RMS errors of peak slip rate. This parameter is the most sensitive dynamic rupture characteristic in our analysis. Smaller RMS errors for peak slip rate in the models including plasticity can, therefore, be explained by a better resolution of Λ, plus a general regularizing effect of plasticity (smoother reference solution).

In contrast, we observe larger RMS errors for the time-dependent quantities (rupture arrival time, peak slip rate time) with plasticity than at equivalent N_c in the elastic case. This suggests that time-dependent measurements are considerably impacted by (mesh-dependent) off-fault plastic energy dissipation. A higher resolution is required to resolve these quantities with the same accuracy as in a comparable elastic simulation.

In summary, the overall minimum required resolution of the cohesive zone width is 34.4–50 per cent lower with plasticity due to the lower resolution required for the sensitive peak slip rate.

5 LANDERS FAULT SYSTEM SCENARIO WITH OFF-FAULT PLASTICITY

We now aim to demonstrate the specific advantages of the modelling framework for large-scale earthquake simulations including various representations of natural complexity. To this end, a dynamic rupture scenario of the branched fault system hosting the 1992 M_w 7.3 Landers earthquake is presented. We analyse source dynamics at 3-D multisegment curved faults with and without plastic material response in terms of rupture transfers via dynamic triggering (‘jumping’) and direct branching. Furthermore, we compare peak slip rate, rupture arrival and total slip distribution.

The Landers earthquake is a prominent example of an earthquake rupturing across a geometrically complex fault network. The cascading event activated at least five fault segments overlapping over approximately 80 km, demonstrating a surprising interconnectivity. The Landers event raised awareness of unexpectedly large magnitude earthquakes enabled by rupture transfer mechanisms still not completely understood.

Based on high-quality observational data, the Landers earthquake has been investigated in great detail: Its slip distribution has been inferred from inversion of seismological and geodetic data (e.g. Wald & Heaton 1994; Cotton & Campillo 1995; Fialko 2004; Xu et al.2016), as well as from geological perspectives (e.g. Wesnousky 2006; Madden & Pollard 2012; Madden et al.2013) and in dynamic rupture simulations (e.g. Peyrat et al.2001; Aochi & Fukuyama 2002; Fliss et al.2005; Tago et al.2012). All these studies assume purely elastic behaviour of the material surrounding the fault.

Including off-fault plasticity in dynamic rupture simulations on a single planar fault revealed that the reduction of slip at shallow depth (‘shallow slip deficit’), which has been observed in the Landers fault system amongst others (Fialko et al.2005), can be partially related to near-surface plastic deformation (Kaneko & Fialko 2011; Roten et al.2017). Accounting for fault complexity in conjunction with off-fault plasticity is expected to alter fault dynamics with respect to stress transfer and dynamic triggering of different segments.

We point out that including off-fault plasticity (using the NB approach) only adds 7 per cent of computational cost compared to a fully elastic production run. In the following, we restrict ourselves to the analysis of on-fault earthquake dynamics. A detailed analysis of induced ground motions and other off-fault observables is beyond the scope of this paper and will be addressed in future work.

5.1 Model setup

We base the dynamic rupture scenario on the fully elastic model presented in Heinecke et al. (2014), which has proven excellent scalability on some of the largest supercomputers worldwide. This allows the analysis of high-resolution rupture behaviour at full geometrical complexity directly coupled to seismic wave propagation. The model includes fault traces by Fleming et al. (1998; white lines in Fig. 11) extended to 16 km in depth. The computational mesh accounts for regional topography and a 1-D subsurface velocity structure, adapted from Graves & Pitarka (2010).

Previous dynamic rupture simulations of the Landers earthquake have failed to dynamically interconnect all fault segments (Pelties et al.2012; Tago et al.2012) or required initial stress or strength heterogeneities, such as principal stress angle rotation or weak step-over faults (e.g. Peyrat et al.2001; Aochi & Fukuyama 2002).

We here prescribe lateral homogeneous, but depth-dependent initial stress conditions throughout the modelling domain, in contrast to previous dynamic rupture scenarios assuming varying principal stress angles (e.g. Aochi et al.2003; Galvez et al.2014, 2016). We assume a maximum principal stress orientation of 11° North for the entire fault system, which has been found to optimally reproduce slip in geomechanical modelling (Madden et al.2013). The lateral homogeneous stress generates a locally heterogeneous stress field due to the geometrical complexity of the fault. The amplitudes of the initial stress components as well as the friction parameters are chosen to match observed macroscopic earthquake characteristics, such as rupture duration and seismic moment, but are laterally homogeneous across the fault system. Stresses are consistent with commonly used dynamic rupture parametrization for strike-slip faults assuming that the axis of the intermediate principal stress is vertical. The amplitude ranges between 0 and 500 MPa, similar to what is assumed in Aochi et al. (2003).

The parametrization of the linear slip-weakening frictional behaviour is constant across all fault segments. Friction parameters are based on laboratory experiments, assuming a static coefficient close to Byerlee’s coefficient and a high stress drop to facilitate rupture transfer. In accordance with Section 4.1, we nucleate by smoothly enforcing rupture in a patch of radius 3450 m at the inferred hypocentre location (Hauksson et al.1993).

All model parameters that are constant across the faults are summarized in Table 7. All depth-dependent initializations including the 1-D velocity structure are visualized in Fig. 10. No stochastic variations nor asperities at the fault are included. This choice allows us to focus on the first-order effects of fault geometries oriented in a regional background stress field on rupture dynamics and transfers, which have been found to be a dominant factor in previous work (Gabriel & Pelties 2014).

Figure 10.

Depth-dependent model parameters on- and off-fault. Principal initial stress components σ_xx, σ_yy, σ_zz, 1-D velocity structure for S and P waves (Graves & Pitarka 2010) and plastic cohesion adapted from Roten et al. (2015).

Table 7.

Constant on-fault model parameters.

	Parameter	Value
μs	Static coefficient of friction	0.55
μd	Dynamic coefficient of friction	0.15
Dc	Critical slip-weakening distance	0.6 m
c0	Frictional cohesion	2 MPa
h	Smallest fault element edge length	100 m
p	Polynomial degree	3
ν	Bulk friction	0.55
T|$v$|	Viscoplastic relaxation time	0.03s

Off-fault plasticity requires a cohesion model of the host rock. Cohesion is material dependent and altered by damage induced by previous earthquakes. Since these conditions are not well constrained, a wide range of damage levels and rock types may be included in earthquake simulations (Ma & Andrews 2010; Roten et al.2014, 2015).

In the Landers region, we may expect rather high cohesion as the main surficial rock type is granodiorite (Dibblee 1967). As a consequence, we base our cohesion model on values for undamaged, granite-type rock (Roten et al.2015). Additionally, cohesion is varied in depth aligned with the major zones of our velocity model (Fig. 10).

The cohesive zone width (Λ) varies greatly across the complex fault segments in this model. Therefore, we consider in the following the minimum cohesive zone width (Λ_min) instead of its average. We find Λ_min = 300 m in the elastic Landers scenario compared to 350 m in the corresponding simulation with plasticity. According to Section 4, Λ_min needs to be resolved by at least 1.5 elements for polynomial order p = 3 in the elastic case.

The mesh is discretized by tetrahedral elements constrained by an edge length of 100 m across the fault. We note that on-fault resolution is increased by a factor of 2 in comparison to the mesh used in Heinecke et al. (2014), while the polynomial degree is decreased from 5 to 3. This discretization results in a sufficiently high resolution of the cohesive zone width (three elements per Λ_min in the elastic case). To improve computational efficiency, the mesh features coarsening away from the fault by a gradation rate of 1.06 (Simmetrix Inc. 2017), which reduces the total mesh size from 191 million elements to 22 million elements.

A fully elastic simulation of a simulated time of 40 s requires 3 hr and 35 min on 240 MPI nodes, each using 28 OpenMP threads, on phase 2 of the supercomputer SuperMUC. The simulation’s output is written in an asynchronous manner, reserving additional 16 nodes for output only (Rettenberger et al.2016). 3-D wavefield and plastic strain are written every second, 2-D fault quantities are written every 0.1 s. Additionally, time-series are recorded at 60 on-fault and 380 off-fault measuring points. The simulations with plasticity use the NB approach, which is beneficial for small fault discretizations (here, h = 100 m) and high polynomial degrees (here, p = 3), as discussed in Section 4.4.1. A full production run accounting for off-fault plasticity increases the computational cost by only 7.2 per cent compared to an elastic simulation. The additional cost agrees well with the estimates given in Section 2.3.3.

5.2 Geometrically complex rupture dynamics with off-fault plasticity

In the following, we discuss the effects of off-fault plasticity on dynamic rupture across a curved, branched and segmented fault system. We compare the key aspects of a fully elastic Landers earthquake scenario to one with plastic yielding. It is important to note that off-fault plasticity does not change the closeness to failure across fault segments, as equivalent initial conditions are assigned in both models. Thus, all observed differences in rupture dynamics are caused by the adjustment of stresses around the fault.

Macroscopic source characteristics of both scenarios are compared in Table 8. The maximum peak slip rate is overall reduced due to off-fault plastic yielding, whereas maximum slip increases by almost 13 per cent. Remarkably, the moment magnitude is almost identical in both scenarios.

Table 8.

Comparison of on-fault macroscopic characteristics of the Landers earthquake scenarios with purely elastic and plastic material response.

Parameter	Elastic	Plastic	Difference
Max. peak slip rate (m s−1)	11.46	9.38	−22.17 per cent
Max. slip (m)	5.54	6.36	+12.9  per cent
Moment magnitude Mw	7.396	7.382	−0.19  per cent

Fig. 11 visualizes the accumulated plastic strain surrounding the fault system. Plastic strain accumulates in the vicinity of geometrical complexity, as fault branching points and upon abrupt changes of fault orientation, that is, when the locally acting fault stresses change direction. Furthermore, plasticity is triggered at fault endings where rupture is not stopping smoothly.

Figure 11.

Map view of the accumulated plastic strain (see eq. 9) in the Landers fault system scenario (fault traces in white). The inset shows a cross-section at a fault branching point marked by the black line in the main figure.

In the following, we discuss differences of the elastic and plasticity scenario in terms of rupture cascading, slip rates and slip distribution. Rupture propagation in terms of slip rate across the full fault system is illustrated in Fig. 12 for both scenarios. Rupture path and rupture timing is very similar at the first segment, the Johnson Valley fault (JVF), up to 6 s of propagation time.

Figure 12.

Temporal evolution of slip rate across the Landers fault segments (JVF, LKF, HVF, the CF segment between the EF and CRF). The purely elastic scenario (left) and the scenario with plastic material response (right) differ in rupture transfer location and mechanisms (dynamic triggering versus branching) and location, marked by red rectangles.

In both simulations, rupture at the Homestead Valley fault (HVF) is triggered dynamically by waves emitted by the failure of the preceding fault segment. However, fault zone plasticity affects not only the timing of rupture transfer, but also its location: In the purely elastic model, we observe rupture transferring from the JVF to HVF after 6.8 s of propagation time. HVF is dynamically triggered (rupture ‘jumps’) very close to the branching point with the Landers Kickapoo fault (LKF; highlighted by a red rectangle in Fig. 12 at 7.1 s snapshot, left). In the simulation with plasticity, this rupture jump is delayed to 7.4 s and shifted to a more distant location on HVF (highlighted by a red rectangle in Fig. 12 at 7.9 s, right). After 8 s of simulation time, we observe the evolution of multiple rupture fronts and back propagation of rupture at HVF in both simulations.

At around 11 s, rupture transfers from the HVF to the Emerson Fault (EF) by different mechanisms: In the elastic simulation rupture is branching directly onto the connecting fault (CF) segment between the HVF and the EF before jumping to the EF (highlighted by a red rectangle in Fig. 12 at 11.3 s, left). However, the accumulation of plastic strain in the vicinity of geometrical complexity, such as fault branching points, suppresses direct branching to CF. EF is instead dynamically triggered when we account for off-fault plastic yielding. At the EF, slip rate is considerably reduced due to plastic yielding (see Fig. 11) and rupture nearly arrests after 12.1 s.

The last fault segment, the Camp Rock Fault (CRF), is ruptured several seconds later compared to the elastic case. As a result of the plastic energy dissipation around the first segments (see JVF, LKF and HVF in Fig. 11), the last two segments (EF and CRF) break only partially, specifically in regions of favourable orientation with respect to the background stress field.

The delay of rupture arrival as well as the damping of peak slip rate in the plasticity scenario become apparent in the time-series of slip rate recorded at on-fault locations at the JVF (r1), the HVF (r2) and the EF (r3) in Fig. 13. Peak slip rates are reduced by 32–48  per cent, rupture arrival is delayed by up to 1 s.

Figure 13.

Time-series of along-strike slip rate recorded at three points of the fault system (marked as r1, r2 and r3 in Fig. 14). Grey denotes purely elastic, orange plastic material response. Similar line styles correspond to equivalent on-fault location.

Overall peak slip rate across the fault system for both simulations are compared in Fig. 14. We observe a general reduction of peak slip rates for plasticity, specifically across the JVF, HVF and EF. The maximum peak slip rate across the entire fault is reduced by 22.17  per cent in the plasticity case (Table 8). In pointwise comparison peak slip rate differs up to 47  per cent.

Figure 14.

Peak slip rate of the Landers earthquake scenarios with purely elastic (top) and plastic (bottom) material response. On-fault recording locations of Fig. 13 are marked by triangles.

A reduction of peak slip rate is consistent with the previous work in 2-D and 3-D (e.g. Andrews 2005; Roten et al.2015). In contrast to a planar fault setup, in which peak slip rate is mainly reduced near the surface (Roten et al.2015), a reduction along the entire fault is possible when accounting for geometrical complexity of the fault. In particular, we observe direct correlation of strong peak slip rate reduction and plastic strain accumulation (Figs 14 and 11). Also, peak slip rate is drastically reduced at localized patches at the HVF and the EM. At such locations, the change of local on-fault stresses due to the change of geometrical orientation alters rupture path and speed that triggers plasticity.

The final slip distribution of both simulations is visualized in Fig. 15. Even though the seismic moment is very similar (Table 8), we observe higher, but more localized maximum slip due to plastic yielding. In both scenarios, maximum slip is occurring at the HVF, where we find 6.36 m in the simulation with plasticity and 5.54 m in the purely elastic simulation. Slip across the JVF and the HVF is overall higher with plasticity, even though peak slip rates and rupture speed are reduced. These findings are consistent with the recent work on planar step-over faults (Nevitt & Pollard 2017) in which a elastoplastic continuum increased the maximum slip on the single fault segments in static FE simulations.

Figure 15.

Accumulated slip of the Landers earthquake scenarios with purely elastic (top) and plastic (bottom) material response.

We note that the observed reduction of slip near the surface is not related to the complex fault geometry but rather to the decreasing bulk cohesion in combination with decreasing stresses towards the free surface: The interaction of rupture with the free surface triggers plastic yielding which, in turn, results in a higher plastic energy absorption with decreasing depth (Roten et al.2015).

In summary, plasticity considerably affects rupture dynamics at complex fault geometries: In addition to a reduction of rupture speed and peak slip rate as observed on planar faults, rupture path, transfer mechanisms and slip distribution are altered. The effect of plasticity is highest wherever rupture is encountering geometrical fault complexity such as bending, branching or abrupt ending. At these locations, initial background stresses are highly altered by locally different fault orientations or rupture is stopped abruptly. These findings compare well to simulations accounting for smaller scale fault roughness (e.g. Dunham et al.2011b; Ulrich & Gabriel 2017), which find plastic strain accumulation even triggered by slight deviations of planarity. Most crucially, we find that the location and timing of rupture transfers may be distinctively affected in large-scale segmented fault systems.

6 DISCUSSION

6.1 Damage rheologies for dynamic rupture problems

This study describes implementation, verification and application of a widely used proxy for plastic material deformation around rupturing faults for a high-order ADER-DG method. The main advantage of the plasticity approximation is to preserve the numerical method solving the linear elastic equations, while reducing potentially unrealistic high stresses around the rupture tip, often observed in dynamic rupture simulations (Andrews 2005; Dunham et al.2011a). Here, we only update the stress state in a separate plasticity corrector step succeeding the elastic algorithm.

Within this approach, however, the material constitutive behaviour encapsulated by the elastic wave equation remains unchanged. Therefore, no direct effect of the damaged material response on seismic wave propagation can be modelled. In contrast, numerical models have been developed to consider the change of elastic moduli surrounding dynamic rupture (Xu et al.2015; Lyakhovsky et al.2016; Okubo et al.2017; Thomas et al.2017). These studies report additional effects, such as slip rate oscillations and as a consequence an increase of near-fault high-frequency radiation. However, continuum damage rheologies are currently restricted to 2-D due to numerical challenges in solving the additional nonlinear equations.

Future work could compare the effects of damage rheology versus plastic yielding, focusing on dynamic triggering and branching in complex fault systems. Such studies could shed light on how rupture transfer may interact with the reduced wave velocities in off-fault damage zones.

6.2 Convergence studies for dynamic rupture problems

Due to the lack of analytical solutions for dynamic rupture problems, we verify the off-fault plasticity implementations presented here in two complementary manners: First by comparing SeisSol’s solution of a well-defined test problem to other state-of-the-art numerical methods (Section 3) and second by conducting h- and p-refinement studies for pointwise on-fault parameters with respect to a given reference solution (Section 4). The latter approach additionally constrains a minimum required resolution on the fault for accurately resolving source dynamics. We point out that the required resolution differs for different numerical approaches: for example, Day et al. (2005) conclude that for a Boundary Integral Method three points and for a second-order FD method five points per cohesive zone width are sufficient to accurately resolve on-fault dynamics in elastic simulations with homogeneous initial stress conditions. However, in setups with heterogeneous on- and off-fault stresses, initial conditions depend on the underlying discretization. The analysis in Section 4.3 indicates that peak slip rate is highly sensitive to the on-fault resolution of initial conditions, while other parameters such as rupture arrival time and peak slip rate time seem more robust. We expect similar effects for other numerical methods in setups with heterogeneous initial conditions.

Since no analytical solution exists for application-based simulations, the reported convergence rates in Section 4 highly depend on the choice of the reference solution. A suitable reference solution is commonly computed at a high mesh resolution and a high polynomial degree of basis functions, if applicable. However, if the reference solution still changes with further h and prefinement, the reported RMS errors may considerably underestimate the true error with respect to a solution with stable on-fault dynamics. As a consequence, the convergence rates (slope of the RMS regression lines with respect to the reference solution) might not be representative for the numerical scheme. The convergence analysis of Pelties et al. (2012) uses a reference solution with 354 m on-fault discretization, corresponding to a minimum cohesive zone width resolution of 0.91 elements. However, in this study we choose a distinctively higher resolution of the minimum cohesive zone width to ensure stable on-fault dynamics of the reference solution under heterogeneous initial conditions (2.28 mesh elements without and 4.58 mesh elements with plasticity).

6.3 Limitation of convergence rates for on-fault dynamic rupture

As discussed in Sections 4.3.2 and 6.2 the convergence rates reported here are representative of analytical convergence of the numerical scheme only if the pre-defined stand-in reference solution is sufficiently close to an unknown exact solution. Specifically, asymptotic convergence behaviour is only expected for sufficiently small discretizations h. Supplemented by the cohesive zone size analysis Appendix  C, we assume in the following discussion that the chosen high-resolution reference solution fulfils this requirement. We then discuss the origin of non-spectral convergence rates by interpreting rupture dynamics simulations in the framework of additional existing theoretical considerations. To our knowledge, a mathematical theoretical framework explaining the observed method-independent restriction of dynamic rupture problems to low-order convergence has not been published.

In general, the theoretical convergence rate of a numerical scheme is bounded by the degree of piecewise polynomials representing the projection space of the exact solution of the problem (Hesthaven & Warburton 2010, Chapter 4.5, p. 87). In problems including high discontinuities or steep gradients, this degree might be distinctively smaller than the polynomial degree p of its numerical representation. In this case, an increase of p may not further increase the convergence rate. However, for realistic dynamic rupture problems, it is not possible to formulate an exact solution to establish a theoretical accuracy barrier for on-fault quantities. In fact, the origin of an eventual underlying general ‘shock-like’ character of the dynamic rupture problem remains a topic of discussion.

Earlier studies suggest that the representation of the friction law causes an inherent non-smoothness of dynamic rupture problems (Kaneko et al.2008; Rojas et al.2008). The often used linear slip-weakening friction law, for example, poses a piecewise linear function. Its derivative exhibits singularities when slip equals zero and D_c (critical slip distance), yielding a discontinuity of slip acceleration despite smooth slip velocity and shear stress across the process zone (Ida 1973). Previous attempts of analysing smoother friction law formulations, such as the rate-and-state friction law, indicate slightly increased convergence rates (Kaneko et al.2008; Rojas et al.2009; Nerger et al.2014). However, for reaching spectral convergence properties fault resolutions several orders of magnitude higher than given by the cohesive zone size may be required.

Additionally, non-smoothness may be induced by the standard modification preventing tensile normal stresses |$\sigma _n^{\mathrm{ mod}} = \min (0,\sigma _n)$| to avoid fault opening. Most dynamic rupture solvers, including SeisSol, allow slip to only occur in shearing directions (mode II and mode III fracture propagation). The formation of tensile cracks away from the prescribed fault surface has been accounted for (Dalguer et al.2003; Dalguer & Day 2009), but is numerically challenging to incorporate and awaiting convergence analysis. Plasticity is expected to prevent fault opening for a wide range of realistic material parameters (e.g. Dunham et al.2011b; DeDontney et al.2012), thus regularizing the solution. However, our convergence analysis with plasticity (Section 4.4) does not result in higher convergence rates than observed in the elastic study. The Return Mapping Algorithm uses a low-order operator splitting approach (e.g. Simo & Hughes 2000) to solve for the elastic and plastic strain successively which may additionally limit the expected convergence to second order. We also note that in the current implementation, the plastic yield criterion is checked after a full time increment Δt, while the dynamic rupture boundary condition is evaluated at each time integration point (Pelties et al.2012).

Low-order convergence of dynamic rupture problems can be furthermore interpreted in the framework of Godunov’s theorem (Godunov 1959), predicting non-monotonicity in the vicinity of steep gradients of the solution for high-order linear solvers. SeisSol uses a linear upwind flux, the Godunov flux, which numerically dissipates oscillations, whether of numerical or physical origin. However, spectral convergence properties might be reduced to low-order accuracy due to the manifestation of the well-known Gibbs phenomena in the vicinity of strong discontinuities (e.g. Hesthaven & Warburton 2010, Chapter 5.6, p. 136). To overcome the limiting characteristics of linear high-order schemes at steep gradients, flux limiters, such as WENO (Weighted Essentially Non-Oscillatory) schemes (Dumbser et al.2007; Krivodonova 2007; Hesthaven & Warburton 2010), are promising approaches to avoid oscillations while preserving the overall high-order accuracy of the underlying method.

In contrast, high-order convergence can be established using the MMS (e.g. Roache 2002). Within this approach, the dynamic rupture problem is transformed into a related problem for which an analytical solution is known a priori (e.g. Kozdon et al.2012). Major simplifications of the boundary and interface conditions and the definition of additional source terms are required. Such solutions have been used to verify stability and convergence of a numerical method (Duru & Dunham 2016; Erickson et al.2017) without relying on numerical reference solutions. Nevertheless, the currently established manufactured solutions are too simplified to shed light on the convergence behaviour of typical dynamic rupture applications, in particular, of ones including off-fault plasticity.

The reported low-order convergence rates of peak slip rate, peak slip rate time and final slip (in Section 4) are consistent with previous dynamic rupture convergence studies (Day et al.2005; Rojas et al.2008; Pelties et al.2012). As a notable exception, Day et al. (2005) and Pelties et al. (2012) report convergence rates higher than 2 for rupture arrival time. However, the therein reported absolute RMS errors of rupture arrival time are considerably higher than in this study: for example, Pelties et al. (2012) report an RMS error of 2.79  per cent for p = 3 and h = 1061 m, while the RMS error of our solution is only 0.65 per cent for the same polynomial order and fault resolution. In previous studies, increasing the fault resolution highly improves the rupture arrival time RMS error and therefore increases the slope of the last square fit, that is, the convergence rates. We conclude that in contrast to the overstressed nucleation patch used in Day et al. (2005) and Pelties et al. (2012), the smoothly initiated rupture here seems to highly decrease rupture arrival time errors, specifically for lower resolutions.

We emphasize that using a high-order scheme results in low dispersion properties (Dumbser & Käser 2006), independently of the achieved order of accuracy. Furthermore, the absolute RMS errors decrease smoothly, even if high-order convergence is not reached.

A potential advantage in terms of computational efficiency when combining small elements with p = 1 at the fault with high-order discretization of the wave propagation part remains for future work: computational cost per element increases roughly with |$\mathcal {O}(p^6)$|⁠. However, hardware components are more efficiently employed at higher GFLOP s⁻¹ counts leading to a nonlinear relationship of p and time-to-solution.

6.4 Cohesive zone width estimation

The analysis of the cohesive zone width Λ in Appendix  C reveals a high spatiotemporal variability of Λ across the fault. The observed decrease of Λ with propagation distance for constant initial stress, as well as larger cohesive zone widths in corresponding plasticity simulations, is consistent with previous studies (Day et al.2005; Duan & Day 2008). Interestingly, our results indicate that an increase of frictional cohesion and decreasing initial stresses lead to an increase of Λ at shallow depths in both elastic and plastically yielding simulations. The interplay of the cohesive zone width with free surface effects remain to be investigated in detail. In general, our results imply that cohesive zone width estimates based on 2-D linear fracture mechanics (e.g. Andrews 1976; Day et al.2005) needs to be treated with care for setups including additional complexities such as heterogeneous stresses, frictional cohesion, 3-D geometries and plasticity.

6.5 Realistic parametrization of the Landers earthquake scenario

The demonstrator dynamic rupture scenario based on the 1992 Landers earthquake reproduces fairly well some of the main observed source characteristics: within the observed total rupture duration of around 24 s the rupturing fault segments produce a magnitude M_w = 7.3 earthquake. Not well constrained initial parameters are prescribed as simple as possible.

We observe that the orientation of the CRF hinders rupture propagation when accounting for plasticity due to high plastic energy dissipation across the previous fault segments. However, observations indicate that rupture reached the surface of most of the northern part (e.g. Wald & Heaton 1994; Hernandez et al.1999; Milliner et al.2015). Including additional fault weakening mechanisms, such as velocity weakening friction laws, thermal pressurization (Noda et al.2009) and thermal decomposition (Platt et al.2015) could help to sustain rupture at this fault segment. An alternative strategy is prescribing a different orientation of the regional stress field or different friction parameters across the CRF (e.g. Aochi & Fukuyama 2002).

The depth-dependent cohesion model is based on a 1-D approximation of granite-type rock (Roten et al.2015). Hoek–Brown parameters are used to calculate equivalent cohesion for a variety of rock types classified by their rock strength and the level of damage. Nevertheless, cohesion remains ill-constrained, as it strongly depends on the damage level of the host rock, which is potentially varying along and away from the fault. We observe that the choice of cohesion highly influences the fault dynamics and the evolution of plastic strain around the fault. Similar effects are reported for ground motions (e.g. Ma & Andrews 2010; Roten et al.2014). Damaged zones near the surface or in the close vicinity of faults may imply very small values of cohesion values, additionally influencing dynamic rupture transfer mechanisms. As lower cohesion values imply a lower yield surface for plasticity, off-fault plastic yielding would be triggered more easily, impeding sustained rupture propagation and transfer to subsequent segments.

7 CONCLUSIONS

We detail two implementation schemes for off-fault plasticity for the 3-D dynamic rupture software package SeisSol. Both approaches incorporate plastic material response at subelement resolution; namely at 3-D QP or at NB coefficients. The presented approaches are applicable to any modal DG method. The algorithms are formulated as matrix–matrix products which can be efficiently calculated, for example, via an offline code generator (Uphoff & Bader 2016). The NB approach is, in general, cheaper by up to a factor of 6, as it uses considerably less points than the QP approach, but allows lower subelemental resolution. Plasticity is numerically regularized by a mesh-independent viscoplastic relaxation time T_|$v$|.

We verify our approaches first by comparison to SCEC community benchmark solutions combining off-fault plastic yielding with realistic fault geometries and heterogeneous initial conditions. Remarkably, on-fault agreement improves distinctively when incorporating plastic material response. Off-fault comparison verifies high-quality solutions calculated on large elements while using high polynomial degrees also for plastically deforming materials.

To reflect the current state-of-the-art of realistic large-scale dynamic earthquake simulations, we provide for the first time 3-D dynamic rupture h- and p refinement studies including depth-dependent initial stresses and off-fault plasticity. The observed general smooth convergence of on-fault dynamics towards the reference solution illustrates the robustness of the numerical method. We establish hands-on guidelines for the minimum required on-fault resolution to accurately resolve source dynamics. For dynamic rupture models with plastic yielding a high resolution at the fault may be less crucial than commonly assumed: In comparison to an elastic simulation, plasticity increases the cohesive zone width due to decreasing rupture speed and increasing plastic yielding close to the free surface. Plasticity additionally regularizes the peak slip rate, which is very sensitive to the cohesive zone width resolution in purely elastic simulations. A too coarse fault resolution seems to highly impact the peak slip rate, which directly affects ground motion assessment. However, the required resolution for time-dependent quantities, such as rupture arrival time and peak slip rate time, increases in setups with plasticity. In summary, the overall minimum required resolution of the cohesive zone width is 34–50 per cent lower with plastic yielding, mainly due to the lower resolution required for the sensitive peak slip rate.

In terms of accuracy, the two presented implementations for plasticity give almost identical results for a cohesive zone resolution of more than 2.04–2.75 mesh elements (10.2–12.24 subelemental points) for p > 2. In these cases, it is beneficial to use the NB approach due to its higher computational efficiency. For polynomial degree p = 2 and lower resolution of the cohesive zone, the QP approach yields higher accuracy since the projection to the polynomial space is more accurate with respect to the L² norm than in the NB approach. On-fault convergence rates with respect to the reference solution do not depict spectral convergence in both cases. We critically discuss a possible theoretical limitation of dynamic rupture to low-order convergence rates.

We demonstrate the influence of non-elastic material behaviour on rupture transfers across a geometrical complex fault system. The earthquake scenario based on the 1992 Landers event reveals that plasticity considerably affects rupture dynamics across complex fault geometries. Plastic strain accumulates at deviations of the rupture path from planarity, for example, at changes of fault strike orientation, branching or segment endings. In particular, we observe direct correlation of strong peak slip rate reduction and plastic strain accumulation. As a consequence, the spatiotemporal distribution of rupture transfers is distinctively altered by off-fault plastic energy dissipation. Off-fault plasticity reduces the peak slip rate by up to 50  per cent and delays rupture arrival across the entire fault and to a larger extent than reported for scenarios on planar faults (Roten et al.2015). In the simulation with plasticity, slip is found to be locally higher, but accumulated across a smaller area. As a result, moment magnitudes are comparable with and without plasticity, even though the rupture path differs dynamically. We find that the cohesive zone width varies considerably across the fault system, implying that a minimum inherent length scale may be sought to be resolved instead of an average.

Our study emphasizes the importance of combining fully 3-D dynamic modelling, complex fault geometries and off-fault plastic yielding to realistically capture dynamic rupture transfers between fault segments and further the understanding of the activation of fault branches and the potential for dynamic triggering of adjacent fault segments.

ACKNOWLEDGEMENTS

The work presented in this paper was supported by the German Research Foundation (DFG; project no. KA 2281/4-1, AOBJ 584936/TG-92); by the Bavarian Competence Network for Technical and Scientific High Performance Computing (KONWIHR), project GeoPF (Geophysics for PetaFlop Computing); by the Volkswagen Foundation (project ASCETE – Advanced Simulation of Coupled Earthquake-Tsunami Events, grant no. 88479); by the European Union’s Horizon 2020 research and innovation program under grant agreement no. 671698 and by King Abdullah University of Sciene and Technology (KAUST) in Thuwal, Saudi Arabia, under grant ORS-2016-CRG5-3027-04 and OSR-CRG2017-3389. Computing resources were provided by the Leibniz Supercomputing Centre (LRZ, projects no. h019z, pr63qo and pr45fi on SuperMUC). The authors thank Dave A. May and Kenneth C. Duru for fruitful discussions. Furthermore, we thank LMU Geophysics' system administrator team for their support (Oeser, 2009). The manuscript benefited from valuable reviews and comments by R. Ando and an anonymous reviewer.

APPENDIX A: DERIVATION OF A CLOSED-FORM UPDATE FORMULA FOR THE RETURN MAPPING ALGORITHM

After solving the dynamic rupture problem for purely elastic material behaviour, we update the elastic trial stress state in case of plastic yielding (F(σ) ≥ 0). We calculate the additional change in stress due to plasticity, assuming the trial stress state |$\sigma _{ij}^{\text{trial}}$| at the beginning of the time step. The goal is to derive the updated stress state |$\sigma _{ij}^{n+1}$| solving the following equation:

\begin{equation*} \dot{{\sigma }}_{ij} = -2\mu \dot{\epsilon }_{ij}^p = -\frac{1}{T_v} (\sigma _{ij}-P_{ij}(\boldsymbol {\sigma })), \end{equation*}

(A1)

where |$\dot{\epsilon _{ij}^p} \ne 0$| is defined by the flow rule given in eq. (10).

We first multiply this equation by the integrating factor exp ( − (t_{n + 1} − t)/T_|$v$|) and then perform integration over the time interval of Δt = t_{n + 1} − t_n to get

\begin{equation*} \begin{split} \sigma _{ij}^{n+1} = & \exp (-\Delta t/T_v) \ \sigma _{ij}^{\text{trial}} \\ & + \int _{t_n}^{t_{n+1}} \exp (-(t_{n+1}-t)/T_v)\frac{1}{T_v} P_{ij}(\boldsymbol {\sigma }) \ \mathrm{ d}t. \end{split} \end{equation*}

(A2)

Note that we make use of the relation

\begin{equation*} \begin{split} \int _{t_n}^{t_{n+1}}(\dot{{\sigma }}_{ij} + \frac{1}{T_v} \sigma _{ij} )\ \exp (-(t_{n+1}-t)/T_v) \ \mathrm{ d}t &= \\ \int _{t_n}^{t_{n+1}} \frac{\partial }{\partial t} [ \sigma _{ij} \exp (-(t_{n+1}-t)/T_v)] \ \mathrm{ d}t &= \\ \sigma _{ij}^{n+1} + \exp (-\Delta t/T_v) \ \sigma _{ij}^{\text{trial}} & \end{split} \end{equation*}

(A3)

to derive eq. (A2). Assuming |$P_{ij}(\boldsymbol {\sigma })$| to be constantly equal to |$P_{ij}(\boldsymbol {\sigma ^{\text{trial}}})$| over the time interval [t_n, t_{n + 1}], the integral in eq. (A2) can be approximated as

\begin{equation*} \begin{split} \int _{t_n}^{t_{n+1}} \frac{1}{T_v} \exp (-(t_{n+1}-t)/T_v) \ \mathrm{ d}t \ P_{ij}(\boldsymbol {\sigma }^{\text{trial}}) = \\ (1-\exp (-\Delta t/T_v))\ P_{ij}(\boldsymbol {\sigma }^{\text{trial}}) &. \end{split} \end{equation*}

(A4)

Combining now eqs (A2) and (A4) and inserting P_ij as defined in eq. (11) results in

\begin{equation*} \begin{split} \sigma _{ij}^{n+1} \approx & \exp (-\Delta t/T_v) \ \sigma _{ij}^{\text{trial}} \\ & + (1-\exp (-\Delta t/T_v))\left(\frac{\tau _c}{\sqrt{I_2}} \ s_{ij}^{\text{trial}} + \sigma _m^{\text{trial}}\delta _{ij}\right) \\ = & ((1-\exp (-\Delta t/T_v))\frac{\tau _c}{\sqrt{I_2}} + \exp (-\Delta t/T_v) ) \ s_{ij}^{\text{trial}} \\ & + \sigma _{m}^{\text{trial}}\delta _{ij}, \end{split} \end{equation*}

(A5)

where we recover the adjustment factor f* applied in eq. (14).

APPENDIX B: VERIFICATION OF QUADRATURE POINTS VERSUS NODAL BASIS PLASTICITY IMPLEMENTATION

We compare the verification results of the two implementation approaches for off-fault plastic yielding presented in Section 2.3 for the benchmark setup TPV27 described in Section 3.1. The QP and the NB implementation are employed to simulate plastic yielding around a strike-slip fault under heterogeneous initial stress conditions using an identical computational mesh and parametrization.

Fig. B1 depicts exemplarily slip rate and shear stress on the fault, as well as horizontal velocity off-fault at the same measuring locations as shown in Figs 3 and 4 (which include the QP approach only). We compare both SeisSol solutions using h = 250 m discretization and polynomial degree p = 4 to the second-order FE method FaultMod using h = 50 m. The near-perfect agreement of the QP and NB approach for on-fault dynamics and off-fault wave propagation include matching of the peak amplitude of the waveform in Fig. B1(c), which varies from the FaultMod solution.

Figure B1.

Comparison between FaultMod with a fault discretization of 50 m and SeisSol QP (black) and SeisSol NB (orange) with fault discretizations of 250 m for the benchmark problem TPV27. (a) Along-strike slip rate and (b) shear stress at a fault location 10 km downdip and 10 km along-strike, (c) horizontal velocity at a receiver 20 km off-fault. SeisSol’s QP and NB approach to implement plastic yielding perform equally well.

Overall, both approaches yield very similar results and are in equally good agreement with FaultMod serving as method-independent comparison solution. A detailed quantitative analysis of the trade-off between computational cost and accuracy of both implementations is presented in Section 4.4.

APPENDIX C: COHESIVE ZONE WIDTH FOR THE CONVERGENCE ANALYSIS

The cohesive zone (also known as breakdown or process zone) is defined as the area behind the rupture front in which shear stress decreases from its static to its dynamic value (Day et al.2005). Across this portion of the fault, the slip rate and shear stress significantly change. Therefore, the width of the cohesive zone represents a dynamic rupture problem inherent length scale which is crucial to accurately resolve. In this appendix, we determine the time- and location-dependent cohesive zone width Λ for the 3-D dynamic rupture setup described in Section 4 for elastic and plastic material properties under depth-dependent initial stresses. All calculations are based on the reference (i.e. high-resolution) solutions employing h = 71 m at the fault and polynomial degree p = 5.

Fig. C1 visualizes contour lines of the rupture front (RF) arrival time and the first point in time at which shear stresses reach their dynamic value (DS) across the fault for the plastically yielding convergence setup. The distance between contours at the same point in time (same colour) represents the cohesive zone width. Below 15 km depth, the shear stress is tapered to stop rupture smoothly at the end of the fault. The rupture front does not reach the free surface due to plastic energy dissipation and increasing frictional cohesion in the shallow part of the fault.

Figure C1.

Contour plot of the rupture front (RF) arrival time (solid line) and the first point in time at which shear stresses (DS) reach their dynamic value (dashed line), plotted at each second of the plasticity reference solution (QP approach) for the convergence test setup defined in Section 4. The star indicates the hypocentre, the triangle the location shown in Fig. 9. The zoom shows the distribution of SeisSol’s integration points (QPs, grey) within one element for polynomial degree p = 3 for a h = 1061 m discretization at the fault. The black points denote the gridpoints of the mesh. The cohesive zone contours in the zoom are plotted at t = 5 s.

SeisSol’s underlying numerical scheme defines shear and normal stress at 2-D QPs located inside each tetrahedral element face which is linked to the fault, as discussed in Section 2.1.2. The grey dots in the zoom-in of Fig. C1 visualize the distribution of QPs exemplarily for the lowest resolution we consider in the convergence analysis. Note that the effective numerical discretization is distinctively higher than the mesh discretization, denoted by black dots. From Fig. C1 we see that the effective discretization at the fault can be approximated by the edge length of the mesh divided by p + 1 additional integration points.

Fig. C2 depicts the spatiotemporal evolution of the cohesive zone width for a fixed depth along-strike and for a fixed along-strike location along-dip (capturing depth-dependent initial conditions) for the elastic and the plasticity reference solutions. The figure is plotted analogous to fig. 3 in Day et al. (2005).

Figure C2.

Cohesive zone width in dependence of rupture propagation distance and time (cf. fig. 3 Day et al. (2005)). Solid lines denote the rupture front (RF) arrival and dashed lines the dynamic shear stress (DS) arrival. Upper panel: along dip evolution at 0 m along-strike, lower panel: along-strike evolution at −10 500 m depth. In grey the cohesive zone width of the purely elastic reference (high-resolution) solution is given, in orange the equivalent with off-fault plasticity.

In the elastic simulation (Fig. C2, grey), the cohesive zone decreases along strike: from 706 m at 2 s simulation time to 214 m once rupture reaches the end of the fault. This reflects the decreasing width of the cohesive zone with increasing rupture velocity, consistent with findings in Day et al. (2005). In contrast to the decreasing cohesive zone width for homogeneous initial stresses without frictional cohesion (Day et al.2005), the cohesive zone width in our setup increases from 850 m to a maximum value of 2204 m at a shallow depth of −1800 m in the up-dip direction . The increasing frictional cohesion towards the free surface slows down rupture and results in a larger cohesive zone.

In the plastically yielding case (Fig. C2, orange), the cohesive zone width Λ decreases with propagation distance along-strike: from around 690 m just after nucleation to a minimum of 460 m. Note that the minimal width reached is twice as large as in the purely elastic simulation. For up-dip rupture propagation the width increases from 901 m after nucleation to a maximum value of 1555 m (at a depth of −3000 m). This reflects a less-pronounced increase than in the elastic setup, while only in the latter rupture reaches the free surface.

The cohesive zone width Λ is overall smaller in purely elastic simulations than when accounting for off-fault plasticity. A natural exception poses the region close to the free surface, which does not rupture in the plastic simulations (and is excluded in the evaluation of convergence characteristics).

To determine the minimum required discretization for accurate (i.e. converged) on-fault solutions, we calculate the cohesive zone width at each on-fault location used in the convergence analysis of Section 4. To this end, the time difference (DS−RF) is multiplied by the rupture speed measured at each fault location under consideration.

Across all fault receivers, the minimum cohesive zone width is evaluated as |$\Lambda _{\mathrm{ min}}^\mathrm{ e}$| = 162 m in the elastic and |$\Lambda _{\mathrm{ min}}^\mathrm{ p}$| = 325 m in the simulation with plasticity. The median cohesive zone width is evaluated as |$\bar{\Lambda }^\mathrm{ e}$| = 583 m and |$\bar{\Lambda }^\mathrm{ p}$| = 722 m, respectively.