Verification of an ADER-DG method for complex dynamic rupture problems

We present thorough benchmarking of an arbitrary high-order derivative Discontinuous Galerkin (ADER-DG) method on unstructured meshes for advanced earthquake dynamic rupture problems. We validate the method in comparison to well-established numerical methods in a series of verification exercises, including dipping and branching fault geometries, heterogeneous initial conditions, bimaterial cases and several rate-and-state friction constitutive laws. We show that the 5 combination of meshing flexibility and high-order accuracy of the ADER-DG method makes it a competitive tool to study earthquake dynamics in complicated setups.


Introduction
The combined simulation of dynamic fault rupture and seismic wave propagation is a useful tool to gain insight into the poorly constrained processes of earthquake faulting.Dynamic rupture modeling aims to understand the underlying physics governing earthquakes and may be incorporated in physics-based seismic hazard assessment and strong motion prediction in preparation of future, possibly devastating, events (Ely et al., 2010;Harris et al., 2011;Roten et al., 2011).
Recent advances in dynamic rupture simulations furthered our understanding of the earthquake cycle in the Parkfield region (Barbot et al., 2012), the influence of low velocity fault zones (Huang and Ampuero, 2011) and off-fault plasticity (Templeton and Rice, 2008;Ma and Andrews, 2010;Gabriel et al., 2013) on the dynamics of the fracture process, mechanisms to generate pulse-like and supershear earthquakes and their consequences (Dunham, 2007;Daub et al., 2010;Gabriel et al., 2012) and the interaction of fault zone branches (Oglesby et al., 2003;Bhat et al., 2007;DeDontney et al., 2011).Dynamically consistent predictions of strong ground motion excitation in earthquake scenarios have been pushing computational performance to the petaflop scale (Cui et al., 2010;Zhou et al., 2013).
Nevertheless, despite the current achievements, the numerical simulation of the rupture process and its implementation into elastodynamics solvers poses challenges.For instance, the natural discontinuity of physical variables across the fault interface has to be represented accurately and accounted for by the computational mesh.Furthermore, the dynamic rupture implementation in many numerical methods, such as finite difference, finite element and spectral element methods, suffers from spurious high-frequency oscillations.These oscillations may counteract the non-linear wave and rupture interaction in dynamic earthquake simulations and potentially contaminate the solution over all space-time scales (Duan and Day, 2008).Thus, numerical regularization, artificial attenuation or smoothing is usually necessary to suppress high-frequency numerical noise with spatial scale near the resolution limit of the mesh.Such artificial damping mechanisms, as for example the deployment of a thin layer of Kelvin-Voigt viscous material surrounding the fault, e.g.Ampuero (2008); Day et al. (2005); Dalguer and Day (2007), are not completely satisfying.The physical solution is not necessarily insensitive to the precise parametrization of the added damping.The damping may interfere with the actual physics of interest, for example by slowing down the rupture propagation (Andrews, 2005) and smoothing out small scale features.The artificial damping may also reduce the time step length and thus increase the computational effort considerably.
Realistic earthquake scenario simulations would ideally cover a frequency range relevant for engineering applications (up to 20 Hz), and include geological models spanning hundreds of kilometers consisting of complicated fault geometries, topography, oceans and low velocity sedimentary basins.
Simultaneously, small-scale high resolution is required to resolve the frictional sliding in the socalled cohesive zone at the rupture tip with a sufficient number of computational nodes or elements, e.g. at scales down to meters in case of slip-weakening constitutive behavior (Day et al., 2005).
The constitutive laws describing the frictional sliding of faults originate in laboratory experiments carried out on much smaller scales than in natural setups (Dieterich, 1978;Ohnaka and Kuwahara, 1990;Di Toro et al., 2005).While it is still uncertain how to extrapolate fault constitutive properties from laboratory to natural scales, testing the large scale implications of laboratory friction poses an enormous computational challenge.Additionally, the propagation of small wavelengths over large distances accumulates numerical dispersion and diffusion errors and favors a high-order accurate discretization.Another challenge for dynamic rupture simulations are large uncertainties in physical initial conditions and fault constitutive parameters.The rupture evolution depends strongly on model parameters such as initial background stresses, nucleation procedure and frictional properties.
In order to avoid additional errors that interfere with the physical problem accurate numerical methods that produce reliable results are desirable.Furthermore, computational efficiency, parallelization and high scalability are crucial demands for numerical methods simulating realistic earthquake scenarios.
In this paper, we present a thorough verification of the software SeisSol (Käser and Dumbser, 2006;Dumbser and Käser, 2006), a high-order derivative Discontinuous Galerkin (ADER-DG) method on unstructured meshes, for advanced dynamic rupture problems (de la Puente et al., 2009;Pelties et al., 2012).In contrast to the well verified and validated simulation of seismic wave propagation, the verification process of spontaneous dynamic rupture simulations suffers from missing analytical reference solutions.Therefore, we verify the performance of the ADER-DG method in the benchmark suite established by the Southern California Earthquake Center (SCEC) (Harris et al., 2009(Harris et al., , 2011)).All simulation results presented here are available at http://scecdata.usc.edu/cvws/.To this end, we cover many aspects important in realistic faulting setups, including dipping and branching fault geometries, bimaterial cases, heterogeneous initial conditions and different formulations of rate-and-state friction laws.
Our results demonstrate the benefits of SeisSol for dynamic rupture and ground motion simulations in realistic settings.Importantly, we confirm the lack of systematic numerical artifacts in advanced faulting setups, as reported for the basic example of a planar fault embedded in a homogeneous full space (SCEC Test Case TPV3) by Pelties et al. (2012).
2 Numerical Method de la Puente et al. (2009); Pelties et al. (2012) presented a new, alternative numerical scheme for the simulation of earthquake faulting.The underlying solver for wave propagation is an ADER-DG method (Käser and Dumbser, 2006;Dumbser and Käser, 2006) with high order accuracy in space and time based on tetrahedral element discretization.Between any two elements information is exchanged via numerical fluxes.We briefly outline the algorithm to evaluate the friction law in appendix A and extend this approach to faulting at dissimilar material contacts.For further details on the concept, the reader is referred to de la Puente et al. (2009); Pelties et al. (2012).The software package SeisSol provides pre-and post-processing tools including interfaces to external mesh generators and a mesh partitioning concept based on METIS (Karypis and Kumar, 1998).
The use of a tetrahedral element discretization leads to rapid and automatized mesh generation that can be readily aligned to complex geometrical features.The possibility of mesh refinement is advantageous for dynamic rupture problems as the mesh resolution can be adapted to areas of interest, such as the fault plane or complex topography.Furthermore, the mesh size can be coarsened with increasing distance to the fault to reduce the computational cost.Note that no artificial reflections due to mesh coarsening are observed (de la Puente et al., 2009).
In contrast to the typically applied traction at split-node approach (Andrews, 1999;Day et al., 2005;Dalguer and Day, 2007), ADER-DG solves the frictional sliding via the inverse Riemann problem (Toro, 1999;LeVeque, 2002), in which the exact solution is modified to incorporate frictional boundary conditions.Solving the inverse Riemann problem inherits the favorable numerical properties from the exact Riemann solver or Godunov flux.The numerical properties of the ADER-DG algorithm are extremely sensitive to the choice of flux function.In particular, our implementation introduces a very selective numerical dissipation by employing an upwind flux, as discussed in Pelties et al. (2012).The dissipation increases with increasing frequencies and is stronger beyond a cut-off frequency that depends on the mesh size h and on the order of accuracy O.The cut-off frequency is expected to be inversely proportional to the travel time of S waves over a typical grid spacing ∼ h/O/c s .Higher frequency modes are subdued while the physically meaningful lower frequencies are minimally affected.This is advantageous for dynamic rupture simulations: spurious high-frequency oscillations are not generated in the first place and, thus, no additional damping procedures need to be applied.
For consistency, we compare our results with results of the well-established software FaultMod (Barall, 2009) throughout this paper.FaultMod is a finite element (FE) code implemented on hexahedral elements and designed specifically for constructing physics-based models of fault systems which involve complex 3D geometry and 3D variation of material properties.The implementation of fault friction is based on the traction-at-split-nodes method.Faults are represented using common and differential nodes, which have a non-diagonal mass matrix.Also, FaultMod uses an implicit time stepping scheme, where displacement, velocities, and acceleration are computed simultaneously.The method implements Newmark damping (Hughes, 2000) and an optional thin viscous layer surrounding the fault zone (Day et al., 2005;Dalguer and Day, 2007) to suppress unwanted high-frequency oscillations.For specific details about damping tuning strategies we refer to (Rojas et al., 2008) or according FaultMod publications (Barall, 2009).In all here presented benchmarks FaultMod applies a grid doubling technique to enable high resolution at the fault as well as at several grid layers around it succeeded by a wider grid spacing at a certain distance from the fault.In section 7 we additionally show a comparison of our results with other high-order dynamic rupture codes based on a variety of numerical methods.
3 Dip-slip fault with depth dependent background stress conditions Many studies of earthquake source physics have focused on purely strike-slip faulting.Nevertheless, many other fault types exhibit dipping components including subduction zone settings, in which most of the seismic energy was released over the last century (Pacheco and Sykes, 1992).The asymmetric geometry of dipping faults is of particular interest, since it results in asymmetric nearsource ground motion (Oglesby et al., 1998).Furthermore, in case of surface rupture, the arrival of the rupture front at the free surface strongly excites seismic waves (Madariaga, 2003) and the dynamic rupture evolution interacts with reflected waves from the free surface (Huang et al., 2013).
Reproducing these effects accurately in numerical simulations poses challenges in terms of mesh generation and numerical stability.
We model nucleation followed by spontaneous rupture on a 60 degree dipping normal fault reaching a free surface in a homogeneous half-space.Subshear (SCEC Test Case TPV10) and supershear (SCEC Test Case TPV11) conditions are parameterized by varying the value of the static coefficient of friction.Rupture is initiated by increasing the background dip-direction shear stress until a value larger than the static yield strength in a pre-defined nucleation patch is reached.For the entire fault plane, including the nucleation patch, both the initial fault normal stress and the initial along-dip shear-stress increase linearly with down-dip distance.In order to accurately sample the initial stress and material parameters, we assign individual values to each Gaussian integration point (GP) used in the friction solver (Pelties et al., 2012).Slip-weakening friction and elastic off-fault yielding are assumed.All simulation parameters are summarized in Table 1.
We compare the ADER-DG O5 solution on an edge length of 200 m to the finite element method FaultMod of Barall (2009) with an edge length of 100 m and O2.Our mesh was gradually coarsened up to maximum edge lengths of 5 km far away from the fault in order to concentrate numerical costs where needed.an excellent agreement with the FEM reference solution, despite the asymmetric unstructured mesh surrounding the fault discontinuity in the ADER-DG case.Nevertheless, the already damped finite element solution still shows signs of high-frequency oscillations in slip-rate amplitudes.
The development of a supershear rupture front in TPV11 is equally well captured, as shown in Fig.
2. The rupture time contour plot in Fig. 2 (a) captures the boost in rupture velocity after supershear transition.We also point out the normal stress variation (Fig. 2 (b)) caused by the interaction of the non-vertical fault with the Earth's free surface as investigated e.g. by Rudnicki and Wu (1995); Dalguer et al. (2001); Ma and Archuleta (2006); Andrews et al. (2007); Ma and Beroza (2008).
Interestingly, we see differences in the evolution of stresses after 10 s: the FEM solution reaches higher normal and along-dip shear stresses leading to a slight difference in slip rate.This might be due to fault and free surface interaction, differently handled by both methods.
Effects of strong ground shaking due to free surface interaction, as well as an asymmetry between foot wall and hanging wall can be observed in Fig. 3. Larger ground motions on the hanging wall than on the footwall are observed in natural earthquakes, e.g.Abrahamson and Somerville (1996), and can be related to normal stress variation inducing strength drop variations.The agreement between ADER-DG and FEM regarding the synthetic ground motions in the vicinity of the fault is near-perfect.

