A unified first order hyperbolic model for nonlinear dynamic rupture processes in diffuse fracture zones
Alice-Agnes Gabriel, Duo Li, Simone Chiocchetti, Maurizio Tavelli, Ilya Peshkov, Evgeniy Romenski, Michael Dumbser
AA unified first order hyperbolic model for nonlinear dynamic ruptureprocesses in diffuse fracture zones
A.-A. Gabriel , D. Li , S. Chiocchetti , M. Tavelli ,I. Peshkov , E. Romenski , M. Dumbser Earthquake fault zones are more complex, both geometrically and rheologically, than an idealised in-finitely thin plane embedded in linear elastic material. To incorporate nonlinear material behaviour, nat-ural complexities and multi-physics coupling within and outside of fault zones, here we present a firstorder hyperbolic and thermodynamically compatible mathematical model for a continuum in a gravi-tational field which provides a unified description of nonlinear elasto-plasticity, material damage andof viscous Newtonian flows with phase transition between solid and liquid phases. The fault geometryand secondary cracks are described via a scalar function ξ ∈ [0, 1] that indicates the local level of mate-rial damage. The model also permits the representation of arbitrarily complex geometries via a diffuseinterface approach based on the solid volume fraction function α ∈ [0, 1]. Neither of the two scalar fields ξ and α needs to be mesh-aligned, allowing thus faults and cracks with complex topology and the useof adaptive Cartesian meshes (AMR). The model shares common features with phase-field approaches,but substantially extends them. We show a wide range of numerical applications that are relevant for dy-namic earthquake rupture in fault zones, including the co-seismic generation of secondary off-fault shearcracks, tensile rock fracture in the Brazilian disc test, as well as a natural convection problem in moltenrock-like material. Multiple scales, multi-physics interactions and nonlinearities govern earthquake source processes, rendering the un-derstanding of how faults slip a grand challenge of seismology [26, 58]. Over the last decades, earthquake rupturedynamics have been commonly modeled as a sudden displacement discontinuity across a simplified (potentially het-erogeneous) surface of zero thickness in the framework of elastodynamics [3]. Such earthquake models are commonlyforced to distinguish artificially between on-fault frictional failure and the off-fault response of rock. Here, we modelnatural fault damage zones [12, 56] by adopting a diffuse crack representation. Ludwig-Maximilians-Universität München, Theresienstr. 41, 80333 München, Germany, email Laboratory of Applied Mathematics, University of Trento, Via Mesiano 77, 38123, Trento, Italy Free University of Bolzano, Universitätplatz 1 , 39100 Bolzano, Italy Sobolev Institute of Mathematics, 4 Acad. Koptyug Avenue, 630090, Novosibirsk, Russia,Accepted to
Phil. Trans. R. Soc. A , doi.org/10.1098/rsta.2020.0130 a r X i v : . [ phy s i c s . g e o - ph ] O c t n recent years, the core assumption that faults behave as infinitely thin planes has been challenged [90]. Effortscollapsing the dynamics of earthquakes to single interfaces may miss important physical aspects governing fault-system behaviour such as complex volumetric failure patterns observed in recent well-recorded large and small earth-quakes [11,76] as well as in laboratory experiments [62]. However the mechanics of fault and rupture dynamics in gen-eralized nonlinear visco-elasto-plastic materials are challenging to incorporate in earthquake modeling. Earthquakespropagate as frictional shear fracture of brittle solids under compression along pre-existing weak interfaces (faultzones), a problem which is mostly unsolvable analytically. For numerical modeling, dynamic earthquake rupture isoften treated as a nonlinear boundary condition in terms of contact and friction, coupled to seismic wave propaga-tion in linear elastic material. The evolving displacement discontinuity across the fault is defined as the earthquakeinduced slip. Typically, the material surrounding the fault is assumed to be linear, isotropic and elastic, with all non-linear complexity collapsed into the boundary condition definition of fault friction [29, e.g.], which take the form ofempirical laws describing shear traction bounded by the fault strength. In an elastic framework, high stress concen-trations develop at the rupture front. The corresponding inelastic off-fault energy dissipation ( off-fault damage ) andits feedback on rupture propagation [45] can be modelled in form of (visco-)plasticity of Mohr-Coulomb or Drucker-Prager type [4, 82], a continuum damage rheology which may account for high strain rate effects [6, 51, 83], or explicitsecondary tensile and shear fracturing [17, 60, 93].Numerical modeling of crack propagation has been a long-standing problem not only in seismology but also incomputational mechanics. Emerging approaches in modeling fracture and rupture dynamics include phase-field andvarifold-based representations of cracks to tackle the major difficulty of the introduction of strong discontinuities inthe displacement field in the vicinity of the crack. Current state-of-the-art methods in earthquake rupture dynam-ics [39] require explicit fracture aligned meshing, thus, generally (with recent exceptions [61]) require fractures to bepredefined, and typically only permit small deformations. Using highly efficient software implementations of thisapproach large-scale earthquake modeling is possible [16, 40, 88]. Alternative spatial discretizations which allow rep-resenting strong discontinuities at the sub-element level, such as the eXtended Finite Element Method (XFEM) [57],introduce singularities when an interface intersects a cell, but are quite difficult to implement in an efficient manner.In distinction, diffuse interface approaches “smear out” sharp cracks via a smooth but rapid transition between in-tact and fully damaged material states [8, 19, 55]. Within various diffuse interface approaches, the most popular oneis the phase field approach , which allows to model complicated fracture processes, including spontaneous crack ini-tiation, propagation, merging, and branching, in general situations and for 3D geometries. Critical ingredients ofthe phase-field formulation are rooted in fracture mechanics, specifically by incorporating a critical fracture energy(from Griffith’s theory [37]), which is translated into the regularized continuum gradient damage mechanics [54]. Sev-eral theoretical methods have been recently proposed for shear fracture ( [78, e.g.] for mode III) which is dominatingearthquake processes. Phase-field models have also been successfully applied for brittle fracture in rock-like materi-als [95] on small scales (mm’s of slip).The material failure model discussed in this paper also belongs to the class of diffuse interface models in whichthe damaged material or a crack is considered as another “ phase ” of the material and represented by a continuousscalar field ξ ∈ [0, 1], called the damage variable . As in phase field approaches, a crack or failure front is representednot as a discontinuity of zero thickness but as a diffuse interface across which ξ changes continuously from 0 (in-tact material) to 1 (fully damaged material) resulting in gradual but rapid degradation of material stiffness. Despitethis conceptual similarity, the model developed here is very different from the phase field models. An importantfeature of the phase field models is the presence of the non-local regularization term ∼ (cid:107)∇ φ (cid:107) in the free energy,with φ being the phase field. Without such a regularization term, the numerical treatment of a phase field modelis problematic due to numerical instabilities and mesh dependency of the numerical solution. This indicates theill-posedness of the underlying governing PDEs, e.g. see [48, 50]. In contrast, the model developed here originating Faults are then idealised as two matching surfaces in unilateral contact not allowed to open or interpenetrate and typically implemented bysplitting the fault interface [18]. and is formulated based on the thermodynamicallycompatible continuum mixture theory [70, 72] which results in a first-order symmetric hyperbolic governing PDE sys-tem and thus is intrinsically well-posed, at least locally in time. Mathematical regularity of the model is supported bythe stability of the hereafter presented numerical results, including a mesh convergence analysis (see Sec. 3). Moregenerally, the developed model belongs to the class of Symmetric Hyperbolic and Thermodynamically Compatible(SHTC) equations [32, 34, 64, 75]. Apart from the PDE type used (the phase field models are formulated as second-order Allen-Cahn-type [1,36] or fourth-order Cahn-Hilliard-type [7,13,25] parabolic PDEs), there is also an importantconceptual difference between the developed mixture type approach and the phase field approaches. In the latter,the phase transformation is entirely controlled by the free energy functional, which usually consists of three terms: Ψ ( ε , φ , ∇ φ ) = Ψ ( ε , φ ) + Ψ ( φ ) + Ψ ( ∇ φ ), where ε is the small elastic strain tensor, Ψ is the elastic energy which com-prises a degradation function, Ψ is the damage potential (usually a double-well potential but also single-well po-tentials are used [47]), and Ψ is the regularization term. In our approach, only an energy equivalent to Ψ ( ε , φ ) isused [74, 81], while the phase-transition is described in the context of irreversible thermodynamics and is controlledby a dissipation potential which is usually a highly nonlinear function of state variables [63,64]. Yet, it is important toemphasize that the irreversible terms controlling the damage are algebraic source terms (no space derivatives), whichdo not affect the differential operator of the model. This greatly simplifies the discretization of the differential terms inthe governing PDE, but nevertheless requires an accurate and robust stiff ordinary differential equation solver [14, 81]for the source terms. Since the governing PDE system of our theory contains only first-order derivatives in spaceand time, it is possible to use explicit time-stepping in the numerical integration [81]. In contrast, the second andfourth-order phase field PDEs require the use of an implicit time discretization [36], which is more difficult to imple-ment and may not have advantage over explicit methods if the time step is dictated by the physical time scales, suchas in strongly time dependent processes, e.g. fracture dynamics and wave propagation. We note that a hyperbolicreformulation of phase-field models is possible as recently proposed in [44].Alternatively, variational views on fracture mechanics can describe crack nucleation intrinsically without a priorifailure criteria [27,53]. Accounting for microscopic surface irregularities or line defects can be achieved by combininga sharp interface approach with advanced tools of differential geometry such as curvature varifolds [31]. These ideascan be seen as a natural extension of the pioneering Griffith’s theory [37] with cracks being represented by almosteverywhere differentiable surfaces and evolving Griffith’s energies to account for curvature effects. In this context, weremark that the model presented here by no means is a complete fracture model. In specific situations requiring avery accurate prediction of the fracture process the merely constitutive capabilities of the present model may not besufficient. Instead, accounting explicitly for the energy accumulating at the irregularities of the crack surface (e.g., atcorners and cusps) or the dynamics of microscopic defects near the crack tip might be required. In the here presentedfirst-order hyperbolic diffuse interface framework, this can be achieved by taking into account higher gradients of thestate variables such as curvature and torsion in the form of independent state variables [15, 66]. The continuum model for damage of solids employed in this paper consists of two main ingredients. The first ingre-dient is the damage model proposed by Resnyansky, Romenski and co-authors [69, 74] which is a continuous damagemodel with a chemical kinetics-type mechanism controlling the damage field ξ ∈ [0, 1] ( ξ = ξ = Non-local terms can be introduced in our theory if it is physically motivated, e.g. see [15, 66] For example, relaxation times may change over several orders of magnitude across the diffuse interface zone.
3y the same constitutive relations (relaxation functions), but have different characteristic time scales, e.g. see [81].The second ingredient is the Eulerian finite strain elastoplasticity model developed by Godunov and Romenski in the1970s [33,35,73]. It was recently realized by Peshkov and Romenski [65] that the same equations can also be applied tomodeling viscous fluid flow, as demonstrated by Dumbser et al in [23] and thus, this model represents a unified formu-lation of continuum fluid and solid mechanics. In the following we shall refer to it as Godunov–Peshkov–Romenski(GPR) model. Being essentially an inelasticity theory, the GPR model provides a unified framework for continuousmodeling of potentially arbitrary rheological responses of materials, and in particular of inelastic properties of thedamaged material. This, in turn, can be used for modeling of complex frictional rheology in fault zones in geomateri-als, see Sec. 3. For further details on the GPR model, the reader is referred to [10, 23, 43, 65, 75].Our diffuse interface formulation for moving nonlinear elasto-plastic solids of arbitrary geometry and at large strainis given by the following PDE system in Eulerian coordinates: ∂ t α + v k ∂ k α = ∂ t ¯ ρ + ∂ k ( ¯ ρ v k ) =
0, (1a) ∂ t ( ¯ ρ v i ) + ∂ k (cid:161) ¯ ρ v i v k + α p δ ik − ασ ik (cid:162) = ¯ ρ g i , (1b) ∂ t A ik + ∂ k ( A im v m ) + v m ( ∂ m A ik − ∂ k A im ) = − θ − ( τ ) E A ik , (1c) ∂ t J k + ∂ k ( v m J m + T ) + v m ( ∂ m J k − ∂ k J m ) = − θ − ( τ ) E J k , (1d) ∂ t ξ + v k ∂ k ξ = − θ E ξ , (1e) ∂ t ( ¯ ρ S ) + ∂ k (cid:161) ¯ ρ Sv k + ¯ ρ E J k (cid:162) = ¯ ρ ( α T ) − (cid:161) θ − E A ik E A ik + θ − E J k E J k + θ E ξ E ξ (cid:162) ≥
0, (1f) ∂ t ( ¯ ρ E ) + ∂ k (cid:161) v k ¯ ρ E + v i ( α p δ ik − ασ ik ) (cid:162) = ¯ ρ g i v i , (1g)where we use the Einstein summation convention over repeated indices and ∂ t = ∂ / ∂ t , ∂ k = ∂ / ∂ x k . Here, (1a) is theevolution equation for the colour function α that is needed in the diffuse interface approach (DIM) as introducedin [10, 80] for the description of solids of arbitrary geometry ( α = α = ρ = αρ and (1a) is the mass conservation law with ρ being the material density; (1b) is the momentum conservation lawand v i is the velocity field; (1c) is the evolution equation for the distortion field A = [ A ik ], which is the main field inthe GPR model and can be viewed as the field of local basis triads representing the deformation and orientation ofan infinitesimal material element [23, 65, 66]; (1d) is the evolution equation for the specific thermal impulse J k , de-scribing the heat conduction in the matter via a hyperbolic (non Fourier–type) model; (1e) is the evolution equationfor the material damage variable ξ ∈ [0, 1], where ξ = ξ = E = E ( ρ , S , v , A , J , ξ ): p = ρ E ρ is the thermodynamicpressure, σ = [ σ ik ] = [ σ eik + σ tik ] is the stress tensor with contributions to the mechanical stress due to tangential[ σ eik = − ρ A ji E A jk ] and thermal stress [ σ tik = ρ J i E J k ] (note that σ e in not necessary trace-free), and T = E S is the tem-perature. The total mechanical stress tensor is defined as Σ = [ Σ ik ] = [ − p δ ik + σ e ik ], where δ ik is the Kronecker delta.With a state variable in the subscript of the energy, we denote partial derivatives, e.g. E ρ = ∂ E / ∂ρ , E A i j = ∂ E / ∂ A i j ,etc. Furthermore, g i is the gravitational acceleration vector. Also, because we are working in an Eulerian frame ofreference, we need to add transport equations of the type ∂ t λ + v k ∂ k λ = E = E ( ρ , S , v , A , J , ξ ). This potential then generates the fluxes (reversible time evolution) and source terms (irreversible Global deformation can not be restored from the local triad since they represent only local deformation and thus, incompatible deformation. E = E + E + E , decomposing the energy into a contribution from the microscale E , the mesoscale E and the macroscale E . The individual contributions read as follows: E = K (cid:161) − ρ / ρ (cid:162) /(2 ρ ) + c v T (cid:161) ρ / ρ (cid:162) (cid:161) e S / c v − (cid:162) + H ( T − T c ) h c , (2)where ρ and T are the reference mass density and temperature, h c is the latent heat, T c is the critical temperatureat which phase transition occurs, H ( T ) is the Heaviside step function, c v is the heat capacity at constant volume.As a proof of concept, we added the last term in (2) and present a demonstration example of the model’s capabilityto deal with solid-fluid phase transition (melting/solidification) in Section 5 of the supplementary material. Yet, thiscorresponds to a simplified (time-independent) modeling of phase transition and will be improved in the future. Also, K ( ξ ) = λ ( ξ ) + µ ( ξ ) is the bulk modulus, λ ( ξ ) and µ ( ξ ) are the two Lamé constants that are functions of the damagevariable ξ specified, following [69], as λ ( ξ ) = K I K D / ˜ K − µ I µ D /(3 ˜ µ ), µ ( ξ ) = µ I µ D / ˜ µ , (3)where the subscripts I and D denote ‘ intact ’ and ‘ damaged ’ respectively, K I = λ I + µ I , K D = λ D + µ D , ˜ K = ξ K I + (1 − ξ ) K D , ˜ µ = ξµ I + (1 − ξ ) µ D , and it is assumed that the elastic moduli of the intact material λ I , µ I and of the fully damagedmaterial λ D , µ D are known.The macro-scale energy is the specific kinetic energy E = v i v i . Finally, E reads E = c s ˚ G i j ˚ G i j + c h J i J i , (4)where c s ( ξ ) = (cid:112) µ ( ξ )/ ρ is the shear sound speed and c h is related to the speed of heat waves in the medium (alsocalled the second sound [67], or the speed of phonons). ˚ G ik = G ik − G j j δ ik is the deviator of the Finger (or metric)tensor G ik = A ji A jk that characterizes the elastic deformation of the medium.The dissipation in the system includes three irreversible processes that raise the entropy: the strain relaxation (orshear stress relaxation) characterized by the scalar function θ ( τ ) > τ , theheat flux relaxation characterized by θ ( τ ) > τ , and the chemical kineticslike process governing the transition from the intact to damaged state and controlled by the function θ in (1e).The main idea of the diffuse interface approach to fracture is to consider the material element as a mixture ofthe intact and the fully damaged phases. These two phases have their own independent material parameters andclosure relations, such as functions characterizing the rate of strain relaxation. The strain relaxation approach in theframework of the unified hyperbolic continuum mechanics model [23, 65] represented by the evolution equation forthe distortion field A allows to assign potentially arbitrary rheological properties to the damaged and intact states. Inparticular, the intact material can be considered as an elastoplastic solid, while the damaged phase can be a fluid, e.g.a Newtonian fluid (see Sec. 33.3) or viscoplastic fluid, which can be used for modeling of in-fault friction, for example.Yet, in this paper, we do not use an individual distortion evolution equation for each phase, but we employ the mixtureapproach [69,74], and we use a single distortion field representing the local deformation of the mixture element, whilethe individual rheological properties of the phases are taken into account via the dependence of the relaxation time τ on the damage variable ξ as follows: τ = ((1 − ξ )/ τ I + ξ / τ D ) − , (5)where τ I and τ D are shear stress relaxation times for the intact and fully damaged materials, respectively, which areusually highly nonlinear functions of the parameters of state. The particular choice of τ I and τ D that is used in thispaper reads τ I = τ I exp( α I − β I (1 − ξ ) Y ), τ D = τ D exp( α D − β D ξ Y ), (6)5here Y is the equivalent stress, while τ I , α I , β I , τ D , α D , β D are material constants. In this work, the stress norm Y iscomputed as Y = A Y s + B Y p + C , (7)where Y s = (cid:112) Σ dev Σ )/2, with dev Σ = Σ − (tr Σ /3) I , is the von Mises stress and Y p = tr Σ /3 accounts for thespherical part of the stress tensor. The choice A = B = C =
0, gives Y = Y s , that is, the von Mises stress, while otherchoices of coefficients in Eq. (7) are intended to describe a Drucker–Prager-type yield criterion.Note that to treat the damaged state as a Newtonian fluid, it is sufficient to take τ D = const (cid:191)
1, see Sec. 33.3 or [23].Non-Newtonian rheologies can be also considered if the proper function for τ D ( Y ) is provided, e.g. see [43]. Thus, thefunction θ = τ c s ( ξ ) /3 | A | − is taken in such a way as to recover the Navier-Stokes stress tensor with the effectiveshear viscosity η = ρ τ c s in the limit τ (cid:191) τ I = ∞ . By this means, all numerical examples presented throughout Sec.3 follow the rheologicalformulation given by θ with varying parametrisation.The transition from the intact to the fully damaged state is governed by the damage variable ξ ∈ [0, 1] satisfyingthe kinetic-type equation (1e), ˙ ξ = − θ E ξ , with the source term depending on the state parameters of the medium(pressure, stress and temperature). In particular, the rate of damage θ is defined as θ = θ (1 − ξ )( ξ + ξ (cid:178) ) (cid:163) (1 − ξ ) ( Y / Y ) a + ξ ( Y / Y ) (cid:164) , (8)where ξ (cid:178) , Y and Y , a are constants. ξ (cid:178) is usually set to ξ (cid:178) = − in order to trigger the growth of ξ with the initialdata ξ =
0. We note that similar to the chemical kinetics, the constitutive functions of the damage process drive thesystem towards an equilibrium that is not simply defined as E ξ =
0, but as θ E ξ =
0, e.g. see [63]. As a result, the overallresponse of the material subject to damage (i.e., its stress-strain relation, see also [81]) is defined by the interplay ofboth irreversible processes; (i) the degradation of the elastic moduli controlled by (8) and (ii) the inelastic processesin the intact and damaged phases controlled by (5) and (6). In the numerical experiments carried out in Sec. 33.2, thedamage kinetics ξ also strongly couple with strain relaxation effects, by means of Eq. (5). The function θ , governingthe rate of the heat flux relaxation, is taken as θ ( τ ) = τ c h ρ T that yields the classical Fourier law of heat conductionwith the thermal conductivity coefficient κ = τ c h in the stiff relaxation limit ( τ → a in Eq. (6). By this means, the rate of growth θ of the damage variable ξ will activateas a switch when Y reaches the Y threshold. In this specific case, Y can be chosen equivalently to a yield stress.Also, the sensitivity to tensile stresses can be modelled by resorting to techniques that are routinely used in scienceand engineering, e.g., using the Drucker–Prager yield criterion to compute Y . In the Brazilian tensile fracture exam-ple in Sec. 33.2, β I , D are set to zero as the complex stress-dependent mechanisms they control are not necessary forachieving the desired material behaviour. Controlling the relaxation time of the damaged state ( τ D ) can be useful formodelling friction within a natural fault zone: if a very low relaxation time is chosen, which can be easily achieved bysetting τ D = − s, α D = β D =
0, the fault will exert no tangential stresses on the surrounding intact rock, as if it werefilled with an inviscid fluid. Specific frictional regimes and (time-dependent) plastic effects can be described by prop-erly choosing the relaxation times τ I , D (via τ I D , α I , D , β I , D ), which in general may require more complex automaticoptimisation strategies. 6 b) intact host rock -σ n τ L fc W fc ψτ max xy fault core host rock 1, cases (i),(ii)fully elastic host rock 2, case (iii)spontaneous off-fault damage low velocity fault rock, cases (ii),(iii) (a) Figure 1: (a) Typical strike-slip fault zone structure showing a multiple fault core with associated damage zone in a quartzofelds-pathic country rock (from [56]). (b) Sketch of the GPR model setup for 2D in-plane right-lateral shear fracture under compressionused throughout Sec.33.1. In light grey we depict the prescribed fault core of length L f c and width W f c which is fully damaged( ξ =
1) and embedded in intact host rock ( ξ = In this section we present a variety of numerical applications of the GPR model relevant for earthquake rupture andfault zones. The governing PDE system (1) is solved using the high performance computing toolkit
ExaHyPE [68],which employs an Arbitrary high order derivative (ADER) discontinuous Galerkin (DG) finite element method in com-bination with an a posteriori subcell finite volume limiter on space time adaptive Cartesian meshes (AMR). For details,the reader is referred to [81] and to [9, 10, 22–24, 89, 94] and references therein.
In the following, we explore the GPR diffuse fault zone approach extending the modeling of dynamic earthquakerupture beyond treatment as a discontinuity in the framework of elastodynamics. Fig. 1 illustrates the model setupcorresponding to the geological structure of a typical strike-slip fault zone. Dynamic rupture within the ’fault core’is governed by a friction-like behaviour achieved by time-dependent modulation of the shear relaxation time τ D ofthe fault core’s fully damaged material. At the onset of frictional yielding, the shear relaxation time ( τ D ) decreasesexponentially as in (6) with a time-dependent β (cid:48) D . The temporal evolution of β (cid:48) D is modulated at a constant rate duringrupture as β (cid:48) D ( t ) = β D min (1, max (0, 1 − C t )) where C and β D are constant. Visco-elastic slip accumulates across thediffuse fault core coupled to either fully elastic wave propagation or Drucker-Prager type damage in the host rock. i) Kinematic self-similar Kostrov-like crack. We first model a kinematically driven non-singular self-similar shearcrack analog to Kostrov’s solution for a singular crack [46] to study the relation between fault slip, slip rate and shearstress in comparison to traditional approaches, while imposing tractions here avoids the full complexity of frictionalrupture dynamics. The 2D setup [20, e.g.] assumes a homogeneous isotropic elastic medium (Table S2, c s = c p / (cid:112) σ n =
40 MPa and shear stress τ =
20 MPa. An in-plane right-lateral shear crack is driven by prescribing the (sliding) friction µ f as linearly time-dependent weakening: µ f ( x , t ) = max{ f d , f s − ( f s − f d )( v r t − | x | )/ R c }, with process zone size R c =
250 m, rupture speed v r = f s = f d = C =
10 reproduces the propagatingshear crack in the reference solution. Thus, β (cid:48) D evolves linearly from β D to 0 during rupture.7 DER-DG (GPR model)SEM σ x y (d) (e)0 km2468(a) (b) (c) ADER-DG (GPR Model) SEM t=4.0 s v e l o c i t y ( m / s ) Figure 2:
Comparison of the self-similar Kostrov-like crack of the diffuse GPR model (ADER-DG, p = W f c = L f c = SEM2DPACK ( p = h = t =
4s (see also Animation S1), and zoom into the rupture tip.
We assume a fully-damaged fault core ( ξ =
1) of prescribed length L f c =
20 km and width W f c =
100 m embeddedin a continuum material resembling intact elastic rock ( ξ =
0) as illustrated in Fig. 2a. Both, the fault core and thesurrounding host rock are treated as the same continuum material besides their differences in ξ . The GPR specificmaterial parameters are detailed as ‘host rock 1’ (here, λ D = λ I , µ D = µ I ) in Table S1 in the supplementary material.The model domain is of size 70 km ×
70 km bounded by Dirichlet boundary conditions and employs a statically re-fined mesh surrounding the fault core. The domain is discretised into hierarchical Cartesian computational grids,spaced h = h =
311 m at the second refinement level (Table S3). We use polynomialdegree p = p + =
13 subcells in each spatial dimension. Fig.2a-ccompares slip, slip rate and shear traction during diffuse crack propagation modeled with the GPR model to a spectralelement solution assuming a discrete fault interface spatially discretised with h =
100 m with
SEM2DPACK [2]. TheGPR model analog captures the kinematics (i.e., stress drop and fault slip) of the self-similar singular Kostrov crack aswell as the emanated seismic waves (Fig. 2d,e and Animation S1), while introducing dynamic differences on the scaleof the diffuse fault (zoom-in in Fig. 2d). Slip velocity (Fig. 2a) remains limited in peak, similar to planar fault mod-eling with off-fault plastic deformation [30]. Fault slip (Fig. 2b) appears smeared out at its onset, yet asymptoticallyapproaches the classical Kostrov crack solution. Similarly, shear stresses (Fig. 2c) appear limited in peak and more dif-fuse, specifically with respect to the secondary peak associated with the passing rupture front. Importantly, (dynamic)stress drops are comparable to the expectation from fracture mechanics for a plane shear crack (even though peak and8ynamic level appear shifted). At the crack tip, we observe an initial out-of-plane rotation within the fault core leadingto a localised mismatch in the hypocentral region and at the onset of slip across the fault. The GPR model approachesthe analytical solution, as illustrated for increasing polynomial degree p in the electronic supplement (Fig. S1). ii) Spontaneous dynamic rupture. We next model spontaneous dynamic earthquake rupture in a 2D version [20]of the benchmark problem TPV3 [39] for elastic spontaneous rupture propagation defined by the Southern CaliforniaEarthquake Center. Our setup resembles the kinematic model (Fig. 1a) including the time-dependent choice of β (cid:48) D ( t )with C =
10 with an important distinction: we assume a low-rigidity fault core (‘low velocity fault rock’ in Table S1) bysetting P-wave and S-wave velocity in the fault core 50% lower, i.e. λ ( ξ ) and µ ( ξ ) are decreased by 30%, with respect tothe intact rock. A 50% reduction of seismic wave speeds matches natural fault zone observations. The thickness of thelow velocity fault rock unit equals the thickness of the fault core itself where ξ =
1. The surrounding country rock isagain parameterised as fully elastic with the ‘host rock 1’ GPR parametrisation (Table S1). The fault core is L f c =
30 kmlong and W f c =
100 m wide, the domain size is 40 km ×
40 km, initial loading is σ y y = −
120 MPa and σ x y =
70 MPa. Thecomputational grid is spaced h = h =
177 m at the second refinement level (Table S3).Fig. 3 compares, similar to the kinematic case, the diffuse low-rigidity fault ADER-DG GPR results to an elastic discretefault interface spectral element solution. Fault slip rates (Fig. 3a) are limited in peak and are now clearly affected bysmaller scale dynamic complexity, e.g. out-of-plane crack rotation and wave reflections, within the fault core. Faultslip (Fig. 3b) asymptotically resembles the discontinuous, elastic solution. Shear stresses (Fig. 3c) are smeared outand shifted, but capture (dynamic) stress drops, similar to the kinematic model in i) . We note that residual shearstress levels remain higher potentially reflecting oblique shear developing within the diffuse fault core and/or viscousbehaviour within the fault core. The diffuse fault core slows down the emitted seismic waves, while amplifying sharpvelocity pulses (Fig. 3d,e and Animation S2) aligning with observational findings [79]. The GPR model successfullyresembles frictional linear-slip weakening behaviour [42] within the fault core by defining: µ f ( x , t ) = max{ f d , f s − ( f s − f d ) δ ( x , t )/ D c }, with slip-weakening distance D c = f s = f d = δ ( x , t ) denotes here the diffuse slip within the fault core and is measured as the difference of displacements at itsadjacent boundaries. Rupture is not initiated by an overstressed patch, which would be inconsistent with deformingmaterial, but as a kinematically driven Kostrov-like shear-crack with v r = t = rheological fault core properties , and not the friction law, control important crack dynamics such as rupture speedin our diffuse interface modeling, cf. [41]. A comparison of results assuming a further reduction of fault rock wavespeeds to 60% is discussed in the supplementary material. iii) Dynamically generated off-fault shear cracks. Localized shear banding is observed in the vicinity of naturalfaults spanning a wide spectrum of length scales [56], and contributes to the energy balance of earthquakes. We modeldynamically generated off-fault shear cracks by combining the spontaneous dynamic rupture model embedded in‘low velocity fault rock’ with ‘host rock 2’ outside the fault core (Table S1, µ D = µ I , λ D = λ I + µ I − µ D ) in(3)). ’Host rock 2’ is governed by Drucker–Prager yielding [21, 82, 92] as given by Eq. (7), with A = (cid:112) B = sin( π /18),and C = − cos( π /18) ·
95 MPa. The model domain size is 20 km ×
15 km spatially discretised with h =
800 m at thecoarsest mesh level (Table S3). We here use dynamic adaptive mesh refinement (AMR) with two refinement levels andrefinement factor r = h =
89 m. Fig. 4a illustrates spontaneous shear-cracking in the extensional quadrants of the main fault core,where the passing rupture induces a dynamic bimaterial effect [84]. While previous models [82] based on ideal plastic-ity without damage accumulation numerically capture the formation of single sets of shear bands in Drucker-Pragertype off-fault material induced by dynamic rupture propagation across a main fault, we here observe the formationof two conjugate sets of shear fractures: Cracks are distributed around two favourable orientations (Fig. 4b). Spacingand length of these shear deformation bands [52, 60] may depend on GPR material parametrization ( Y , β D , cohe-sion, internal friction angle, etc., see Table S1 and [81]) as well as on the computational mesh and will motivate future9 x y v e l o c i t y ( m / s ) ADER-DG (GPR Model)SEM4 km6810 (a) (b) (c)(d) t=3.1 s
ADER-DG (GPR Model) SEM x(km) x(km) y ( k m ) (e) Figure 3:
Computational results for the 2D TPV3 dynamic rupture problem. Comparison of the diffuse interface GPR model(ADER-DG, p = W f c = L f c = SEM2DPACK solution ( p = h = t = analysis, as in Sec. 33.2. High particle velocity is associated with the strong growth of off-fault shear stresses near thefault tip shifting from the propagation direction of the main rupture [28]. We observe the dynamic development ofinterlaced conjugate shear faulting (Animation S3) resembling recent high-resolution imaging of earthquakes [76]. The GPR framework can be applied to capture tensile fracture, important for earthquake nucleation processes andthe microscale of fault zone fracture and damage. We now show that our model is able to correctly describe thefracture mechanisms observed in laboratory settings. Specifically, we reproduce the experimental results of [38] whichinvolve the compression of a rock disc along its diameter (a so-called Brazilian test). In this case the disc presents acentral slit with a given orientation, which influences the early stages of the failure of the rock sample. The test iscarried out in two space dimensions on a square computational domain centered at the origin and with sidelength2.2 units. The interface of the disc is defined by setting α = ρ = , µ I = λ I = µ D = λ D = a) (b) σ % x (km) y ( k m ) t= 5.0 s - - - Figure 4:
Off-fault shear cracks spontaneously generated in the extensional quadrants of dynamic earthquake rupture (TPV3) inthe GPR model (ADER-DG, p = W f c = L f c = t = ξ illustrating the fault core initialized as fully damaged (cf. Fig. 2a) and the propagating secondary off-fault cracks. (b)Polar diagram of the statistical orientation of off-fault shear cracks. The two dominant orientations are ≈ ◦ and ≈ ◦ . Themaximum compressive stress ( σ ) has an orientation angle of 65.3 ◦ . θ = Y =
10 MPa, Y = a = τ I = s, τ D = − s, β I = β D =
0. For | y | > Y = Y =
100 TPa to model unbreakable clamps. Thermal effects are neglected. For this test, thecoefficients of the Drucker–Prager equivalent stress formula (7) are A = B = C = − p =
3) scheme on a uniform Cartesian mesh of 192 by 192 cells,showing good agreement with the experimental data. For a detailed mesh refinement study, see the supplementarymaterial.
Seismic fault slip velocities and low thermal conductivity of rock can lead to the formation of veins of molten rock(pseudotachylytes), which are thought of as an unambiguous indicator of earthquake deformation, however, are notcommon features of active faults [77]. In our model, the phase transition between solid and liquid occurs simply viathe definition of the total energy by adding the contribution of the latent heat for T > T c , see (2), and by modifyingthe relaxation time for T > T c . More precisely, in this example, the relaxation time τ is not computed according to(5) and (6) but is considered constant (time-independent) in the solid state and is computed in terms of the dynamicviscosity η as τ = ηρ c s for the molten state ( T > T c ) treated as a Newtonian fluid. Also, in this example, θ has tobe taken as θ = τ c s ( ξ ) /3 | A | − , see the result of the asymptotic analysis presented in [23]. In the supplementarymaterial of this paper we validate our simple approach for phase transition for a well-known benchmark problemwith exact solution, namely the Stefan problem, see [5]. The obtained results clearly show that the proposed modelcan properly deal with heat conduction and phase transition between liquid and solid phases.Next, we show the capability of the GPR model to describe also the motion of viscous fluids under the influenceof gravity. The stresses acting on faults are key initial conditions for earthquakes and seismic fault dynamics, but are11oorly known. At very long time scales, these initial conditions are governed by plate tectonics and mantle convec-tion, which is included in the GPR model as a special case. We therefore simulate a rising bubble in molten rock-likematerial. In the following, we use SI units. The critical temperature is set to T c = h c = g = (0, − η =
20. We furthermore set theremaining parameters to ρ = γ = p = · , c v = c s = α = λ = T = v i = A = I , J = p = − (cid:107) g (cid:107) ρ y and a hot circular bubble of radius R = x c = (0, 0) with atemperature increase of ∆ T =
200 for (cid:107) x − x c (cid:107) ≤ R . The domain is Ω = [ −
2, 2] × [ −
1, 3] and simulations are carried outuntil t = p =
3) scheme on a mesh of 200 ×
200 elements. For comparison, we run two simula-tions, one with the GPR model presented in this paper and another simulation with the compressible Navier-Stokesequations, which serves as a reference solution for the GPR model in the viscous fluid limit. The computational resultsare depicted in Fig. 6, where we can observe an excellent agreement. This demonstrates that the model presented inthis paper is able to correctly describe also natural convection in molten material when T > T c . We have presented a unified hyperbolic model of inelasticity that incorporates finite strain elastoviscoplasticity andviscous fluids in a single PDE system, coupled with a hyperbolic model for continuous modeling of damage, includingbrittle and ductile fracture as particular cases. The governing equations are formulated in the Eulerian frame and via adiffuse interface approach permit arbitrary geometries of fractures and material boundaries without the necessity ofgenerating interface-aligned meshes. We emphasize that the presented diffuse interface approach is not merely a wayto regularize otherwise singular problems as posed by earthquake shear crack nucleation and propagation along zero-thickness interfaces, but potentially allows to fully model volumetric fault zone shearing during earthquake rupture,which includes spontaneous partition of fault slip into intensely localized shear deformation within weaker (possiblycohesionless/ultracataclastic) fault-core gouge and more distributed damage within fault rocks and foliated gouges.The model capabilities were demonstrated in several 2D examples related to rupture processes in earthquake faultzones. We compare kinematic, fully dynamic and off-fault damage GPR diffuse rupture to models employing the tra-ditional elasto-dynamic viewpoint of a fault, namely a planar surface across which slip occurs. We show that the con-
Figure 5:
Crack formation in a rock-like disc under vertical load (Brazilian test) for different angles of the pre-damaged area.Comparison of the contour colors of the damage variable ξ obtained in the numerical simulations of the GPR model with thecracks observed in experiments. The simulation results are overlaid on top of the photographs from [38]. From left to right: 45 ◦ ,60 ◦ and 90 ◦ . Only the regions of the disc where α > igure 6: Temperature contours for the rising bubble problem in molten rock-like material at time t =
4. Solution obtained withthe GPR model (left) and Navier-Stokes reference solution (right). The melting temperature is set to T c = tinuum model can reproduce and extend classical solutions, while introducing dynamic differences (i) on the scaleof pre-damaged/low-rigidity fault zone, such as out-of-plane rupture rotation, limiting peak slip rates, non-frictionalcontrol of rupture speed; and (ii) on the scale of the intact host rock, such as conjugate shear cracking in tensile lobesand amplification of velocity pulses in the emitted wavefield. A natural next step is to combine the successful microfracture laboratory scale Brazilian tests with dynamic rupture to span the entire scales of fault zone fracture. The GPRparameters for the host rock and fault zone rock material can also be calibrated to resemble natural rock, as e.g. West-erly granite [49]. Also, using the computational capabilities of the model’s ExaHyPE implementation, one can studyrelated effects on ground shaking (see [80, 81] for GPR modeling of 3D seismic wave propagation with complex to-pography) and detailed 3D fault zone models [85, 86, 91, cf.] including trapped/head waves interacting with dynamicrupture [41]. Inelastic bulk processes are important during earthquake rupture [87, e.g.,], but also in between seismicevents, including off-fault damage and its healing, dynamic shear localization and interseismic delocalization, andvisco-elasto-plastic relaxation. Since the unified mathematical formulation proposed in this paper is able to describeelasto-plastic solids as well as viscous fluids, future work will also concern the study of fully coupled models of dy-namic rupture processes triggered by mantle convection and plate tectonics. Extensions to non-Newtonian fluids willbe considered, as well as to elasto-plastic saturated porous media, see e.g. the recent work presented in [43, 71]. Wealso plan more detailed investigations concerning the onset of melting processes in shear cracks. Finally, we note thatthe material failure is due to the accumulation of microscopic defects (micro-cracks in rocks or dislocations in crys-talline solids). It is thus interesting to remark that the distortion field being the field of non-holonomic basis triadsprovides a natural basis for further development of the model towards a micro-defects-based damage theory. Thiscan be achieved via concepts of the Riemann-Cartan geometry, such as torsion discussed in [66].
Data Accessibility.
ExaHyPE is free software hosted at . The presented numerical examples willbe accessible and reproducible at https://gitlab.lrz.de/exahype/ExaHyPE-Engine and https://github.com/TEAR-ERC/GPR2DR . 13 unding.
This research has been supported by the European Union’s Horizon 2020 Research and Innovation Pro-gramme under the projects
ExaHyPE , grant no. 671698, ChEESE, grant no. 823844 and TEAR, grant no. 852992. MDand IP also acknowledge funding from the Italian Ministry of Education, University and Research (MIUR) via the De-partments of Excellence Initiative 2018–2022 attributed to DICAM of the University of Trento (grant L. 232/2016) andthe PRIN 2017 project
Innovative numerical methods for evolutionary partial differential equations and applications .SC was also funded by the Deutsche Forschungsgemeinschaft (DFG) under the project DROPIT, grant no. GRK 2160/1.ER was also funded within the framework of the state contract of the Sobolev Institute of Mathematics (project no.0314-2019-0012). AG also acknowledges funding by the German Research Foundation (DFG) (grants no. GA 2465/2-1,GA 2465/3-1), by KAUST-CRG (grant no. ORS-2017-CRG6 3389.02) and by KONWIHR (project NewWave). Computingresources were provided by the Institute of Geophysics of LMU Munich [59] and the Leibniz Supercomputing Centre(project no. pr63qo).
In memoriam.
This paper is dedicated to the memory of Anne-Katrin Gabriel (*March 7, 1957 †July 25, 2020) whosecreativity and curious mind will live on – in science with her daughter Alice.
References [1] M. Ambati, T. Gerasimov, and L. De Lorenzis. A review on phase-field models of brittle fracture and a new fasthybrid formulation.
Computational Mechanics , 55(2):383–405, 2015. (Cited on page 3.)[2] J.-P. Ampuero. SEM2DPACK, a spectral element software for 2D seismic wave propagation and earthquake sourcedynamics, v2.3.8, 2012. https://sourceforge.net/projects/sem2d/. (Cited on page 8.)[3] D.J. Andrews. Rupture velocity of plane strain shear cracks.
Journal of Geophysical Research , 81(32):5679–5687,1976. (Cited on page 1.)[4] D.J. Andrews. Rupture dynamics with energy loss outside the slip zone.
Journal of Geophysical Research: SolidEarth , 110(B1), 2005. (Cited on page 2.)[5] H. D. Baehr and K. Stephan.
Heat and Mass Transfer . Springer, 2011. (Cited on pages 11 and 3.)[6] H.S Bhat, A.J. Rosakis, and C.G. Sammis. A micromechanics based constitutive model for brittle failure at highstrain rates.
Journal of Applied Mechanics , 79(3), 2012. (Cited on page 2.)[7] M. J. Borden, T.J.R. Hughes, C.M. Landis, and C.V. Verhoosel. A higher-order phase-field model for brittle fracture:Formulation and analysis within the isogeometric analysis framework.
Computer Methods in Applied Mechanicsand Engineering , 273:100–118, 2014. (Cited on page 3.)[8] B. Bourdin, G.A. Francfort, and J.J. Marigo. Numerical experiments in revisited brittle fracture.
Journal of theMechanics and Physics of Solids , 48(4):797–826, 2000. (Cited on page 2.)[9] H.J. Bungartz, M. Mehl, T. Neckel, and T. Weinzierl. The PDE framework Peano applied to fluid dynamics: Anefficient implementation of a parallel multiscale fluid dynamics solver on octree-like adaptive Cartesian grids.
Computational Mechanics , 46:103–114, 2010. (Cited on page 7.)[10] S. Busto, S. Chiocchetti, M. Dumbser, E. Gaburro, and I. Peshkov. High order ADER schemes for continuummechanics.
Frontiers in Physiccs , 8:32, 2020. (Cited on pages 4 and 7.)[11] Y. Cheng, Z.E. Ross, and Y. Ben-Zion. Diverse volumetric faulting patterns in the San Jacinto fault zone.
Journalof Geophysical Research: Solid Earth , 123(6):5068–5081, 2018. (Cited on page 2.)1412] F.M. Chester, J.P. Evans, and R.L. Biegel. Internal structure and weakening mechanisms of the San Andreas fault.
Journal of Geophysical Research: Solid Earth , 98(B1):771–786, 1993. (Cited on page 1.)[13] L.R. Chiarelli et al. Comparison of high order finite element and discontinuous Galerkin methods for phase fieldequations: Application to structural damage.
Computers & Mathematics with Applications , 74(7):1542–1564,2017. (Cited on page 3.)[14] S. Chiocchetti and C. Müller. A Solver for Stiff Finite-Rate Relaxation in Baer-Nunziato Two-Phase Flow Models.
Fluid Mechanics and its Applications , 121:31–44, 2020. (Cited on page 3.)[15] S. Chiocchetti, I. Peshkov, S. Gavrilyuk, and M. Dumbser. High order ADER schemes and GLM curl cleaning for afirst order hyperbolic formulation of compressible flow with surface tension.
Journal of Computational Physics ,2020. (Cited on page 3.)[16] Y. Cui et al. Physics-based seismic hazard analysis on petascale heterogeneous supercomputers. In
SC ’13:Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analy-sis , pages 1–12, 2013. (Cited on page 2.)[17] L.A. Dalguer, K. Irikura, and J.D. Riera. Simulation of tensile crack generation by three-dimensional dynamicshear rupture propagation during an earthquake.
Journal of Geophysical Research: Solid Earth , 108(B3), 2003.(Cited on page 2.)[18] S.M. Day, L.A. Dalguer, N. Lapusta, and Y. Liu. Comparison of finite difference and boundary integral solutionsto three-dimensional spontaneous rupture.
Journal of Geophysical Research: Solid Earth , 110(B12), 2005. (Citedon page 2.)[19] R. de Borst, J.J.C. Remmers, A. Needleman, and M.A. Abellan. Discrete vs smeared crack models for concretefracture: bridging the gap.
International journal for numerical and analytical methods in geomechanics , 28(7-8):583–607, 2004. (Cited on page 2.)[20] J. de la Puente, J.-P. Ampuero, and M. Käser. Dynamic rupture modeling on unstructured meshes using a dis-continuous Galerkin method.
Journal of Geophysical Research: Solid Earth , 114(B10), 2009. (Cited on pages 7and 9.)[21] D.C. Drucker and W. Prager. Soil mechanics and plastic analysis or limit design.
Quarterly of Applied Mathemat-ics , 10(2):157–165, 1952. (Cited on page 9.)[22] M. Dumbser, F. Fambri, M. Tavelli, M. Bader, and T. Weinzierl. Efficient implementation of ADER discontinuousGalerkin schemes for a scalable hyperbolic PDE engine.
Axioms , 7:63, 2018. (Cited on page 7.)[23] M. Dumbser, I. Peshkov, E. Romenski, and O. Zanotti. High order ADER schemes for a unified first order hy-perbolic formulation of continuum mechanics: Viscous heat-conducting fluids and elastic solids.
Journal ofComputational Physics , 314:824–862, 2016. (Cited on pages 4, 5, 6, 7, and 11.)[24] M. Dumbser, O. Zanotti, R. Loubère, and S. Diot. A posteriori subcell limiting of the discontinuous Galerkin finiteelement method for hyperbolic conservation laws.
Journal of Computational Physics , 278:47–75, 2014. (Cited onpage 7.)[25] L. O. Eastgate et al. Fracture in mode I using a conserved phase-field model.
Physical Review E , 65(3):1–10, 2002.(Cited on page 3.) 1526] D.W. Forsyth, T. Lay, R.C. Aster, and B. Romanowicz. Grand challenges for seismology.
EOS, Transactions Ameri-can Geophysical Union , 90(41):361–362, 2009. (Cited on page 1.)[27] G.A. Francfort and J.-J. Marigo. Revisiting brittle fracture as an energy minimization problem.
Journal of theMechanics and Physics of Solids , 46(8):1319–1342, 1998. (Cited on page 3.)[28] L.B. Freund.
Dynamic Fracture Mechanics . Cambridge University Press, 1990. (Cited on page 10.)[29] A.A. Gabriel, J.-P. Ampuero, L.A. Dalguer, and P.M. Mai. The transition of dynamic rupture styles in elastic mediaunder velocity-weakening friction.
Journal of Geophysical Research: Solid Earth , 117(B9), 2012. (Cited on page 2.)[30] A.A. Gabriel, J.-P. Ampuero, L.A. Dalguer, and P.M. Mai. Source properties of dynamic rupture pulses with off-faultplasticity.
Journal of Geophysical Research: Solid Earth , 118:4117–4126, 2013. (Cited on page 8.)[31] M. Giaquinta, P.M. Mariano, G. Modica, and D. Mucci. Ground states of simple bodies that may undergo brittlefractures.
Physica D: Nonlinear Phenomena , 239(15):1485–1502, 2010. (Cited on page 3.)[32] S.K. Godunov. An interesting class of quasilinear systems.
Dokl. Akad. Nauk SSSR , 139(3):521–523, 1961. (Citedon page 3.)[33] S.K. Godunov.
Elements of mechanics of continuous media . Nauka, 1978. (in Russian). (Cited on page 4.)[34] S.K. Godunov, T.Y. Mikhailova, and E.I. Romenskii. Systems of thermodynamically coordinated laws of conser-vation invariant under rotations.
Siberian Mathematical Journal , 37:690–705, 1996. (Cited on page 3.)[35] S.K. Godunov and E.I. Romenski. Nonstationary equations of nonlinear elasticity theory in Eulerian coordinates.
Journal of Applied Mechanics and Technical Physics , 13(6):868–884, 1972. (Cited on page 4.)[36] H. Gomez and K.G. van der Zee. Computational Phase-Field Modeling. In
Encyclopedia of Computational Me-chanics Second Edition , pages 1–35. John Wiley & Sons, Chichester, UK, 2017. (Cited on page 3.)[37] A.A. Griffith. The phenomena of rupture and flow in solids.
Philosophical Transactions of the Royal Society ofLondon. Series A, Containing Papers of a Mathematical or Physical Character , 221(582-593):163–198, 1921. (Citedon pages 2 and 3.)[38] H. Haeri, K. Shahriar, M. Fatehi Marji, and P. Moarefvand. Experimental and numerical study of crack prop-agation and coalescence in pre-cracked rock-like disks.
International Journal of Rock Mechanics and MiningSciences , 67:20 – 28, 2014. (Cited on pages 10 and 12.)[39] R.A. Harris et al. A suite of exercises for verifying dynamic earthquake rupture codes.
Seismological ResearchLetters , 89(3):1146–1162, 2018. (Cited on pages 2 and 9.)[40] A. Heinecke et al. Petascale high order dynamic rupture earthquake simulations on heterogeneous supercom-puters. In
SC’14: Proceedings of the International Conference for High Performance Computing, Networking,Storage and Analysis , pages 3–14. IEEE, 2014. (Cited on page 2.)[41] Y. Huang, J.-P. Ampuero, and D. V. Helmberger. Earthquake ruptures modulated by waves in damaged fault zones.
Journal of Geophysical Research: Solid Earth , 119(4):3133–3154, 2014. (Cited on pages 9 and 13.)[42] Y. Ida. Cohesive force across the tip of a longitudinal-shear crack and griffith’s specific surface energy.
Journal ofGeophysical Research , 77(20):3796–3805, 1972. (Cited on page 9.)1643] H. Jackson and N. Nikiforakis. A numerical scheme for non-Newtonian fluids and plastic solids under the GPRmodel.
Journal of Computational Physics , 387:410–429, 2019. (Cited on pages 4, 6, and 13.)[44] D. Kamensky, G. Moutsanidis, and Yu. Bazilevs. Hyperbolic phase field modeling of brittle fracture: Part I - Theoryand simulations.
Journal of the Mechanics and Physics of Solids , 121:81–98, 2018. (Cited on page 3.)[45] M. Kikuchi. Inelastic effects on crack propagation.
Journal of Physics of the Earth , 23(2):161–172, 1975. (Cited onpage 2.)[46] B.V. Kostrov. Self–similar problems of propagation of shear cracks.
Journal of Applied Mathematics and Mechan-ics , 28(5):1077–1087, 1964. (Cited on page 7.)[47] V.I. Levitas, H. Jafarzadeh, G.H. Farrahi, and M. Javanbakht. Thermodynamically consistent and scale-dependentphase field approach for crack propagation allowing for surface stresses.
International Journal of Plasticity ,111:1–35, 2018. (Cited on page 3.)[48] T. Liebe and P. Steinmann. Theory and numerics of a thermodynamically consistent framework for geometricallylinear gradient plasticity.
International Journal for Numerical Methods in Engineering , 51(12):1437–1467, aug2001. (Cited on page 2.)[49] D.A. Lockner. A generalized law for brittle deformation of westerly granite.
Journal of Geophysical Research: SolidEarth , 103(B3):5107–5123, 1998. (Cited on page 13.)[50] E. Lorentz and A. Benallal. Gradient constitutive relations: numerical aspects and application to gradient dam-age.
Computer Methods in Applied Mechanics and Engineering , 194(50-52):5191–5220, dec 2005. (Cited onpage 2.)[51] V. Lyakhovsky, Y. Ben-Zion, and A. Agnon. A viscoelastic damage rheology and rate-and state-dependent friction.
Geophysical Journal International , 161(1):179–190, 2005. (Cited on page 2.)[52] X. Ma and A. Elbanna. Strain localization in dry sheared granular materials: A compactivity-based approach.
Physical Review E , 98(2):022906, 2018. (Cited on page 9.)[53] P. M. Mariano. Physical significance of the curvature varifold-based description of crack nucleation.
Rend. LinceiMat. Appl. , 21:215–233, 2010. (Cited on page 3.)[54] C. Miehe, L.M. Schänzel, and H. Ulmer. Phase field modeling of fracture in multi-physics problems. Part I. Bal-ance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids.
Computer Methodsin Appl. Mechanics and Engineering , 294:449–485, 2015. (Cited on page 2.)[55] H. Mirzabozorg and M. Ghaemian. Non-linear behavior of mass concrete in three-dimensional problems using asmeared crack approach.
Earthquake engineering & structural dynamics , 34(3):247–269, 2005. (Cited on page 2.)[56] T.M. Mitchell and D.R. Faulkner. The nature and origin of off-fault damage surrounding strike-slip fault zoneswith a wide range of displacements: A field study from the atacama fault system, northern chile.
Journal ofStructural Geology , 31(8):802 – 816, 2009. (Cited on pages 1, 7, and 9.)[57] N. Moës, J. Dolbow, and T. Belytschko. A finite element method for crack growth without remeshing.
Interna-tional Journal for Numerical Methods in Engineering , 46(1):131–150, 1999. (Cited on page 2.)1758] S. Nielsen. From slow to fast faulting: recent challenges in earthquake fault mechanics.
Philosophical Transac-tions of the Royal Society A: Mathematical, Physical and Engineering Sciences , 375(2103):20160016, 2017. (Citedon page 1.)[59] J. Oeser, H.P. Bunge, and M. Mohr. Cluster design in the Earth Sciences – Tethys. In
Int. Conference on HighPerformance Computing and Communications , pages 31–40. Springer, 2006. (Cited on page 14.)[60] K. Okubo et al. Dynamics, radiation, and overall energy budget of earthquake rupture with coseismic off-faultdamage.
Journal of Geophysical Research: Solid Earth , 124:11771–11801, 2019. (Cited on pages 2 and 9.)[61] K. Okubo, E. Rougier, Z. Lei, and H.S. Bhat. Modeling earthquakes with off-fault damage using the combinedfinite-discrete element method.
Computational Particle Mechanics , 2020. (Cited on page 2.)[62] F.X. Passelègue et al. Frictional evolution, acoustic emissions activity, and off-fault damage in simulated faultssheared at seismic slip rates.
Journal of Geophysical Research: Solid Earth , 121:7490–7513, 2016. (Cited on page 2.)[63] M. Pavelka, V. Klika, and M. Grmela.
Multiscale Thermo-Dynamics . De Gruyter, Berlin, 2018. (Cited on pages 3and 6.)[64] I. Peshkov, M. Pavelka, E. Romenski, and M. Grmela. Continuum Mechanics and Thermodynamics in the Hamil-ton and the Godunov-type Formulations.
Continuum Mechanics and Thermodynamics , 30:1343–1378, 2018.(Cited on page 3.)[65] I. Peshkov and E. Romenski. A hyperbolic model for viscous newtonian flows.
Continuum Mechanics and Ther-modynamics , 28:85–104, 2016. (Cited on pages 4 and 5.)[66] I. Peshkov, E. Romenski, and M. Dumbser. Continuum mechanics with torsion.
Continuum Mechanics andThermodynamics , 31(5):1517–1541, 2019. (Cited on pages 3, 4, and 13.)[67] V. Peshkov. Second sound in helium ii.
Journal of Physics (USSR) , 8:381, 1944. (Cited on page 5.)[68] A. Reinarz et al. Exahype: An engine for parallel dynamically adaptive simulations of wave problems.
ComputerPhysics Communications , 254:107251, 2020. (Cited on page 7.)[69] A. D. Resnyansky, E. I. Romensky, and N K Bourne. Constitutive modeling of fracture waves.
Journal of AppliedPhysics , 93(3):1537–1545, 2003. (Cited on pages 3 and 5.)[70] E. Romenski, D. Drikakis, and E. Toro. Conservative Models and Numerical Methods for Compressible Two-PhaseFlow.
Journal of Scientific Computing , 42(1):68–95, 2010. (Cited on page 3.)[71] E. Romenski, G. Reshetova, I. Peshkov, and M. Dumbser. Modeling wavefields in saturated elastic porous mediabased on thermodynamically compatible system theory for two-phase solid-fluid mixtures.
Computers & Fluids ,206:104587, 2020. (Cited on page 13.)[72] E. Romenski, A.D. Resnyansky, and E.F. Toro. Conservative hyperbolic formulation for compressible two-phaseflow with different phase pressures and temperatures.
Quarterly of Applied Mathematics , 65(2):259–279, 2007.(Cited on page 3.)[73] E. I. Romenski. Dynamic three-dimensional equations of the Rakhmatulin elastic-plastic model.
Journal ofApplied Mechanics and Technical Physics , 20(2):229–244, 1979. (Cited on page 4.)1874] E. I. Romenski. Deformation model for brittle materials and the structure of failure waves.
Journal of AppliedMechanics and Technical Physics , 48(3):437–444, 2007. (Cited on pages 3 and 5.)[75] E.I. Romenski. Hyperbolic systems of thermodynamically compatible conservation laws in continuum mechan-ics.
Mathematical and computer modelling , 28(10):115–130, 1998. (Cited on pages 3 and 4.)[76] Z.E. Ross et al. Hierarchical interlocked orthogonal faulting in the 2019 ridgecrest earthquake sequence.
Science ,366(6463):346–351, 2019. (Cited on pages 2 and 10.)[77] R.H. Sibson. Generation of pseudotachylyte by ancient seismic faulting.
Geophysical Journal International ,43(3):775–794, 1975. (Cited on page 11.)[78] R. Spatschek, E. Brener, and A. Karma. Phase field modeling of crack propagation.
Philosophical Magazine ,91(1):75–95, 2011. (Cited on page 2.)[79] P. Spudich and K. B. Olsen. Fault zone amplified waves as a possible seismic hazard along the Calaveras Fault incentral California.
Geophysical Research Letters , 28(13):2533–2536, 2001. (Cited on page 9.)[80] M. Tavelli et al. A simple diffuse interface approach on adaptive cartesian grids for the linear elastic wave equa-tions with complex topography.
Journal of Computational Physics , 386:158–189, 2019. (Cited on pages 4 and 13.)[81] M. Tavelli, E. Romenski, S. Chiocchetti, A.-A. Gabriel, and M. Dumbser. Space-time adaptive ADER discontinuousGalerkin schemes for nonlinear hyperelasticity with material failure.
Journal of Computational Physics , page109758, 2020. (Cited on pages 3, 4, 6, 7, 9, and 13.)[82] E.L. Templeton and J.R. Rice. Off-fault plasticity and earthquake rupture dynamics: 1. dry materials or neglect offluid pressure changes.
Journal of Geophysical Research: Solid Earth , 113(B9), 2008. (Cited on pages 2 and 9.)[83] M.Y. Thomas and H.S. Bhat. Dynamic evolution of off-fault medium during an earthquake: a micromechanicsbased model.
Geophysical Journal International , 214(2):1267–1280, 2018. (Cited on page 2.)[84] M.Y. Thomas, H.S. Bhat, and Y. Klinger.
Effect of Brittle Off-Fault Damage on Earthquake Rupture Dynamics ,chapter 14, pages 255–280. American Geophysical Union, 2017. (Cited on page 9.)[85] T. Ulrich et al. Coupled, physics-based modeling reveals earthquake displacements are critical to the 2018 Palu,Sulawesi Tsunami.
Pure and Applied Geophysics , 176(10):4069–4109, 2019. (Cited on page 13.)[86] T. Ulrich, A.-A. Gabriel, J.P. Ampuero, and W. Xu. Dynamic viability of the 2016 Mw 7.8 Kaik¯oura earthquakecascade on weak crustal faults.
Nature communications , 10(1):1–16, 2019. (Cited on page 13.)[87] T. Ulrich, A.-A. Gabriel, and E.H. Madden. Stress, rigidity and sediment strength control megathrust earthquakeand tsunami dynamics.
EarthArXiv , 2020. submitted. https://eartharxiv.org/s9263/. (Cited on page 13.)[88] C. Uphoff et al. Extreme scale multi-physics simulations of the tsunamigenic 2004 sumatra megathrust earth-quake. In
Proceedings of the international conference for high performance computing, networking, storage andanalysis , pages 1–16, 2017. (Cited on page 2.)[89] T. Weinzierl and M. Mehl. Peano - A traversal and storage scheme for octree-like adaptive Cartesian multiscalegrids.
SIAM Journal on Scientific Computing , 33:2732–2760, 2011. (Cited on page 7.)[90] C.A.J. Wibberley, G. Yielding, and G. Di Toro. Recent advances in the understanding of fault zone internal struc-ture: a review.
Geological Society, London, Special Publications , 299:5–33, 2008. (Cited on page 2.)1991] S. Wollherr, A.-A. Gabriel, and P.M. Mai. Landers 1992 ’reloaded’: Integrative dynamic earthquake rupture mod-eling.
Journal of Geophysical Research: Solid Earth , 124(7):6666–6702, 2019. (Cited on page 13.)[92] S. Wollherr, A.-A. Gabriel, and C. Uphoff. Off-fault plasticity in three-dimensional dynamic rupture simulationsusing a modal Discontinuous Galerkin method on unstructured meshes.
Geophysical Journal International ,214:1556–1584, 2018. (Cited on page 9.)[93] T. Yamashita. Generation of microcracks by dynamic shear rupture and its effects on rupture growth and elasticwave radiation.
Geophysical Journal International , 143(2):395–406, 2000. (Cited on page 2.)[94] O. Zanotti, F. Fambri, M. Dumbser, and A. Hidalgo. Space-time adaptive ADER discontinuous Galerkin finiteelement schemes with a posteriori subcell finite volume limiting.
Computers and Fluids , 118:204–224, 2015.(Cited on page 7.)[95] X. Zhang, S.W. Sloan, C. Vignes, and D. Sheng. A modification of the phase-field model for mixed mode crackpropagation in rock-like materials.
Computer Methods in Applied Mechanics and Engineering , 322:123–136, 2017.(Cited on page 2.) 20 upplementary Material for:
A unified first order hyperbolic model for nonlinear dynamicrupture processes in diffuse fracture zones material Y (Pa) Y (Pa) a α I β I α D β D µ I , D λ I , D host rock 1 1.8e22 1.0e20 32.5 36.25 0.0 36.25 1.0e-6 1.0 0.0host rock 2 1.8e8 1.0e10 32.5 36.25 0.0 36.25 1.0e-6 0.8571 0.6667low-velocity fault rock 1.8e22 1.0e20 32.5 36.25 0.0 12.687 1.75e-7 1.0 0.0 Table S1:
GPR material parameterisation for earthquake shear rupture models in section 3(a). Host rock 1 is the fully elastic ma-terial used in the full model domain of case (i), the kinematic Kostrov-like crack. The fully elastic host rock 1 is also surroundingthe fault core in case (ii), the spontaneous dynamic rup-ture example.
Host rock 2 is the Drucker-Prager type off-fault materialof case (iii) in which off-fault shear cracks nucleate and propagate. The (visco-)elastic low-velocity fault rock material is usedwithin the fault core in cases (ii) and (iii) and corresponds to the fully elastic host rock 1 with 50% decreased seismic velocities.Base relaxation times have been set as τ I = − · µ I and τ D = − · µ D . case ρ (kg/m ) c p (m/s) c s (m/s) f s f d D c or R c (m) σ y y (Pa) σ x y (Pa)(i) 2500 4000 2309 0.4 0.2 250 -40e6 20e6(ii) 2676 6000 3464 0.677 0.525 0.4 -120e6 70e6(iii) 2676 6000 3464 0.677 0.525 0.4 -120e6 70e6 Table S2:
Effective GPR material properties and reference friction parameters for section 3(a). Case (i) specifies the parameters ofthe self-similar Kostrov-like crack model. Cases (ii) and (iii) specify the parameters of the spontaneous dynamic rupture in fullyelastic and Drucker-Prager type plastic material. Note that the low-velocity fault rock material of the fault core in cases (ii) and(iii) has c p and c s
50% lower than host rock 1. case domainlength (m) domainwidth (m) fault corelength L f c (m) fault corewidth W f c (m) coarsest meshsize h (m) static refinementwidth (m) mesh refinementlevel; factor r (i) 70000 70000 20000 100 2800 1400 2; 3(ii) 40000 40000 30000 100 1600 800 2; 3(iii) 20000 15000 20000 100 800 400 2; 3 Table S3:
Geometry and computational mesh used in section 3(a). The domain is discretised into hierarchical Cartesian com-putational grids of two mesh refinement levels and using mesh refinement factor r =
3. In case (i) it is spaced h = h min =
311 m at the second refinement level in in the vicinity of the fault core, in case (ii) h = h min =
177 m, respectively, and in case (iii), which uses dynamic AMR, h =
800 m and h min =
89 m. Earthquake shear fracture animations
Supplementary Animation S1:
A map-view animation of the p=6 GPR Kostrov-like crack and seismic wave radia-tion model (case i) in terms of velocity is provided as supplementary material and accessible here: https://youtu.be/CBcbLeqaB_A
Supplementary Animation S2:
A map-view animation of the p=6 GPR dynamic rupture and seismic wavefieldTPV3 model (case ii) in terms of velocity is provided as supplementary material and accessible here: https://youtu.be/wx6-m6XS8C0
Supplementary Animation S3:
A map-view animation of the p=6 GPR co-seismic off-fault shear cracks sponta-neously generated in the extensional quadrants of dynamic earthquake rupture TPV3 model and seismic wavefield(case iii) in terms of velocity is provided as supplementary material and accessible here: https://youtu.be/95oEDZBqIiE
Figure S1:
Fault slip rate of the self-similar Kostrov-like crack (case i)modeled with the diffuse ADER-DG GPR method under varying poly-nomial degree p at 4 km hypocentral distance. As refer-ence the dis-crete fault spectral element SEM2DPACK solution is given. Low-rigidity fault zone effects on spontaneous dynamic rupture (TPV3model)
Figure S2:
Further reduction of P-and S-wave velocity of the rheol-ogy of the ‘low-velocity fault rock’for GPR TPV3 (case ii) results. (a,c)are taken from Fig. 3 (a,c) inthe main manuscript. (b,d) illus-trate a decrease in rupture speedand increase in peak slip ratewhile the stress drop is equiva-lent. All computational results aresolving the 2D TPV3 dynamic rup-ture problem of the diffuse inter-face GPR model using an ADER-DG ( p =
6, static AMR) schemewith the discrete fault spectral ele-ment SEM2DPACK solution ( p = h =
100 m), with (a,b) slip rate and(c,d) shear stress measured at in-creasing hypocentral distance.
Here we validate this simple approach on a well-known benchmark problem with exact solution, namely the Stefanproblem, see [5]. We consider the computational domain Ω = [0, 2] × [ − + y di-rection. The model parameters are set to γ = p = , ρ = α = , λ = , c s =
0. The initial densities arechosen as ρ L = x < ρ R = x ≥
0. Velocity, J and the pressure are initially set to zero and A is ini-tialized with the identity matrix. The critical temperature is chosen as T c =
900 and the latent heat is h c = . We nowsolve the GPR model until time t = × at constant density, seeFig. S3. The agreement is good, in particular considering the fact that in the GPR model density and pressure do notremain constant in time, as they are instead in the exact solution of the Stefan problem. The obtained results clearlyshow that the proposed model can properly deal with heat conduction and phase transition between liquid and solidphase H.D. Baehr and K. Stephan. Heat and Mass Transfer. Springer, 2011 T Exact solution of Stefan problemADER DG O4 (GPR model)
Figure S3:
Computational resultsobtained with the GPR model forthe Stefan problem at time t = T c = Here we report the results of a mesh convergence study carried out on the Brazilian disc test. We aim to show thatthe main directions in which cracks develop do not depend on the mesh resolution used in the computations. Thesimulations are carried out with a fourth order ADER-DG scheme with MUSCL-Hancock Finite Volume subcell lim-iting. We would like to stress that both the intensity of material damage and the problem geometry, as well as thelocal material properties, are represented by means of scalar fields, so that arbitrarily complex configurations can besimulated without requiring to use any ad-hoc meshing strategies. The setup for the test problem is shown in Fig.S4,where we plot the fields associated with the geometrical field α (which identifies the presence of solid material or ofair/vacuum) and with a material- type indicator. The use of two different sets of material properties is due to the needto model a clamping apparatus that we approximate as unbreakable. Note that no feature jump coincides with a cellboundary and in general the geometry can be skewed with respect to the computational grid, or even curved. Figure S4:
Computational resultsobtained with the GPR model forthe Stefan problem at time t = T c = In Fig.S5 we observe that a three-fold increase in the mesh resolution (from 64 cells per space dimension, to 96,to 128, to 192) allows to achieve better sharpness and definition of the cracks, but leaves their position essentiallyunaltered. For the last two steps, the crack topology converges to a stable configuration also at the points of contactbetween the clamps and the sample, which appear to be more sensitive to grid resolution with respect to the interiorof the disc. Furthermore, in no case we are able to observe any Cartesian mesh-imprinting artefacts in our solutions.4 igure S5:
Mesh convergence study for crack formation in a rock-like disc under vertical load (Brazilian test). We also comparethe contour colours of the damage parameter ξ obtained in the numerical simulations of the GPR model with the cracks observedin experiments. From left to right, the inclination angles of the pre-damaged slit are: 45 ◦ , 60 ◦ and 90 ◦ . From top to bottom, thenumber of grid cells in each dimension is: 64, 96, 128, and 192.. From top to bottom, thenumber of grid cells in each dimension is: 64, 96, 128, and 192.