Heterogeneous background stress fields
Tectonic loading plays a fundamental role in earthquake source dynamics and controls the size of an earthquake.Stress heterogeneities could potentially influence nucleation and arrest of a rupture (Day, 1982;Boatwright and Quin, 1986;Oglesby and Day, 2002;Ampuero et al., 2006).An open question is how small scale fluctuations of prescribed initial stress fields modulate the rupture process.As stresses in the Earth's interior can not be measured directly, dynamic rupture simulations represent a proper tool to analyze the impact of stress variations.
The SCEC benchmark problems TPV16 and 17 focus on the modeling of dynamic rupture under heterogeneous initial stress conditions.Here, we consider only TPV17, as TPV16 is very similar.
The randomly generated heterogeneous initial stress conditions and the frictional parameters are provided as pre-defined input.The setup contains a planar strike-slip fault embedded in a linear elastic medium with wave speeds specified in Table 2. Linear slip-weakening governs the frictional sliding.The nucleation method combines high stress values with an approximated radius of 1 km and low D c values with an approximated radius of 4 km in the hypocentral region.Furthermore, the friction coefficient is lowered to its dynamic value following a time dependent function to create an  expanding circular region of forced rupture propagation off the hypocenter.
The effect of heterogeneous initial stresses on dynamic rupture propagation needs particular scrutiny in methods based on high-order discretizations, which work more efficiently on large elements (Käser et al., 2008).Even if the dispersion requirements are sufficiently addressed by large elements and a high-order approach, yet another issue is the correct sampling of the given initial stress and friction data on the fault.Pelties et al. (2010) attempted to define rules to respect material properties of a complex geological medium correctly, but we are not aware of a study addressing this issue for the initial stress conditions of a dynamic rupture simulation.As introduced in section 2 and appendix A, in the case of slip occurring on the fault we modify flux functions in accordance with the Coulomb failure criterion.The fault plane is thereby always located at the interface between two adjacent elements.These flux functions are integrated on (N +2) 2 GPs irregularly distributed across each triangular element face, where N is the polynomial degree of the basis functions.We assign initial values of friction coefficient, critical slip distance, stress or other parameters to each GP at the fault interface individually.
In this particular example we use a trilinear interpolation algorithm to map the given gridded data on the irregularly distributed GPs.In this way, the smallest possible scale of the numerical method is exploited without decreasing the element size and thus sub-element resolution is enabled.A typical resolution of element edge length h = 200 m and O4 is applied in order to test if our proposed sub-sampling scheme is sufficient to capture the initial stress values and frictional parameters.Note again that mesh coarsening with increasing distance to the fault is applied.rate.With the aid of Fig. 7, we can identify this signal as the passage of an interface wave (Dunham, 2005).Reflection at the free surface causes a healing front which passes the receiver after the rupture front and is followed immediately by a secondary slip rate peak.Note that the amplitude of the latter peak even exceeds the peak slip rate of the primary rupture front.At other positions on the fault, a healing effect caused by interface waves is not as pronounced.We argue that this effect is physically correct within the limits of the provided model as it is also resolved by other numerical methods at high resolution (see the data provided by the SCEC Code Comparison tool).We emphasize the remarkable accuracy of the ADER-DG solver allowing to observe such small scale features with a relatively large element edge length of h = 200 m.We conclude that our approach to sample a given background model by the Gaussian integration points leads to sufficient conformity with FaultMod results.

Fault branching
Fault branching has been observed in natural events, e.g.Schwartz et al. (2012).Incorporating the possibility of earthquake rupture propagating from one fault to another may lead to earthquake scenarios involving very different magnitudes and dynamics, which would in turn be of crucial importance for seismic hazard analysis.However, the role of fault branching in earthquake dynamics is poorly understood.Realistic simulations need reliable information on fault geometries on the one hand, on the other hand the theoretical analysis of fault branching mechanisms is limited to simplified scenarios.Defining mathematical predictions for dynamic rupture propagation at the junction point of main fault and branch turns out to depend on the specific boundary conditions assumed at the junction (DeDontney et al., 2011).
To verify the performance of the ADER-DG algorithm on a branching fault system we discuss here SCEC's Test Cases TPV14 and TPV15.The geometrical setup contains two vertical, planar strike-slip faults; a main fault and a branch fault intersecting at an angle of 30 degree.The faults   The elliptical shape of the interface wave is clearly visible, which consists of a healing phase followed by a high slip rate front.
reach the Earth's surface.The initial stress conditions determine the fault system as right-lateral in TPV14 and as left-lateral in TPV15.Below we focus on TPV15.
The line where main fault and branch fault intersect is referred to as junction point and defines along-strike distance 0 km.However, as specified in the benchmark description, the branch fault should not fully reach the junction point.A small gap with the length of one grid spacing between the two faults is prescribed.This way, a rupture nucleated on the main fault can propagate freely onto the right side of the main fault passing the junction point, but the rupture must jump the distance of one grid spacing to propagate onto the branch fault.Since our method is able to model both the gap between main and branch fault as well as the setup without gap, we address both cases and discuss the differences in the following.However, in the modal ADER-DG formulation there is no solution defined exactly at the junction point, as physical solutions are only computed within the entire element and at the element edges, with exception of the vertices.We note that the branching benchmark problems TPV24/25 without a gap at the junction were more recently formulated and are not considered here.
It is important to notice that we use the same mesh for the simulations with gap and without gap.
Technically, this is achieved by simply locking a part of the branch (gap) in distinction to allowing the full branch to rupture (no gap).This way, all differences in the results are due to the influence of the gap exclusively.We choose the recommended node spacing of 100 m as gap length, although the generally applied mesh resolution at the fault is h = 300 m.This results in a slightly higher mesh resolution of up to h = 100 m in the direct vicinity of the junction point.We illustrate this in a zoomed view of the junction point in a 2D setup in Fig. 8.The different colors represent different MPI partitions for parallel computing.Fault and MPI domains can interact in arbitrary ways, e.g.intersecting perpendicular or being aligned parallel and combinations of both.Again, mesh coarsening with distance to the fault is applied.The order of accuracy in space and time is for all simulations O4.We summarize the complete parameters in Table 3.In Fig. 10 we present time series of shear stresses and slip rates in strike and dip direction of a station on the branch in 9 km along-strike distance to the junction point and 7.5 km depth.In general, we observe a good agreement between the four solutions.Only small differences in the dip components (Fig. 10 (b) and (d)) are visible at ∼ 7.5 s.Since the dip component is one order of magnitude smaller than the strike component, respectively, these differences are considered negligible.
Larger discrepancies can be observed at receivers closer to the junction point, for example at the main fault in 2 km along-strike disctance to the junction point and in a depth of 7.5 km as shown in Fig. 11.Whereas the shear stress in dip direction (Fig. 11 (b)) and the normal stress (Fig. 11 (c)) are very similar, the shear stresses along-strike differ clearly in all four simulations.In this case, we can observe that the 100 m results of FaultMod tend to be closer to the ADER-DG solution with gap and FaultMod's 50 m results tend to be closer to the ADER-DG without gap solution.We argue that this is expected as the geometrical setups are more similar.In general, all receivers far away from the junction point including the off-fault receivers (not shown here) match very well and only receivers in the direct vicinity of the branching point show discrepancies that can be traced back to the gap between main and branching faults.We would like to point out that such small differences in the geometrical model might lead to a different rupture behavior if the setup is more sensitive to the geometry or the initial stress values than prescribed in this benchmark.In agreement with results by DeDontney et al. (2011), an insensitivity of simulations of realistic branching scenarios to gaps implemented at branching points cannot be guaranteed and we advise modelers to handle branching geometries with caution.

Bimaterial contrast
Natural faults often separate rocks with different rheologies, e.g.Thurber et al. (2006), which leads to normal stress variations during faulting in turn influencing dynamic rupture propagation (Harris and Day, 1997;Ampuero and Ben-Zion;Brietzke et al., 2009).Under certain circumstances such stress perturbations can generate additional propagation modes of rupture, for example selfsustaining pulses, and are therefore of particular interest.The instantaneous response of shear stress to normal stress changes at bimaterial interfaces leads to an ill-posed problem for a wide range of elastic material contrasts causing an instability which can excite waves of all wavelengths (Adams, 1995).Furthermore, convergence through grid size reduction cannot be achieved in ill-posed setups and regularization needs to be applied.Thus, we address regularization concerns in appendix B.
Here, we concentrate on the well-posed bimaterial benchmark case TPV6, which shall be performed without any form of regularization.We define a near side and a far side of the fault plane.
The far side has a wave speed reduction of 60 % and a density reduction of 20 %.The fault is a planar strike-slip fault that reaches the free surface.Frictional sliding is governed by linear slip-weakening.
The exact values of the setup can be found in Table 4.The ADER-DG results are computed with   Fig. 12 shows the rupture time contour plots illustrating rupture front evolution.We find very good agreement of rupture speeds between the two methods.Only small differences in the arrival time occur in dip direction.These differences can be noted as well in the time series shown in Fig. 13.This receiver is located at the near side of the fault above the nucleation zone at the free surface.
Whereas the initial waves arrive at matching times for both methods, the FEM solution seems to evolve slower to the peak values.Thus, the fault fails later than in the ADER-DG solution, which would explain the delay in the rupture arrival time.This effect might be caused by the damping algorithm or the viscosity layer surrounding the fault implemented in FaultMod.The general shape of stress and particle velocity time series are, however, very similar.Only the relatively small, vertical components of shear stress and velocity show larger deviations as observable in Fig. 13 (b,e).The match of the normal components is again very good, although FaultMod computes slightly smaller peak values (Fig. 13 (c,f)).Note that the ADER-DG solution exhibits only small oscillations in all produced time series.
For the sake of completeness, time series recorded at the far side of the fault, larger distance of the nucleation zone, and at same depth as the previous station are presented in Fig. 14.At this station, rupture arrival time and slip rate peak seem to match better than for the near side station above the  nucleation zone.We conclude that for the well-posed bimaterial problem a very good agreement between ADER-DG and FaultMod is reached.
7 Rate-and state-dependent friction 330 The appropriate form of the constitutive law which describes the relationship between fault stress and slip along a fault plane is topic of intense research.A constitutive law should provide a correct description of the spatio-temporal variability of parameters involved in the earthquake rupture process, which proves to be difficult.Widely applied empirical friction laws are derived from smallscale experiments (Brace and Byerlee, 1966;Ruina, 1983;Ohnaka and Kuwahara, 1990;Di Toro et al., 2005;Niemeijer et al., 2010), although the question of a proper scaling of the parameters to seismic faulting is evoked.Developed friction laws differ in correlating weakening of tractions on    the fault surface to the slipped distance (slip-weakening), the particle velocity on the fault (ratedependent) and/or the history of the microstructural contacts (state-dependent).At high slip rates, non-linear physical weakening processes, as for example thermal pressurization of pore fluids, are thought to cause an observed extra-ordinary decrease of effective friction, referred to as fast velocityweakening.
The implementation of frictional sliding governed by slip-weakening friction into the ADER-DG algorithm has been presented in de la Puente et al. (2009); Pelties et al. (2012).Here we advance to the following implementations of rate-and state-dependent constitutive relationships.

Slow velocity friction
Let the frictional strength of the fault evolve as follows: where σ is the (constant) effective normal stress, L is a characteristic slip scale and µ 0 is the reference value of the friction coefficient at steady-state slip at reference sliding velocity V 0 .A rate-dependent part of the strength is proportional to the logarithm of the sliding velocity V and the frictional parameter a and is thought to reflect a thermally-activated Arrhenius process involving the breaking of atomic bonds at contact points bridging the sliding surface (Rice et al., 2001).A state-dependent part of the strength is proportional to the logarithm of the state variable Θ and the frictional parameter b and is thought to reflect the product of the true contact area and the intrinsic strength of those contacts.In the ageing law formulation the state variable evolves as whereas in the slip law formulation it is The implementation of rate-and state-dependent friction follows roughly Kaneko et al. (2008).A 5 stage Newton-Raphson algorithm is used to determine the value of the slip rate in every sub-time step obtained by the ADER integration scheme starting from the slip rate of the previous time step.
The state variable of the current time step is evaluated from the updated slip rate.Consecutively, the algorithm is repeated once.
We benchmark the ADER-DG scheme performance for the ageing law in the SCEC test case TPV101.We model a planar fault in an isotropic, linear elastic half-space.The model parametrization is given in Table 5.A transition layer of 3 km width in which the frictional properties continuously change from velocity-weakening to velocity-strengthening surrounds the central velocityweakening region of the fault.Outside of the transition region, the fault is velocity-strengthening.
The initial friction law parameter a (and thus, for self-consistency, the initial state variable) is space dependent, but the initial velocity and normal and horizontal shear stresses are uniform along the fault plane.The medium is initially moving with equal and opposite horizontal velocities of V 0 /2 on the two sides of the fault.Rupture is nucleated by imposing a horizontal shear traction perturbation that grows mathematically smoothly in time and space to its maximum amplitude ∆τ 0 over a finite time interval T in a region of the fault of radius R.

Fast velocity-weakening
In the same manner, we implement a rate-and state-dependent friction law with severe velocityweakening at slip rates faster than a characteristic velocity, as adopted by Ampuero and Ben-Zion; Dunham et al. (2011);Gabriel et al. (2013).Strong velocity-weakening has been proposed to fit results of laboratory experiments at fast slip velocities (see e.g.Di Toro et al. (2011) and references therein) and is predicted by a flash heating model (Rice, 2006).
The frictional strength is determined by the slip velocity (V ) and a state variable (Θ) as: where a is a positive coefficient quantifying a direct effect and V 0 is a reference slip rate.The state variable has units of slip and obeys the following evolution equation: Following Noda et al. (2009), we regularize the steady-state friction coefficient µ ss (obtained if Θ = 0) in the framework of rate-and-state friction in the slip law form: with µ s being the static friction coefficient, µ 0 a reference friction coefficient, V w a weakening velocity scale, µ w the fully weakened friction coefficient and b a positive coefficient quantifying an evolution effect.In this formulation, the transition between low velocity friction and strongly velocity-weakened friction is relatively smooth, which is favorable for numerical accuracy (Dunham et al., 2011).
We perform SCEC Test Case TPV103 which differs from slow velocity friction Test Case TPV101 as summarized in Table 6.In Fig. 16 we show near-perfect agreement of ADER-DG with other high-order dynamic rupture software packages based on a variety of numerical methods: the multidimensional spectral boundary integral code MDSBI (Dunham, 2008), the three-dimensional spectral elements code SPECFEM3D (Kaneko et al., 2008) and the three-dimensional spectral boundary integral element (SBIE) method implementation by Lapusta and Liu (2009).The agreement of the modeled dynamic rupture processes is depicted by the measured along-strike shear stress and slip rate on the fault outside the nucleation zone (Fig. 16    unknowns is used.The most important benefits of this method compared to existing methods are the absence of spurious modes in the slip rate and the tetrahedral discretization enabling geometrical complexity of fault surfaces.
Here, we briefly review a general algorithm for implementing fault dynamics in schemes based on flux functions, neglecting the specifics of the ADER-DG method.This way, our implementation might be readily applied to a wide range of related numerical schemes.In particular, the algorithm is independent of the chosen time integration method.However, the impact of the implementation of fault dynamics on numerical properties as dispersion, diffusion and accuracy is expected to differ.
In general, we assume that the Coulomb failure criterion needs to be fulfilled at any point on the fault, e.g. at an integration point of the DG method: with τ G being the shear traction and τ S the fault strength.We define the fault plane to be situated in the Cartesian y, z-plane and denote σ i as the stress tensor, with i = xx, yy, zz, xy, xz, yz as the normal and shear stresses, respectively.u, v, w are the particle velocities in x, y, z direction, respectively.An arbitrary fault side is defined as the +-side, the opposing one as the −- As long as the failure criterion is not fulfilled, the fault is locked and the physical stress and velocity variables, grouped into a vector Q, are assigned to be the solution of the standard elastic wave propagation equations obtained from the Godunov state leading to Q = Q G .
The set of equations to compute the Godunov states accounting for a material contrast across the fault, based on Toro (1999); LeVeque (2002); de la Puente et al. (2009); Pelties et al. (2012), is given by: where ρ denotes density, c p P-wave velocity and c s S-wave velocity.Note that σ yy , σ zz and σ yz are associated to zero wave speeds and do not contribute to the solution of the Riemann problem.
Thus, we do not compute them.
In case |τ G | = τ S the fault fails and slip is declared.Then, we impose the shear stresses to be and use Eq.A3 to obtain the slip rates In the case of on-going slip, we impose v and w individually for the + and − side as than the latter.This explains why the effect of the Prakash-Clifton regularization is much less apparent if L = 0.01 m.We note that, because ADER-DG's intrinsic numerical dissipation depends on resolution (h and O), the frequencies that are well regularized for a given resolution can be more unstable at a finer resolution.Our analysis of regularization time scales provides a criterion to set regularization and resolution parameters in ADER-DG simulations of bimaterial rupture in the ill-posed regime.
Fig. 1 (a) depicts the asymmetric unstructured mesh discretizing the computational domain and a snapshot of the absolute particle velocity on the dipping fault plane during the simulation.The ADER-DG on-fault receivers, such as the station shown in Fig. 1 (b) and (c), show generally

Fig. 1 .
Fig. 1.(a) Asymmetric unstructured mesh discretizing the dipping fault setup of benchmark problems TPV10 and TPV11 for SeisSol, and absolute slip rate |v| on the fault plane at t = 7.4 s of TPV10 (subshear scenario).(b) Along-dip slip rate and (c) along-dip shear stress of an on-fault receiver at 7.5 km along-strike and 12 km along-dip of TPV10.The ADER-DG solution is shown in black, the FEM comparison solution in red.

Fig. 2 .
Fig. 2. Supershear and free surface effects in TPV11.(a) Rupture front contours on the fault plane every 0.5 s.(b) Normal stress, (c) along-dip slip rate and (d) along-dip shear stress at an on-fault receiver located 0 km along-strike and 1.5 km along-dip.The ADER-DG solution is shown in black, the FEM comparison solution in red.

Fig. 3 .
Fig. 3. (a) Asymmetric ground motion in the surrounding of the dipping fault at t ≈ 8 s in TPV10.The color scale indicates the magnitude of absolute particle velocity at the surface, arrows and surface morphology reflect its direction.(b) and (c) are horizontal ground velocity time series at seismic stations located 3 km away from the fault trace and 12 km along-strike on the foot wall and hanging wall, respectively.(d) and (e) are the vertical ground velocities at the same stations.The ADER-DG solution is shown in black, the FEM comparison solution in red.

Fig. 4
Fig. 4 demonstrates the complexity of the rupture propagation caused by the small-scale heterogeneous background stress.Considering the given complexity of the model we find good agreement between the two compared methods.The good match can also be confirmed by Figs. 5 and 6 where the time series of two receivers located on the fault plane are shown.Fig. 7 depicts a snapshot of the horizontal component of slip rate at 5.5 s across the entire fault plane and indicates the location of the receiver in Fig. 5.In Fig. 5 (c) we observe a sharp feature at ∼ 5.5 s in the horizontal slip

Fig. 4 .
Fig. 4. Rupture front every 0.5 s of problem TPV17 with heterogeneous initial stress conditions.The ADER-DG rupture front is indicated in black, the FEM comparison solution in red.

Fig. 5 .
Fig.5.Shear stresses and slip rates for TPV17 at a receiver located at 9 km depth and −9 km along-strike distance.A sharp feature caused by an interface wave can be identified in the along-strike slip rate (c) at t ∼ 5.5 s.The ADER-DG solution is shown in black, the FEM comparison solution in red.

Fig. 6 .
Fig.6.Shear stresses and slip rates for TPV17 at a receiver located at 9 km deth and 9 km along-strike distance.The ADER-DG solution is shown in black, the FEM comparison solution in red.

Fig. 7 .
Fig. 7. Along-strike slip rate in m/s for TPV17 at t = 5.5 s, the instance in time at which the interface wave which originated at the free surface approaches the fault receiver indicated by an arrow and shown in Fig. 5.

For
comparison we choose the results of FaultMod with 100 m and 50 m node spacing.As expected, we see that the 100 m results of FaultMod are closer to the ADER-DG solution with gap and FaultMod's 50 m results are closer to the ADER-DG without gap solution.First, we discuss the rupture time contour plots of both methods and resolutions on the main fault in Fig.9(a) and on the branch in Fig.9 (b).We find very good agreement of early rupture evolution departing from the nucleation area.Clear differences occur along the junction point at along-strike distance 0 km.Whereas both FaultMod simulations stop shortly after the junction point, except some slip at the fault bottom and at the free surface, the ADER-DG simulation with gap continues further along the main fault.The ADER-DG solution without gap continues only for a short distance along the main fault, stops there a bit earlier than the three other solutions, but exhibits early rupture near the junction point.This observation can be linked to earlier rupture initiation at the branch (see Fig.9 (b)).The branch rupture of the ADER-DG solution with gap starts slightly later similar to the FaultMod solutions as the gap needs to be overcome first.With increasing distance along strike, the ADER-DG solution with gap propagates slower than all other solutions, however, only along the branch.The concentration of rupture fronts along-strike distance > 0 km on the main fault for ADER-DG without gap (concentrated blue lines in Fig.9 (a)) is simply the result of a smooth, spontaneous rupture arrest in the branch (as opposed to abrupt arrest by a barrier).Further, minor differences in the compared rupture time contour plots occur at the end of the seismogenic zone and at the free surface along the branching fault.

Fig. 8 .
Fig. 8. Zoom view of the junction point of TPV15 (2D example) to illustrate parallelization concept and discretization strategy.

Fig. 9 .
Fig. 9. Rupture contours of benchmark TPV15 on (a) the main fault and (b) the branch fault.The junction point is located at along-strike distance 0 km.ADER-DG solutions are shown in black (with gap) and blue (without gap), the FEM comparison solutions in red (50 m discretization) and green (100 m discretization).

Fig. 10 .
Fig. 10.Shear stresses and slip rates for TPV15 at a receiver on the branch fault at 9 km along-strike and 7.5 km depth.ADER-DG solutions are shown in black (with gap) and blue (without gap), the FEM comparison solutions in red (50 m discretization) and green (100 m discretization).

Fig. 11 .
Fig. 11.Shear and normal stresses of TPV15 recorded at the main fault in 2 km along-strike distance to the junction point and 7.5 km depth.Differences are due to differences in geometrical setups as discussed in the main text.ADER-DG solutions are shown in black (with gap) and blue (without gap), the FEM comparison solutions in red (50 m discretization) and green (100 m discretization).

Fig. 12 .
Fig. 12. Rupture front every 0.5 s in problem TPV6 with bimaterial conditions.The ADER-DG solution is shown in black, the FEM comparison solution in red.

Fig. 13 .
Fig. 13.Stresses and velocities of TPV6 at a surface station located on the near side of the fault at km along-strike.The ADER-DG solution is shown in black, the FEM comparison solution in red.

Fig. 14 .
Fig. 14.Stresses and velocities of TPV6 at a station located on the far side of the fault at −12 km along-strike and 7.5 km depth, which corresponds to the hypocenter depth.The ADER-DG solution is shown in black, the FEM comparison solution in red.
Fig. 15 a) illustrates this initial setup and the lateral mesh coarsening around the fault plane.Figs. 15 b) and c) illustrate the nearperfect agreement of ADER-DG and FEM in along-strike shear stress and slip rate measured at the hypocenter throughout the nucleation period.

Fig. 15 .
Fig. 15.(a) Nucleation and initial frictional parameters of SCEC Test Case TPV101.A velocity-weakening fault smoothly transitions to a surrounding velocity-strengthening material, by adapting the frictional parameter a. Nucleation is achieved by prescribing a space and time dependent circular stress perturbation.Along-strike (b) shear stress and (c) slip rate during nucleation at the hypocenter, located at 0 km along-strike and 7.5 km along-dip.The ADER-DG solution is shown in black, the FEM comparison solution in red.

Fig. 16 .
Fig. 16.Fast velocity-weakening benchmark TPV103 results.Comparison of ADER-DG (black) with the highorder dynamic rupture software packages MDSBI, SPECFEM3D and a SBIE implementation.Along-strike (a) shear stress and (c) slip rate on the fault at 12 km along-strike and 3 km along-dip.Along-strike (b) shear stress and (d) slip rate during nucleation at the hypocenter at strike 0 km, dip 7.5 km.The ADER-DG solution is shown in black, the comparison solutions of MDSBI in green, SPECFEM3D in blue and SBIE in purple.
side.Corollary, ∆v = v + − v − and ∆w = w + − w − are the slip rates.Eq.A1 prescribes that theshear traction |τ G | = (σ G xy + σ 0 xy ) 2 + (σ G xz + σ 0 xz ) 2 is always lower or equal to the fault strength τ S = µ f (σ G xx + σ 0 xx ), where µ f denotes the friction coefficient, which may be a function of slip, slip rate, state evolution variables and other frictional parameters.Variables denoted by superscript 0 are the corresponding initial values.

Fig. 17 .
Fig. 17.Influence of Prakash-Clifton regularization on ADER-DG solution of the ill-posed bimaterial problem TPV7.Time series at a surface receiver located on the far side off the fault above the nucleation zone for different values of the characteristic regularization length scale L with unit meters.

Table 6 .
SCEC Test Case TPV103 simulation parameters which differ from Test Case TPV101.