A Diffuse Interface Method for Earthquake Rupture Dynamics Based on a Phase-Field Model

.

The spectral element method (SEM) has been a method of choice in the computational seismology community for simulating wave propagation in heterogeneous and homogeneous media (Komatitsch & Tromp, 1999).It aims to combine the geometrical flexibility of the FEM with the accuracy of spectral methods interpolating with high-order basis functions (e.g., Igel, 2017).The SEM is well suited for highly non-linear problems with non-smooth solutions, including simulations of dynamic rupture (Festa & Vilotte, 2006;Kaneko et al., 2008) using a split-node approach (Day et al., 2005) and hexahedral spectral elements for example, Galvez et al. (2014).SEM allows using non-linear off-fault plasticity (Gabriel et al., 2013) and continuum damage (Xu et al., 2015) but requires, similar to other established dynamic rupture modeling methods, to explicitly discretize fault discontinuities.
An alternative approach for representing a fault as a material discontinuity is the inelastic zone or stress glut method.Backus and Mulcahy (1976) termed stress glut the stress mismatch of a purely elastic medium and the true physical stress.The stress glut method was first established within the context of kinematic earthquake source descriptions, where the region of non-zero stress glut is the internal source (Bukchin, 1995;Chen et al., 2005;Clévédé et al., 2004;Jordan & Juarez, 2019, 2020, 2021;McGuire, 2017;McGuire et al., 2001McGuire et al., , 2002)).The stress glut (Andrews, 1999) and the thick fault zone (Madariaga et al., 1998) methods in the context of dynamic rupture modeling have been implemented in the finite difference method (Andrews, 1976(Andrews, , 1999;;Dalguer & Day, 2006;Herrendörfer et al., 2018;Madariaga et al., 1998;Preuss et al., 2019).There, the stress glut approximates the fault-jump conditions through inelastic increments to the stress components in an inelastic zone that is one grid cell wide.While the thick fault method leads to qualitative disagreement, the stress glut method produces qualitatively consistent results with the discontinuous reference solutions.However, it features inherently poor convergence with mesh refinement irrespective of the order of the finite difference approximation that was used (Dalguer & Day, 2006).While these inherent challenges have somewhat damped interest in the stress glut method for dynamic rupture earthquake modeling, we demonstrate in this study that a novel stress glut phase-field adaption can yield quantitatively consistent results to discrete fault reference solutions and empirical earthquake frictional behavior.
A classical phase-field approach has not yet been applied to fully dynamic earthquake rupture modeling.Phase-field approaches introduce a scalar phase-field, which varies between 0 and 1, to represent the degree of damage of the material (e.g., Bourdin et al., 2000).One major advantage of "field-based" approaches is that 10.1029/2023JB027143 3 of 40 fractures do not need to be explicitly meshed-thus enabling the simulation of spontaneous fracture development (e.g., Bourdin et al., 2008).Critical ingredients of the phase-field formulation are rooted in fracture mechanics, specifically by incorporating a critical fracture energy, which is translated into the regularized continuum sense of gradient damage mechanics (Miehe & Schänzel, 2014).For shear fracture, which is dominating earthquake processes, theoretical methods have been proposed (e.g., Spatschek et al., 2011) and applied for brittle fracture in rock-like materials under constant normal pressure (e.g., Fei & Choo, 2020, 2021;X. Zhang et al., 2017).Recently, this work has been extended to incorporate a rate-and state-dependent friction law in a promising antiplane quasi-dynamic phase-field model of fault growth and off-fault damage (Fei et al., 2023).
Here, we modify the concept of stress glut and apply it for the first time in a SEM.We combine the method with a spatially smooth and mesh-independent fault representation in a steady-state phase-field approach.We represent the fault geometry as the zero level set of a signed distance function (SDF).We demonstrate that in the phase-field framework, dynamic crack propagation can be handled as a standard multi-field problem by using conventional finite element methods.We note that the methodology described in this paper is not strongly tied to a continuous Galerkin or SEM but is, in principle, applicable to a discontinuous Galerkin approach for wave propagation (Dumbser & Käser, 2006;Wilcox et al., 2010).Our approach and its numerical implementation are explained in Sections 2 and 3. We verify our approach in Section 4 by performing kinematic and dynamic rupture benchmarks and by comparing our diffuse fault results to those of discrete fault modeling.We explore the flexibility of modeling dynamic rupture mesh independently by generalizing the fault geometry to inclined and curved planes not aligning with the prescribed computational mesh.In Section 5, we compare our approach against alternative diffuse crack models, discuss limitations, and anticipate future developments.

A Phase-Field Modified Stress Glut Approach
In this section, we formulate an implicit description of a diffuse fault geometry by means of the SDF.This description enables us to construct an in-fault reference frame that defines an embedded subdomain.In this sub-region, inelastic deformation can take place as a consequence of frictional yielding, in which the friction coefficient is a function of time or displacement.Using these ingredients, we present an extension of the stress glut method using the phase-field mathematical notion.

A Diffuse Fault Representation Using the Signed Distance Function
In this paper, we use the term "diffuse fault" to refer to a fault description of finite thickness.All applications developed in this study model earthquake slip on diffuse faults that are resolved by at least two spectral elements in width.Given a description of a fault as, for example, a parametric curve   =  (),  ∈ ℝ 2 ,  ∈ ℝ embedded in the 2D space  Ω ⊆ ℝ 2 , we construct an implicit model of the same geometry by defining a field  () ∈ ℝ,  ∈ Ω that satisfies the following properties: 1.At each point x ∈ Ω, |φ(x)| measures the Euclidean distance to the point on the curve x f (a) that is nearest to it in the same Euclidean sense, that is,: 2. The sign of the field φ(x) (denoted via sgn(φ)) is informally given by the side on which the coordinate x is with respect to the curve x f .More formally, given a value for the parameter a * = a * (x) that minimizes the Euclidean distance between the points x and x f (a * ), we can define a fault-normal vector  v =  −  ( * ) and a fault-tangent vector  ŵ =  ′  ( * ) , and arbitrarily but consistently assign  sgn(()) = sgn( ŵ1 v2 − ŵ2 v1) , the sign of the rotation from  v into  ŵ .A field that has these properties is called a SDF of the curve x f , and the original curve is partially recovered as the unordered level set  Γ = { ∶ () = 0} .In the following discussion, we will assume to only have access to the SDF and will forgo reference to the parametric curve x f (a) and its parameter a.In this regard, the method readily generalizes to a three-dimensional space embedding a fault as a two-dimensional manifold.
We define a right-handed orthonormal fault-local reference frame that is spanned by the normal vector n(x) = −∇φ(x) (which is a unit vector since |∇φ(x)| = 1 by definition) and a tangential vector t = [n 2 , −n 1 ].In this case, n points from the negative side of the fault to the positive side of the fault, but this is a rather arbitrary convention, much like the handedness of the fault-local coordinate system.
Here we use the SDF to extend the lower-dimensional fault interface to a finite sub-region Σ ⊂ Ω that is delineated by the ±δ level sets of the SDF, that is, Σ = {x ∈ Ω : −δ ≤ φ(x) ≤ +δ}.On account of the smooth nature of the SDF, it can represent a fault in a manner independent of the mesh resolution and orientation as long as the local curvature of the fault plane itself is well resolved.
To implement slip or slip rate dependent friction laws or evaluate results of time-dependent source descriptions, we project material displacement (u) and velocity (v) vectors onto the ±δ level sets.For a given coordinate x ∈ Σ, we compute two related coordinates on opposing sides of the fault as follows: the effective slip S is then calculated as and the effective slip rate  Ṡ is calculated similarly.This procedure generalizes the mesh-aligned stress glut implementation of Andrews (1999).The magnitudes of shear and normal tractions τ and σ n on the fault are expressed throughout the diffuse fault subdomain Σ as where   is the Cauchy stress tensor and the normal stress σ n is negative under compression.Figure 1 illustrates the geometric quantities introduced above, which are associated with the diffuse fault representation.
Note that the slip direction is derived from the evolving displacement field as a consequence of the embedded fault and its conditions.The displacement field evolves relative to the modified stress, where the shear direction and sign in fault local coordinates are inherited from the shear stress component of the stress state outside of the yield envelope.Yielding in our approach is described in the next section.

Yielding and Friction in the Diffuse Fault Stress Glut Approach
In this work, we assume a friction coefficient  =  (    Ṡ . . . ) that is a function of time t, position, slip, slip rate, and potentially other variables as well.Such a general description of the friction law encompasses the time-dependent Kostrov crack model and the linear slip weakening law that we use in this work to verify our method.We note that other frictional constitutive equations will be supported as well as, for example, the rate and state friction law (Dieterich, 1981;Ruina, 1983).A cohesionless frictional yield criterion τ c is stated as The truncation to negative values of normal stress effectively means that free slip (zero shear stress) conditions are applied under tension, in line with the fracture mechanics theory of Palmer and Rice (1973).Note that Day et al. (2005) describes an alternative treatment of jump conditions for tensile stresses, which may be explored in future developments of our approach.
A stress state outside of the yield envelope can be relaxed back onto it in the given direction of shear stress by applying a plastic correction, which we write as 10.1029/2023JB027143 5 of 40 (5) where the subscript f on the modified stress tensor denotes fault or friction, and the subscript e denotes elastic.
The map   →  for stress states that are outside of the yield envelope is achieved in plasticity models in a more subtle way by introducing a plastic strain increment of unknown magnitude and solving for said magnitude such that stress equals strength.This subtlety is not needed in the stress glut approach, as will be further clarified in the following.Under shear failure, the introduction of the stress limiter develops a transversely isotropic constitutive behavior with a plane of isotropy perpendicular to the direction n.See Sharples et al. (2016) for an extensive examination of the formulation and behavior of transversely isotropic materials in failure.In addition, we introduce an additional term in the yielding criterion in Equation 4. This change aims to include the implicit assumptions in traditional planar interface models.Its effects are further described in Section 4 and analyzed in Section 5.1.
We consider a modified yielding criterion that omits the contribution from the term G (∇(u • n)) • t.This simplification is motivated by 2D shear-driven deformation in planar Couette flow solutions (White, 2006).Importantly, this alteration is only applied when evaluating the yielding criterion, not during the elasticity update (see Algorithm 1).This formulation aims to emulate fault normal continuity at interfacial node pairs that exists in traditional dynamic rupture implementations for comparability of our diffuse fault representation to established models for mode II dynamic rupture.Then, such criterion becomes if YieldCriterionType = Volumetric then 10: else if YieldCriterionType = Interface then 12: end if 14: if F yc ≥ 0 then ⊳ Failure criterion reached 15: where G is the shear modulus.This interface yielding criterion may be interpreted as a modified Hooke's law that includes rotations in addition to strains within an infinitesimal continuum volume or as an alternative constitutive regularization of stresses within the inelastic zone, which limits a part of the shear stress components.
In the following, we refer to the yielding criterion introduced in Equation 4 as "volumetric yielding," while we refer to the inequality criterion in Equation 6as "interface yielding."The simplifying assumption of negligible change in fault normal displacements of the interface yielding may introduce numerical artifacts within a wide diffuse fault zone, which we further discuss in Appendix C. We will find that the volumetric yielding criterion is preferred for continuous fault zone representations throughout Section 4.
A major challenge associated with the classical stress glut method is the inherently sharp transition between on-fault and off-fault rheologies, which can lead to poor convergence properties and spurious oscillations, especially if the boundaries of the fault zone Σ intersect with grid cells (Dalguer & Day, 2006;Madariaga et al., 1998).This situation frequently occurs when modeling the fault independent of the mesh through level sets of the SDF, as described in this work.To address this difficulty, we define a time-invariant and smooth parameter ϕ ∈ [0, 1] based on the SDF, that is, ϕ = ϕ(φ), with ϕ(0) ≈ 1 and lim φ≫δ ϕ(φ) = 0. We suggest here to take a function ϕ(φ) of the form: where A, φ c are positive, nonzero parameters that influence the nature of the smooth transition from within the inelastic continuous fault zone to the elastic matrix of the host rock.In the following, we refer to the parameters A, φ c as "blending parameters."Equation 7 is motivated by the steady-state equilibrium profile obtained in thermodynamically derived phase-field models (Beckermann et al., 1999), where it describes the phase-field parameter variation normal to a given interface (Sun & Beckermann, 2007).
A stress tensor that is smoothly distributed over the domain Σ but approximately satisfying the yield limit Equation 4 everywhere can be redefined as The continuity conditions for both the traction components of stress and the fault normal displacement are implicitly enforced as they are integral parts of the continuum problem formulation.The shear stress correction is continuous by the phase-field approach.

Elastodynamics of Dynamic Rupture
The elastic stress tensor is given by the constitutive relation where   is the fourth order constitutive tensor, composed of the Lamé parameters λ, G;   is the second order unit tensor, and   is the strain tensor defined as the symmetric gradient of the displacement u: The dynamic momentum balance governs the wave-mediated evolution of friction on the fault and is expressed as where ρ is the density.The problem is closed and applying boundary conditions on the fault, further explained in Section 3 and model-specific initial conditions that are given in Section 4. In all models presented, we impose a free-surface boundary condition along the entire boundary of Ω, that is we enforce   =  on ∂Ω.

Spectral Elements for a Phase-Field Method
We use se2dr, a rupture dynamics extension of the stress glut method of Andrews (1999), implemented in the 2D wave propagation SEM se2wave using the high-level library PETSc (Abhyankar et al., 2018;Balay et al., 1997Balay et al., , 2019Balay et al., , 2020)), as our linear algebra backend.
Our implementation uses a structured quadrilateral mesh to discretize the domain Ω.The SEM nodal basis is given by a Lagrange polynomial, which in combination with a Gauss-Legendre-Lobatto quadrature rule, the discretization results in a diagonal mass matrix M. By construction, the SEM discretization allows for the flexibility of having locally (element-wise) defined material coefficient (ρ, λ, G) over the domain and also localized stresses element-wise.se2wave wave propagation functionality has been previously applied in Yuan et al. (2021).
We use an explicit Newmark method as the time integration scheme, a conventional choice for wave propagation problems in SEM (Komatitsch & Tromp, 1999, 2002;Peter et al., 2011), which allows the direct solution of a system of second-order differential equations.Within the Newmark family, we adopt the explicit central differences rule scheme (Hughes, 2000).The computation of the internal forces and the application of the dynamic fault constraints in the procedure are further described below.
For the calculation of the internal forces step, we compute the stress tensor at each quadrature point by using the discrete version of where y is an arbitrary vector.Using Voigt notation, the divergence of stress shown in Equation 12 is given by where   = (, , )  is the Voigt representation of the stress tensor   .Similarly, the strain is described as   = (, , 2)  , which we calculate from a displacement field as ɛ = Bu.We then relate both stress and strain vectors under the linear, isotropic relation in the same notation as We damp spurious oscillations generated along fault by using viscous Kelvin-Voigt damping.For this, we add a term  ̇ to the calculation of the elastic stress field following Day and Ely (2002), and thus apply viscous behavior to both volumetric and deviatoric deformations.There, the viscous relaxation time η = 0.3 Δt, and Δt is the simulation time step, inspired after the Kelvin-Voigt damping parameters in Galvez et al. (2014).Without Kelvin-Voigt damping, spurious oscillations arise in the velocity field, as shown in Figure G1.It will be useful to develop a deeper understanding of the stability of our method in future work, for example, on the basis of a semi-discrete energy balance (e.g., Kozdon et al., 2013).
To implement our stress glut extension, we modify the Voigt stress vector within Σ according to a friction law under a yield criterion.In our case, we use the stress components in a fault-local orientation σ n and τ to evaluate the yield criterion.
The stress modification is summarized in Algorithm 1.The slip and slip rate are updated in accordance with the displacement and velocity fields derived from the modified stress field.Our approach does not require a nonlinear or iterative solve in each time step.Alternative friction laws may require additional steps to update their dependent variables (such as the state variable following Kaneko et al. ( 2008)).

Numerical Discretization
Numerical modeling for seismic wave propagation typically acts as a low-pass filter, accurately propagating low frequencies through the mesh, whilst high frequencies undergo undesired alteration due to numerical dispersion and dissipation (Marfurt, 1984;Seriani & Priolo, 1994).The upper limit of the resolved frequency, conventionally called f max (Hanks, 1982), can be quantified in terms of a number of grid points or elements per shortest wavelength.In the context of SEM, the number of nodes per minimum wavelength follows where p is the polynomial degree to represent the basis functions within a   element of size h.We use these parameters to define the spatial resolution of our SEM simulations.The minimum wavelength is defined as For all simulations shown here (unless otherwise stated), we use  3 elements.Each element contains 4 × 4 G-Legendre-Lobatto integration points with an average spacing of Δx = h/3.

Kinematic and Dynamic Rupture Earthquake Modeling
We have introduced our steady-state phase-field stress glut method as a diffuse interface approach.Here, we apply this approach to earthquake modeling.We explore two well-defined problems: Kostrov's kinematically driven self-similar crack (Kostrov, 1964) and the spontaneous dynamic rupture Southern California Earthquake Center (SCEC) benchmark TPV3 (Harris et al., 2009).Both these problems consider in-plane, mode II rupture propagation.We use SEM2DPACK (Ampuero, 2012) to provide discrete fault reference solutions and compare results by evaluating time series of slip and slip rate at specific points along the fault for SEM2DPACK and our se2dr implementation.We set the reference solutions with elements of polynomial order 6 and a cell refinement h = 100 m for both kinematic and dynamic models.The phase-field smoothing parameters are set to A = 12/δ, φ c = 0.65 δ based on manual calibration.

Kinematic Self-Similar Kostrov Crack
In the following, we vary fault geometry by first considering a straight, mesh-aligned but diffuse fault.Then, we rotate this diffuse fault to not align with the computational mesh, expecting comparable results in important aspects.Finally, we perturb the straight diffuse fault geometry to achieve a curved, sigmoidal geometry.This last model configuration deviates from the reference benchmarks and solutions but provides important information on the geometrical flexibility of our approach.Kostrov's non-singular self-similar shear crack is driven by a time-weakening friction law where rupture evolves under a prescribed constant rupture propagation velocity V r and L is the model characteristic distance.The friction coefficient decreases from a static friction coefficient μ s to a dynamic friction coefficient μ d .This model assumes that the rupture starts from the origin and propagates self-similarly along the fault defined as the arc length integral  Γ = ∫ Γ √ 1 + ( ()∕) 2  as a measure of the accumulated length along the prescribed zero level set geometry x f , parameterized by the variable a.Our model assumes a homogeneous isotropic elastic medium and a predefined fault interface loaded by background normal and shear tractions as defined by Madariaga et al. (1998).The setup allows for analysis of the phase-field relations between fault slip, slip rate, and shear stress under imposed reactions, which avoids the full complexity of spontaneous rupture dynamics.The model parameters are summarized in Table 1.We solve our problem in a domain that spans a 20 × 20 km area, with a fault length spanning throughout the domain.In this first model, we demonstrate that a phase-field simulation can resemble a discrete fault solution in difference to previous findings analyzing thick fault or stress glut fault approaches (Dalguer & Day, 2006;Preuss et al., 2019).Figure 2 summarizes the setup and results of a horizontal Kostrov-like kinematic crack simulation performed with squared  3 and cell size h = 25 m, using the volumetric yielding criterion (see Section 2.1), and fault zone half-width δ = h (see Figure 2a), plotted alongside the discrete SEM split-node reference solution.The phase-field solutions are computed in the diffuse interface model using Equation 2.
We observe close agreement between the diffuse interface and reference models in the time series of slip and slip rate, shown for 5 receiver-pairs located along the fault zone.Phase-field fault slip appears slightly smeared out at its onset and asymptotically very slightly underestimates the classical Kostrov crack solution (Figure 2b).In the diffuse model, the slip rate peak is also slightly delayed and lower in amplitude with respect to the discrete fault reference.Analogous to the reference, slip rates asymptotically fall off after the rupture front has passed (Figure 2c).The snapshot of particle displacement at 4 s simulation time (Figure 2d) illustrates the smooth, well-resolved solution everywhere in our domain.The corresponding velocity and shear stress fields are equally well resolved (Figures 2e and 2f).The zoom-in to the fault zone reveals no out-of-plane rotation of the rupture tip.In general, the phase-field model does not introduce dynamic differences on the scale of the diffuse fault, in difference to what was reported in alternative diffuse interface simulations of the same benchmark (cf. Figure 2 in Gabriel et al., 2021).Changing the yielding criterion (Equation 4 or 6) will lead to only minor differences.The results using the volumetric yielding are smoother in comparison to the diffuse interface yielding criterion as shown in Appendix C, Figure C1a.
In our next example, we first demonstrate the mesh independence of our method.Second, we show that the increased demands on the accuracy of mesh-independent simulations can be addressed by using more elements to resolve the fault zone.We rotate the phase-field and stress tensors that constitute the fault geometry and initial conditions by 20° counter-clockwise from the first Kostrov-crack example.Although the computational mesh is not aligned with the fault, the stress background conditions and model assumptions continue to be the same as in the horizontal configuration.For our tilted configuration, we use the volumetric yielding criterion and a fault zone consisting of a total of five elements, δ = 2.5h.Again, we use  3 -elements, and an element size h = 25 m.
Figure 3 shows slip (a) and slip rate (b) time series recorded along the fault and the x-components of the displacement (c) and velocity fields (d) in the domain.We illustrate that the stress glut phase-field model captures the kinematics, that is, the fault slip (a) and slip rate (b), of the now mesh-independently evolving self-similar Kostrov crack.The slip and slip rate amplitudes are slightly reduced compared to the split-node reference solution.The slip rate time series shows secondary complexities developing within the fault zone after the main rupture front has passed (also visible in d) that do not appear in the reference solution.The emanated seismic waves in terms In our first example, the diffuse fault was perfectly aligned with the element edges.Our smooth phase-field function defined in Equation 7was orthogonal to the element edges and reproduced the split-node reference solution using only two high-order elements, δ = h.The tilted model using this minimal fault zone half-with δ = h (low opacity lines in Figures 3a and 3b), however, produces significantly reduced slip rate amplitude.The slip rate profiles do not show the correct asymptotic behavior after the peak slip rate compared to the reference model and as was observed in the mesh-aligned configuration for the same fault zone resolution.We show in Figures 3a  and 3b that the additional challenges of resolving crack propagation now not orthogonal to the element edges require higher accuracy, which can be achieved by using more elements to resolve our stress glut phase-field fault zone.Earlier smeared crack models (e.g., Rots & Blaauwendraad, 1989) have also considered resolving the crack thicknesses with more than 1-2 elements in their models.However, the stress glut approach has been restricted to using δ = h in earlier work.
Increasing δ/h for a given polynomial order and thus increasing the number of elements that describe the fault zone inelastic rupture kinematics in the case of the tilted Kostrov Model leads to the expected asymptotic behavior.Figure D2 shows h-refinement while keeping δ/h fixed to 2.5.Choosing a larger δ/h leads to better convergence of the numerical solution.In the case of dynamic rupture in the TPV3 model, finding an appropriate δ/h for a given polynomial degree is more complex.In our TPV3 simulations, the width of the nucleation patch is fixed to equal the thickness of the fault zone, which challenges convergence analysis due to the sensitivity of nucleation of spontaneous dynamic rupture to the size (and shape) of the nucleation patch (e.g., Gabriel et al., 2012;Galis et al., 2015).An accurate representation of deformation within the fault zone is governed by the interplay of the δ/h ratio, the polynomial order of the elements, and their alignment with the grid.Furthermore, the ratio between δ and the cohesive zone size characterizes the accurate representation of the deformation at the rupture tip stress transition (see Appendix E).These factors are useful to characterize resolution requirements that lead to accurate fault zone modeling.Additionally, the nucleation zone size should be carefully chosen due to the rupture sensitivity to it.
In Figure C1b, we show the results of the same model using the interface yielding criterion.In comparison to using a volumetric yielding criterion, we see small-scale deviations from the reference slip-rate time series.As we will discuss in the following dynamic rupture examples, these may result from physical fault zone effects.
We conclude that for non-mesh-aligned phase-field models, more elements resolving the diffuse fault zone and using the volumetric yielding criterion are beneficial for quantitatively resembling discrete kinematic rupture propagation.
To further evaluate the geometrical flexibility of the method, we distort the planar Kostrov crack into a sigmoidal curve.The zero level set is parameterized as with parameters k ∈ (−1, 0) and   ∈ ℝ , which control the curvature and the scale of the function, respectively.We make the particular choice to set k = −2 × 10 −4 and A s = 2, which results in the sigmoidal fault configuration shown in Figure 4.In our model, such a curve is prescribed as a discrete set of 4 × 10 5 points.By performing a nearest neighbor search, we identify the closest point on the curve to the quadrature points and use it to initialize the phase-field throughout the domain.In future developments, this can be replaced by, for example, a fast marching method approach to enable a re-initialization at every time step.Such flexibility is advantageous in the context of studies involving time-dependent evolution of fault geometries, where it would be necessary to recalculate the SDF at every time step.
Figure 4 shows the result for the sigmoid configuration using the volumetric yielding criterion, with blending parameters held equal to the tilted configuration and a fault thickness of δ = 2.5h.Its results lead to slip and slip rate profiles well comparable to the discrete fault reference solution, slightly reduced in amplitude, similar to the tilted configuration.Again, our slip rate shows small oscillations behind the rupture front.Note that here, the kinematic model defines the background stress components in fault local coordinates, which implies a spatially heterogeneous background stress for a curved geometry.For this reason, metrics based on the sampling of the near field of a kinematic model are comparable to the metrics obtained from a planar simulation in Figure 3.Further away from the fault, the shear stress wavefield mapped on fault local coordinates shows larger regions of lowered differential stress located at the convex side of the curved fault.

Spontaneous Dynamic Rupture
We model dynamic earthquake rupture in the 2D version of the SCEC/USGS community benchmark problem TPV3 for elastic spontaneous rupture propagation (Harris et al., 2018).Our TPV3 configuration extends the kinematic Kostrov models to spontaneous dynamic rupture propagation.This model uses a linear slip weakening friction law (Ida, 1972) given by where D c is the critical slip distance, and S is the slip, which we extract from the displacement field, as in Equation 2. The model contains a sharp overstressed nucleation patch that initiates self-sustained dynamic rupture.The patch is defined by a length of 3 km, and it is located within a fault of 30 km length as defined in Harris et al. (2009), within a 60 × 60 km domain, and the conditions depicted in Table 2.We compare this model first in the mesh-aligned configuration in Figure 5.We find that the results are quantitatively comparable with the discrete fault reference solution of the TPV3 benchmark calculated with SEM2D-PACK.In this example, we use the interface-yielding criterion.At around t = 0.6 s, the rupture front leaves the overstressed nucleation region, propagating spontaneously along the planar fault.We observe in our results a slight delay in rupture speed compared to the reference solution by comparing the arrival times of the peak slip rate at the receivers along the fault.Comparable to the reference, the fault slip rate approaches an asymptotic falloff behavior after the rupture front has passed.The arrival of the stopping phases when the propagating rupture front reaches the fault edges is observed as an abrupt reduction of the slip rate magnitude after 6 s of simulation time near the end of the profiles.Note that the fault-limiting edges are located well within the simulation domain, far from the limiting boundaries.The model domain is chosen large enough to avoid wave reflections from the domain boundaries.
The mesh-aligned dynamic rupture solution using the volumetric yielding criterion shows a small secondary pulse in the slip rate profiles that originates at the transitional interface between the nucleation patch and the remainder of the fault zone.We find that the fault zone shear stress distribution at the sharp edge of the nucleation patch defined in the benchmark causes fault-oblique yielding across the full fault zone width.The resulting stress shadowing temporarily counteracts local yielding in the phase-field model, while a non-disturbed single spontaneous rupture front develops in the discrete fault reference.Later, with the continuous development of the stress field through time, this location also eventually reaches the yield surface, generating delayed reactivation at the hypocenter, causing a secondary slip rate pulse.Delayed, co-seismic fault reactivation has been reported in real earthquakes, such as during the 2019 northern Peru intraslab earthquake (Vallée et al., 2023), the 2011 Tohoku earthquake (Lee et al., 2011) and the 1984 Morgan Hill earthquake (Beroza & Spudich, 1988).Fault reactivation in discrete fault dynamic rupture simulations can be caused by several model complexities, including pulse-like rupture growing stresses after the passage of its healing front (Gabriel et al., 2012;Nielsen & Madariaga, 2003) and the presence of a fault damage zone, approximated as a low rigidity layer surrounding a discrete fault (Huang et al., 2014;Idini & Ampuero, 2020).
Figure 6a shows slip rates using the volumetric yielding criterion.In comparison, introducing the interface yielding produces a solution closer to the discontinuous reference, free of secondary slip pulses, as seen in Figure 6b.Our analyzed metrics of interest, the slip and slip rate, are extracted from the difference between the displacement field components along the fault parallel direction at ±δ.As a result, asymmetries in the rupture front may introduce a brief fluctuation in the slip rate metric before the main slip rate peak arrival.We note that this fluctuation can result in negative values of slip rate when we use the "interface yielding criterion" (e.g., in Figures 5c and 6b).A comparison of these results against solutions using both yielding criteria is given in Appendix C.
As before, for the Kostrov-like crack, we rotate the dynamic rupture model into a configuration that is out of alignment with the computational mesh (Figure 7).Our experiments show that our choice of fault thickness affects the rupture initiation process in our adaptation of the TPV3 benchmark with a fixed overstress as a direct consequence of setting the nucleation zone width equal to the variable width of the fault zone (i.e., 2δ) within our fault representation.Figure 7 shows a numerical example that uses a fault half-thickness parameter δ = 1.43h, which leads to the qualitatively expected behavior of the rupture.This half-thickness is chosen based on the diagonal length of a square element to ensure that a whole element falls within half of the inelastic zone.When δ = 4.0h, we observe the development of small-amplitude slip pulses in the form of reverberating fault zone waves within the nucleation patch.Later, after the rupture has successfully initiated, the velocity wavefield follows the expected overall behavior in line with the mesh-aligned configuration.At small values of δ = 1.0h, the nucleation size is smaller, and we observe complete dissipation of the rupture front over time, leading to dying (unsustained) rupture.For higher values of δ, the fault zone half-thickness relative to the element width (e.g., δ = 4.0h), trapped waves develop within the fault zone.As detailed in Appendix A, our approach results in an effective modification of the stiffness tensor, leading to a transversely isotropic material within the fault zone.This leads to a locally modified shear modulus relative to the rest of the domain.Constructively interfering fault zone waves later form a coherent rupture front, producing a complex wavefield in the interior of the fault zone and exciting high-frequency seismic radiation visible in the velocity field in the vicinity of the nucleation patch in the fault normal direction (Figure 7, panels with δ = 4h).Since we do not allow material failure outside the  prescribed finite thickness fault zone, which is stationary in time, newly yielding localization of secondary faults, continuously growing away from the inelastic subdomain of the main fault, are not expected as a contributing source to high frequency in the model.
The effects of different choices of δ in our non-mesh-aligned numerical tests agree with previous findings (e.g., Gabriel et al., 2012;Galis et al., 2015;Huang et al., 2014): slight variations of the nucleation size, for a fixed prescribed overstress, can lead to unsuccessful initiation of the rupture process on the lower end of the parameter space allowing for self-sustained rupture, implying that the initiation is not sufficiently strong for supporting rupture to spontaneously propagate and develop into a self-sustained propagating rupture.In the overcritical limit, larger patches introduce changes in rupture dynamics, including changes in rupture style and speed, such as super-shear transitions and higher slip-rate amplitudes.

Spectral Properties of the Modeled Seismic Wavefield
For the element choice in our mesh, the shear wave velocity assumed for the TPV3 model, and assuming several integration points per minimum frequency as ζ min = 30 due to the relatively low polynomial order used, we can assess the upper cut-off frequency as approximately 13 Hz.We chose this value for ζ min from the suggested range of 10-30 for low (<4) element order after Seriani and Priolo (1994).For the mesh settings of the reference solution and an assumed ζ min = 10, the reference solution has an upper cut-off frequency of 20 Hz.With this information, we limit the upper band of the frequency spectra in Figure 8d to the cut-off frequency.
Rupture acceleration and deceleration generate high-frequency seismic wave radiation (Madariaga, 1977).For this reason, we expect a roughly flat signal in the spectra of the acceleration records for the kinematic model, where the rupture velocity is constant.The amplitude spectrum from receivers in the Kostrov kinematic model solution is shown in Appendix B. For the dynamic model, we extract the along-fault (a) and fault-normal (b) accelerograms in Figure 8 from receivers at different distances normal to the fault.The accelerogram spectral amplitudes of the dynamic model in Figure 8d at two receiver locations; both 6 km along the fault from the fault center, and respectively at a distance of 25 and 500 m normal to the fault, are increasing until just above 1 Hz, with amplitudes systematically higher than the reference solution.
Discrete coseismic off-fault damage has been considered to enhance high-frequency radiation in acoustic recordings during stick-slip events (Castro & Ben-Zion, 2013;Hanks, 1982).K. Okubo et al. (2019) finds significant high-frequency radiation caused by secondary discrete fractures in simulations compared to the no-off-fault damage case.In our diffuse simulation using an interface-yielding criterion, we overall observe a similar trend as in the reference solution, with no significant shift toward the high frequencies.Analogous to the reference solution, the frequency content decays rapidly with the fault-normal distance, roughly reduced by one order of magnitude near the upper cut-off frequency.However, the fault-reactivation slip pulse (observed in Figure 5, using the volumetric yielding criterion) contributes to the higher frequency content, shifting the amplitude spectra upwards toward the upper cut-off frequency in comparison to the reference solution.This high-frequency contribution can be seen in Figure C3.

Resolution Refinement Analysis
A formal convergence analysis requires an analytical solution that is not available for our steady-state phase-field diffuse fault approach.Instead, we present refinement analysis by means of comparison to a high-resolution reference solution (Pelties et al., 2012;Wollherr et al., 2018) and illustrate convergence toward the reference solution under h-and p-refinement and variation of blending parameters.On the right column, the subset illustrates two end-members of the δ parameter variation: on the top, the dissipation of the rupture front, and below, the formation of small amplitude resonance in the velocity field within the nucleation zone.In our numerical simulations, we link the size of the nucleation zone to the corresponding thickness of the fault zone.As a consequence, the behavior of the numerical simulation is sensitive to the chosen fault width 2δ.
We here choose a steady-state phase-field description of the diffuse transition from the yielded, inelastic region into the pure elastic media, which offers a flexible numerical approach for the mesh-independent representation of a discontinuity.The selection of the blending parameters influences the accuracy of our metrics of interest, the slip and slip rate, against the high-resolution reference solution.Steeper transitions lead to spurious oscillations behind the rupture front, a product of Gibbs phenomena, with strong signal amplitudes, while smoother transitions lead to reduced signal amplitude.Our choice, although not optimized, is intended to balance the trade-offs of end-member choices in the blending parameter space.Future work may explore physic-based considerations to inform these choices.
The blending variables may be interpreted as additionally introduced degrees of freedom to fit against a reference solution; alternatively, they can be further constrained under phase-field theory and be considered a proxy for material damage.Variation of the blending parameters influences the accuracy against a reference solution.Figure 9 (top row) shows the solution to the TPV3 benchmark problem by using a set of blending parameters A = 18/δ, φ c = 0.65δ, that is, a steeper blending than in the result from Figure 5, also using the interface yielding criterion.
The results approach the reference solution, reflected in reduced differences in peak slip rate amplitude and timing.For this choice of blending parameters, we perform mesh refinement analysis and a variation of the simulation polynomial degree for the dynamic rupture TPV3 benchmark, which we compare against the well-defined high-resolution reference solution to the benchmark problem.We use a Butterworth filter with a cut-off frequency of 10 Hz for solutions with different mesh size and order of polynomial refinements (h-and p-refinements respectively) and compare how the peak slip rate of our results differ from the reference in terms of amplitude and timing, as seen in Figure 9 (bottom row).We evaluate the differences for the receivers located at 4, 6, and 8 km along the fault.
Our results show systematically lower amplitudes for our  3− element solutions.However, we also report on results when further increasing the resolution in the simulations.p-refinement leads to faster growth of the peak slip rate amplitude toward the reference and reduced timing differences for the different receivers along the fault.The systematic delay in timing gets reduced by cell refinement, which also affects the inelastic fault zone width.The amplitude of the peak slip rate depends on the discretization within the fault zone and the blending parameters used, as it describes the offset from the elastic stress response and its transition into the elastic media at the quadrature points within the subdomain.Also, using  2− element fails to reproduce the asymptotic fall-off behind the peak slip rate arrival (Figure D7), denoting a requirement for higher spatial resolution.
Our method reproduces the reference solution at a relatively low polynomial refinement for a given phase-field distribution choice.Its accuracy is affected by the resolution of the yielding elements, implying that adaptive 10.1029/2023JB027143 19 of 40 mesh refinement can be applied to the method for future optimization purposes.Higher polynomial orders can be tested after optimizing our implementation to establish that the method maintains the same numerical accuracy (convergence order) as the classical SEM.

Discussion
Dalguer and Day (2006) assessed the accuracy of the Traction at Split Node (TSN) method, the thick fault proposed by Madariaga et al. (1998), and the stress glut from Andrews (1976) in a staggered grid finite difference discretization.The explicit incorporation of discontinuous velocity nodes in the TSN method allows for a natural partition of the equations of motion at each side of the fault surface on which the split nodes are located.
In this context, the stress glut and the thick fault methods require a fixed fault thickness with respect to the computational grid resolution.The stress glut method's fault zone width is two halves of contiguous elements, that is, an inelastic zone of one grid step, with the fault center defined as the common border defining a plane of shear-stress grid points.The thick fault places the fault halfway between two shear-stress planes, and the fault zone thickness is two grid steps.Under these assumptions, the stress glut method reproduces well the qualitative features of the discrete fault reference solution but with quantitative deficiencies.It was reported that the rupture velocity remained systematically low.The thick fault reported a misfit in rupture time of 30% and higher, which then failed to match the rupture behavior of the reference solution qualitatively.One of the main conclusions of Dalguer and Day (2006) was that the formulation of the jump condition mainly controls the solution accuracy in finite-difference spontaneous rupture simulations.
In this work, we have modified the stress glut approach and have improved the solution accuracy by using a steady-state phase-field model, enabling a smooth transition between the yield stress and the elastic shear stress.
Our stress glut extension to the framework of SEM in the mesh-aligned configuration shows overall qualitative and quantitative agreement with 2D benchmarks of kinematic and dynamic rupture problems when verified against a split-node SEM reference from Ampuero (2012).In general, the solutions show an expected systematic delay of the rupture front arrival (i.e., slower rupture speed), depending on the prescribed half-thickness of the fault zone δ, relative to the element dimensions.Introducing the diffuse interface description reduces fault zone spurious oscillations introduced by Gibbs phenomena due to the stress discontinuity from the imposed limiter on the stress.Similar to the typical employed visco-elastic Kelvin-Voigt damping in split-node SEM dynamic rupture modeling, and equivalent to introducing off-fault plasticity or damage (Andrews, 2005;Day et al., 2005;Gabriel et al., 2013;Xu et al., 2015), our diffuse fault zone introduces reduction of amplitude in both slip and slip rate metrics as well as in rupture speed.A higher p-refinement level combined with h-refinement and our proposed blending function (e.g., Figure 5) approach the reference solution in the mesh-aligned case, with reduced spurious oscillations behind the rupture front.

Physical Interpretation of the Stress Field of a Diffuse Fault
Andrews (1999) pointed out that embedding a crack through the stress glut formulation affects the neighboring stress in an irregular way that can be compared to the Eshelby inclusion problem (Eshelby & Peierls, 1957).The complexity of the stress field incurred by such inelastic "inclusions" increases when it interacts with the evolving stress field around the dynamically propagating rupture and may prevent locations within the fault zone from reaching the yielding surface at a specific time.
This increased complexity in the stress field directly affects our dynamic rupture results, as before the onset of yielding and development of the fully spontaneous rupture front, distributed shear stress locally shadows the fault zone regions at the vicinity of the transition between the nucleation patch and the rest of the fault zone.The incipient rupture develops asymmetrical dynamic normal stress evolution, which leads to fault-oblique yielding within the fault zone.Note that this oblique geometry characterizes shear-driven deformation between two surfaces undergoing relative motion at each side of the embedded fault zone.The fault geometries used in this study are prescribed and do not evolve in time, hence, by construction, our models do not permit the development of spontaneous Riedel-type shear structures.However, the observed emergence of oblique yielding within the fault surface may be interpreted as an evolving fault-zone internal shear band.This phenomenon at the vicinity of the boundary between the prescribed nucleation patch and the remainder of the prescribed fault zone alters the timing of the onset of rupture at locations neighboring the nucleation patch when using a volumetric yielding criterion in mesh-aligned numerical simulations.These temporally "dynamically locked" patches are later reactivated by evolving stresses in their vicinity, producing a measurable, small amplitude propagating secondary slip-pulse.We note that the TPV3 model setup includes a challenging characteristic, as it contains a sharp transition between the nucleation asperity and the rest of the fault (Galis et al., 2015); however, we also observe dynamically impacted fault zone yielding in alternative descriptions of the nucleation patch (see in Appendix C, Figure C2).
We designed the interface yielding criterion as defined in Equation 6 to suppress the observed tendency for fault-oblique yielding at asperities within the fault zone.As further discussed in Appendix C, this alternative yielding criterion implies the yielding of all regions behind a fully spontaneous fault zone rupture front.It is a justifiable assumption applied to our diffuse fault when our goal is merely to emulate the results from planar, discrete fault representations.However, for thicker fault zones, the interface yielding assumption introduces perturbations on the stress field within the fault zone, especially at the vicinity of the rupture front, which in turn introduce spurious oscillations of the shear stress component propagating in fault-strike direction (as we show in kinematic and dynamic rupture simulations).This simplification can also be perceived as a variation of Hooke's law that integrates rotations along with strains in evaluating the yield criterion without affecting the elasticity update.
Distinct features are evident throughout the fault-parallel internal deformation within the inelastic zone, as shown in Figure 10.These features include slip-strengthening behavior, double slip-weakening, and nonlinear weakening behavior with long tails.Such behaviors have been assumed to reflect true frictional behavior (McLaskey et al., 2015;Ohnaka & Kuwahara, 1990;Ohnaka & Yamashita, 1989;P. G. Okubo & Dieterich, 1981;Rubino et al., 2017).However, Xu et al. (2019) argued that they might instead capture rupture behavior in off-fault locations.Hence, the latter perspective supports the development of indirect approaches to estimate rupture properties.

Mesh Independence
Our steady-state phase-field approach does not require the mesh to be designed to align with a pre-existing fault.We show that kinematic and dynamic rupture can evolve independently of the spectral element boundaries.
Results obtained with the mesh-aligned relative to the dynamic rupture problem lead to a close match with the reference solution.Non-aligned mesh configurations using  3 elements show a general amplitude reduction in the slip and slip rate metrics and grow closer to the reference solution-as the mesh-aligned case-when increasing the fault zone width.Alternatively, the non-aligned mesh solutions require an increased spatial resolution of the diffuse fault zone to reach the same level of agreement with the fault-interface reference solution as in the mesh-aligned case.This effect becomes especially apparent for the spontaneous dynamic rupture models, where a low integration point density of the elements constituting the nucleation patch may prevent self-sustained spontaneously propagating rupture due to the sensitivity of dynamic rupture to nucleation size, shape and procedure (e.g., Bizzarri, 2010;Festa & Vilotte, 2006;Gabriel et al., 2012Gabriel et al., , 2013;;Galis et al., 2015;Lu et al., 2009;Shi et al., 2008).Alternative nucleation strategies that are smooth in space and/or time can reduce numerical artifacts in spontaneous dynamic rupture problems, such as stress localization in the surroundings of a sharply defined nucleation patch.Such smooth approaches also minimize the influence of a potentially ill-constrained nucleation procedure on the subsequent stages of realistic earthquake scenario simulations (e.g., Biemiller et al., 2022;Harris et al., 2021).We observe comparable numerical artifacts, apparent especially in the stress fields, at quadrature points of element cells located only partially within the diffuse fault zone, for example, in the mesh-independent fault configurations using low polynomial order (see Figures 4,7,and C3).Increasing fault zone width can mitigate this issue.However, at the same time, too high nucleation overstresses should be avoided.

Alternative Smeared Crack Models
Smeared crack models have been applied within the framework of finite element methods to model the fracture mechanics of mode I, II, and mixed-mode cracks in concrete.Thereby, the so-called "stress locking" (Rots, 1988;Rots & Blaauwendraad, 1989) phenomenon refers to spurious stress build-up around the cracking elements, which may pollute the numerical results and lead to an overestimated energy dissipation and non-zero residual strength of a cracked structure.The cause of this spurious stress transfer has been related to a poor kinematic representation of the discontinuous displacement field in the vicinity of the macroscopic crack (Jirásek & Zimmermann, 1998a, 1998b).Unless the direction of the macroscopic crack (represented by a band of cracking HAYEK ET AL.

10.1029/2023JB027143
21 of 40 elements) is parallel to the edges of finite elements, the directions of maximum principal strain determined from the finite element interpolation at individual integration points of the element deviate from the normal to the crack band.The principal lateral stress has a non-zero projection on the crack-band normal, generating spurious cohesive forces.However, higher-order elements offer better kinematic flexibility, such that they can partially relax the spurious stresses by adjusting the interpolated displacement field.Jirásek andZimmermann (2001a, 2001b) deals with stress locking and instability via a varying scalar damage parameter as a function of the crack deformation of a mixed-mode embedded crack.
The diffuse thick fault approach presented here keeps the crack width fixed throughout the simulation, which may perturb the stress component at elements crossing the fault interface.This argues in favor of considering non-local approaches in the future or extending our approach into a non-steady-state phase-field method.Non-local gradient models enrich the local constitutive relations with first-order or higher-order gradients of a state variable or according to associated thermodynamic forces (e.g., Marotti de Sciarra, 2009).Non-local integral models involve weighted averages of a state variable around that point (e.g., Lyakhovsky et al., 2011).Such models may rely on an invariant-based formulation, as it is also used in inelastic yielding models (e.g., Templeton & Rice, 2008), to evaluate and update an energy function to condition the evolution of a damage parameter.Peerlings et al. (1996) describe the responses obtained from non-local and gradient damage approaches as qualitatively similar.Gradient-dependent formulations include all non-local model features essential for describing localization phenomena, with non-negligible quantitative differences arising from the absence of high-order derivatives in the gradient formulation.We may also want to choose a functional based on a non-local plasticity formulation describing the transition from the inelastic to the pure elastic domain within our fault-local compact support as an alternative to our blending approach to reinforce the diffuse character of the method and avoid localization.This treatment may be approached by defining the yielding condition and metrics in an average sense, which offers an alternative treatment to the stress localization described in Section 5.1.Future work may explore the physical constraints to inform the blending parameters that mathematically define our diffuse fault model.
Recently, Fei and Choo (2019, 2020, 2021) have described phase-field models of geological motivated rock fracture that incorporate pressure dependence and frictional contacts for mode I, II, and mixed crack modes.Their approach is formulated as a set of governing equations for different contact conditions in the FEM framework where frictional energy dissipation emerges in the crack driving force during slip.Their method is proposed to allow arbitrary interface geometry representation without an explicit function or enriched basis, an advantage of phase-field methods.It can also accommodate contact constraints without a dedicated algorithm.Their approach ensures that the nonslip direction's stress tensor component is compatible between the interface and bulk regions.This results in modifying the separation between the volumetric-deviatoric stress decomposition approach proposed by Amor et al. (2009).Our formulation of the modified stress tensor in Equation 5resembles theirs for a prescribed shape of stationary crack interface for the phase evolution equation.Our blending function and SDF variables play an equivalent role to the degradation function on their phase-field variable.After modifying the stress in their approach, they use Newton's method to solve the discretized momentum balance equation, and the nodal increment in displacement requires explicitly solving for an updated stiffness tensor to calculate the element-wise Jacobian at the end of each time step.In this context, the crack driving force is calculated from the change of the plastic strain and used to calculate the updated phase-field variable.
A diffuse description motivated by steady-state phase-field profiles allows us to explore the yielding surface's transition into elastic media as a distribution across the fault representation while keeping its logical simplicity in the formulation.This allows the method to be ported into alternative numerical frameworks.Development of our representation into the framework of phase-field requires critical ingredients of the phase-field formulation (based on the theory of fractures of Griffith and Taylor (1921)), such as introducing a phase-field evolution equation and the incorporation of a critical fracture energy which translates into the regularized continuum gradient damage mechanics (Miehe et al., 2015).This increases the complexity in its formulation and introduces parameters to solve for.However, the evolution of the phase-field, and thus the fault-normal growth at different distances along the embedded inclusion of the yielding surface, is pertinent to natural observations.Fault lateral growth is observed in nature as changes in the structural fault complexity along the propagation direction of the parent fault (Perrin et al., 2016).Such variation may avoid accumulating localized stress components throughout time at the inter-element boundary within the numerical grid.In addition, such a hypothetical spectral-element-based phase-field method would avoid explicitly calculating an updated stiffness tensor at the end of each time step.
In contrast to our diffuse fault representation with constant blending and respective parameters, alternative diffuse fault models incorporating increased thermomechanical complexities have been developed.A contemporaneous diffuse crack representation incorporating finite strain nonlinear material behavior and multi-physics coupling into dynamic earthquake rupture modeling is described by Gabriel et al. (2021) using the Godunov-Peshkov-Romenski (GPR) model (Resnyansky et al., 2003;Romenskii, 2007).The model uses a first-order hyperbolic model of inelasticity coupled to finite strain elasto-visco-plasticity of continuum mechanics (Tavelli et al., 2020) and is extended for dynamic rupture using a high-order Discontinuous Galerkin scheme and the ExaHype partial differential equation (PDE) engine (Reinarz et al., 2020).Their model also permits the representation of arbitrarily complex geometries via a diffuse interface approach.In neither of their two scalar fields, the local material damage describing the fault geometry and secondary cracks and the solid volume fraction function need to be mesh-aligned, allowing faults and cracks with complex topology and using adaptive Cartesian meshes (AMR).However, the problem of parameter selection for their unified model of continuum mechanics is a non-trivial task due to the large amounts of parameters and may require numerical optimization algorithms applied to data obtained from observations and laboratory experiments.Our method also requires locally high resolution to describe the diffuse fault zone and would benefit from the future implementation of fault zone AMR, building upon previous implementations of SEM with AMR (e.g., Rudi et al., 2015Rudi et al., , 2017;;Tanarro et al., 2020).

Outlook
Here, we explore simple 2D benchmarks for kinematic and spontaneous dynamic earthquake rupture, including geometrical complexities on a structured mesh.The next step can involve exploring the method's potential for modeling branching and crossing faults.In future work, fault junctions and geometrical complexities such as sharp bends may be implemented using hierarchy levels of fault entities with their respectively defined independent fault zone characteristics and updating the stress field in an iterative manner.In that way, for example, handling different thicknesses of fault zones per hierarchy level would be possible.However, stress concentrations associated with sharp bends or junctions may require careful analysis (Andrews, 1989;Uphoff et al., 2022).Fault intersections and dynamic fault interactions alter the spatial distribution of stress concentrations (e.g., Taufiqurrahman et al., 2023) as well as influence the earthquake energy budget (e.g., K. Okubo et al., 2019), thus, directly affect earthquake rupture dynamics.Future extension of our approach to 3D is essential for direct observational constraints and verification studying real earthquakes.3D unit elements are well established in spectral element methods applied to seismic wave propagation (e.g., Komatitsch & Tromp, 1999) and rupture dynamics (Galvez et al., 2014) and Andrews (1999) demonstrates the feasibility of strategies to evaluate the slip and slip rate from the shear traction components in 3D calculations.Thus, we expect that modifying our approach via a diffuse description of the stress glut should readily be extendable to 3D.Applications of our method may help to further our understanding of fault zone evolution and the effects of internal rheology distribution at coseismic time scales.SEM is a volumetric method that allows for variable material parameters and mesh independence.In principle, this will allow our method to model time-evolving fault geometries, for example, to capture the coupling between different physical mechanisms on-fault and within the bulk and evaluate their relative importance for rupture dynamics.
Thorough additional analysis will be required to extend our approach to rate-and-state friction.While such implementation and analysis are outside the scope of the paper, we envision that the main changes required to incorporate rate-and-state friction within our method are: • Define a state variable at points living along the zero-level set contour (not defined within the volume).
• Change the method to evaluate the friction (between lines 7 and 8 in Algorithm 1).
• Modify the time integrator to use adaptive time-stepping.
• Add the evolution of the ordinary differential equation (ODE) for each state variable within the exiting time loop used to advance the displacement solution.

Conclusion
In this work, we present a novel steady-state phase-field model for rupture dynamics that extends the stress glut approach (Andrews, 1999).Using the high-order accurate and geometrically flexible framework of spectral elements, our diffuse fault zone formulation results in comparable kinematic and dynamic rupture propagation to the conventional planar TSN SEM for dynamic rupture simulations.Our approach supports a general description of the evolution of the effective friction coefficient, which dictates fault yielding and sliding as a function of time, location, slip, slip rate, and potential additional variables.To verify our approach, we first compare mesh-aligned kinematic and dynamic rupture model solutions.Our stress glut spectral element implementation aligns well with the discrete fault split-node spectral element reference solutions.Moving beyond the conventional planar interface, we introduce a diffuse fault zone description.This novel representation condenses fault volumetric complexities into a distribution within a compact support.This diffuse fault description follows a prescribed blending function that characterizes the transition from the inelastic state of the embedded crack to the pure elastic state of the surrounding rock matrix.Our model demonstrates its versatility using mesh-independent planar and curved fault geometries, simplifying the often tedious task of mesh generation.Importantly, our steady-state

Appendix C: On the Yielding Criterion Applied to the Rupture Model
Our setup applied to the kinematic Kostrov crack produces slight differences between the model with the interface yielding criterion and the model with the volumetric yielding criterion.In a mesh-aligned geometrical configuration, the main difference is a smooth overshoot slip rate prior to the peak slip rate arrival, as seen in Figure C1a, for the volumetric criterion, while the interface criterion instead contains a sudden reduction at this position.In the mesh unfitted geometrical configuration (Figure C1a), using the interface criterion with a thick fault geometry introduces numerical oscillations to the trailing signal behind the rupture front.As indicated in Section 5.1, our setup applied to the TPV3 model generates a fault-oblique yielding which, in the transition between the nucleation and the rest of the fault zone, leaves an unyielded location within the fault zone, behind the rupture front.As the simulation progresses, this location reaches the yielding surface due to the evolution of the stress field, producing a small pulse in the slip rate profile when using the volumetric yielding criterion, while introducing the interface yielding produces a solution closer to the reference, free of secondary pulses, as seen in Figure 6.Taking a look into the stress field, the shear stress centered around the transition between the nucleation and the rest of the fault zone for the same model configuration under a reduced element width h = 50 m is depicted in Figure C2a for the volumetric yielding criterion.We extract a transect crossing such a non-yielded location and sample the shear and fault normal stresses, shown in Figure C2b.The same is done for a simulation that uses the interface yielding for Figures C2c and C2d.Note the asymmetry of the shear stress profile in Figure C2b, with positive values of the differential shear stress due to the unyielded location, while Figure C2c shows the state of such stresses, and no unyielded location is left behind the developed fully propagating rupture, which advances past this point.Note that the solution using the volumetric criterion also delivers a smooth stress field within the fault zone outside the nucleation patch.In contrast, the interface yielding solution generates a slight oscillatory perturbation.For this reason, we consider the interface yielding criterion a gateway to emulate planar interface solutions.
The small pulse in the slip rate profile decays fast with a distance normal to the fault as observed in Figures C3a and C3b and contributes to high-frequency signal as observed in the spectrogram in Figure C3c, producing a high-frequency content in the amplitude spectrum toward the cut-off frequency.
The fault internal deformation for the volumetric yielding approach (Figures 10 and C4) depicts how our method handles explicitly the internal deformation in terms of the displacement, velocity and stress components within a finite-thickness zone.A close look into the fault-parallel velocity components in the transects of Figure C4 shows that in the neighborhood of the rupture front, the velocity field can behave in a skewed manner relative to the zero level set for the Q 3 TPV3 simulation in the horizontal configuration.

Figure 1 .
Figure 1.Schematic of the diffuse fault representation using the signed distance function.The mesh independent fault indicator φ(x) is defined within an inelastic zone width of 2δ and acts in the subdomain Σ ⊂ Ω.The blue and yellow circles indicate the projected coordinate pairs x + and x − at a distance δ from the zero level set on opposite sides of the fault, as described in the text of Section 2.1.Each circle includes the fault local orientation axis (t, n).

Figure 2 .
Figure 2. Phase-field stress glut model results for a kinematic Kostrov crack with the mesh-aligned configuration using our diffuse fault zone approach.The model portrays an in-plane right-lateral shear fracture under compression using the volumetric yielding criterion.The structured mesh is composed of square  3 elements with a width of h = 25 m.The fault zone half thickness equals one element width, δ = h.The model evolves for 4 s of simulation time.(a) Illustrates the embedded fault in the mesh and the distribution of the receiver pairs at increasing along-strike distance from the hypocenter indicated by color, where the slip and slip rates are extracted.The next subfigures include our adopted metric of (b) slip and (c) slip rate profiles in continuous lines, compared against an spectral element method (SEM) split node discrete fault reference solution as dashed lines.We also show the corresponding x-component snapshots at t = 4 s for the (d) displacement, (e) velocity, and (f) shear stress.

Figure 3 .
Figure 3. Phase-field stress glut model: Kostrov-like crack model under a tilted, mesh independent geometry.The model uses the same mesh and model parameters as in Figure 2, with a 20° counterclockwise tilting, relative to the mesh axis, and a wider fault zone (δ = 2.5h).The model uses the volumetric yielding criterion.The figure depicts (a) the slip profile and (b) the slip rate profile compared against an spectral element method (SEM) split node discrete fault reference solution at receiver pairs with increasing along-strike distance from the hypocenter indicated by color.Additionally, the figure contains the x-component of the (c) displacement and the (d) velocity field with an inlet focused on the propagating front.In (a) and (b), we highlight the effect of choosing fewer elements to resolve the fault zone width δ by plotting the slip and slip rate results of the same tilted model using δ = h in lower opacity.

Figure 4 .
Figure 4. Phase-field stress glut model: Kostrov-like crack with a sigmoid shape.The model uses the same mesh parameters and thickness of the diffuse zone (δ = 2.5h) as in Figure 3.This configuration depicts a sigmoid shape with a zero-level set following Equation 18, using the volumetric yielding criterion.The figure contains (a) the slip profile and (b) the slip rate profile compared against an spectral element method (SEM) split node discrete fault reference solution, with receiver pairs at increasing along-strike distance from the hypocenter indicated by color, and (c) the x-component of the velocity field and (d) the shear component of the stress in fault local coordinates.

Figure 5 .
Figure 5. Phase-field stress glut model: Mesh-aligned TPV3.The mesh is composed of squared  3 elements and element width h = 25 m.This model uses the interface-yielding criterion.The fault zone parameter δ = h, and the system evolves for 7 s of simulation time.(a) Illustrates the model configuration, including the location of the nucleation within the fault zone and the receiver pairs at increasing along-strike distance from the hypocenter indicated by color.Next to it, the figure includes the (b) slip and (c) slip rate profiles compared against the reference solution profiles.The figure also includes the result's corresponding x-component snapshot at t = 3 s for the displacement (d), velocity (e), and shear stress (f).

Figure 6 .
Figure 6.Comparison of the 2D TPV3 dynamic rupture solution yielding criteria, with square  3 elements of width h = 25 m, and blending parameters as in Figure 5. (a) Shows the simulation results at receiver pairs with increasing along-strike distance from the hypocenter indicated by color, using the volumetric yielding criterion, while (b) depicts the results using the interface yielding criterion.

Figure 7 .
Figure 7. Phase-field stress glut model: Variation of the δ parameter for tilted TPV3 simulations.The mesh is composed of square  3 elements of width h = 25 m.The figure includes three sets of frames per time step.Each frame set contains the x-component of the velocity field and showcases three values of the half inelastic zone parameter δ; h, 1.43h, and 4h.The models here use the volumetric yielding criterion.Each figure subset depicts on the left column the expected qualitative behavior for the intermediate half inelastic zone thickness.On the right column, the subset illustrates two end-members of the δ parameter variation: on the top, the dissipation of the rupture front, and below, the formation of small amplitude resonance in the velocity field within the nucleation zone.In our numerical simulations, we link the size of the nucleation zone to the corresponding thickness of the fault zone.As a consequence, the behavior of the numerical simulation is sensitive to the chosen fault width 2δ.

Figure 8 .
Figure 8. Phase-field stress glut model, TPV3 mesh-aligned model: Variation of the (a) x-and (b) y-components of synthetic accelerograms, at stations located at 6 km along the fault, and varying distances normal to the fault for the simulation in Figure 5. (c) Spectrogram extracted from the y-component of the acceleration from a receiver at the coordinates (6 km, 25 m).(d)Amplitude spectra of the fault-normal accelerograms at two receivers at 6 km along the fault and, 25 and 500 m normal to the fault, simulated with se2dr (continuous and dotted lines for each yielding criterion used) and the split-node discrete fault approach in SEM2DPACK (dashed lines).All spectra are shown to the numerically resolved f max = 13 Hz.

Figure 9 .
Figure 9. Filtered (Butterworth with f c = 10 Hz) slip rate profiles.The top row shows the simulation slip rate results to the TPV3 mesh-aligned model with alternative blending parameters A = 18/δ, φ c = 0.65δ.These results use  3 elements with h = 25 m and δ = h and an interface-yielding criterion color-coded by station location along the fault.The second row depicts the difference between the peak slip rate and timing between the reference and our simulation results from filtered slip rate profiles from three receivers at 4, 6, and 8 km along the fault.The scatter plots include such differences for the blended

Figure 10 .
Figure 10.Internal deformation of the horizontal TPV3 model using Q 4 square elements of width 25 m, δ = 12.5 m, and the volumetric yielding criterion, embedded in a 20 × 20 km domain size.(a) Depicts the snapshot of the x-component of the velocity field, with a zoom-in on the rupture tip and an indication of the extracted transects.(b) Contains an extract from Xu et al. (2019) indicating the apparent friction coefficient constructed from experimental data and from numerical synthetics computed by shear-to-normal stress ratio, using different vertical offsets from the fault for different rupture velocities.Fault-parallel transects extracted at the (c, d) 25 m, (e, f) 12.5 m, (g, h) 0 m, (i, j) −12.5 m, and (k, l) −25 m level sets after 1.08 s of simulation time.The transects are equidistantly sampled, amounting to 10,000 points from the center of the domain to 10 km at the end of the domain.Each row axis pair contains first a plot of the x-component of the displacement (solid lines) and velocity (dashed lines), respectively, and second, the shear and normal stress component (color-coded) sampled as a function of distance along the fault per level set of interest.The profiles are only shown from 0 to 7 km distance from the epicenter along the fault.

Figure C1 .
Figure C1.Comparison of the 2D Kostrov crack solution yielding criteria.The figure contains (a) the mesh-aligned setup shown in Figure 2 with δ = h using the volumetric yielding criterion (dotted colored lines) and the interface yielding (continuous colored lines).The same comparison for both yielding criteria is shown for (b) the tilted setup shown in Figure 3 with δ = 2.5h.

Figure C2 .
Figure C2.Comparison of the 2D TPV3 dynamic rupture solution.The mesh is composed of  3 square elements of width h = 50 m.The solutions shown here use the volumetric yielding (top row) and the interface yielding (bottom row).The blending parameters and the ratio between the fault inelastic zone width relative to the element width are the same as in the results from Figure 5. (a) Depicts the shear stress field zoomed at the transition between the nucleation and the rest of the fault zone.Superposed is the location of the transect extracted in (b), sampling the fault normal and shearing components of the stress field.The center of the transect is located at (1,513.8 m, 0 m) so that it crosses the location left unyielded at the time behind the rupture front.Likewise, (c) and (d) show the shear stress field and the transect, respectively.

Figure C3 .
Figure C3.Phase-field stress glut model, TPV3 mesh-aligned model: Variation of the (a) x-and (b) y-components of synthetic accelerograms, at stations located at 2 km along the fault, and varying distances normal to the fault for the simulation in Figure 6 employing the volumetric yielding criterion in horizontal configuration.(c) Spectrogram extracted from the y-component of the acceleration from a receiver at the coordinates (2 km, 25 m).(d) Amplitude spectra of the fault-normal accelerograms at two receivers at 2 km along the fault and, 25 and 500 m normal to the fault, simulated with se2dr (continuous lines) and the split-node discrete fault approach in SEM2DPACK (dashed).

Figure D2 .
Figure D2.h-refinement test for Kostrov's model in the tilted (20°) geometrical configuration.The models depicted use Q 3 elements with element width h = 25, 50, and 100 m, and δ = 2.5h.We impose a volumetric yielding criterion within the fault zone.

Figure D3 .
Figure D3.h-refinement test for Kostrov's model in the sigmoid geometrical configuration.The models depicted use Q 3 elements with element width h = 25, 50, and 100 m, and δ = 2.5h.We impose a volumetric yielding criterion within the fault zone.

Figure D4 .
Figure D4.h-refinement test for TPV3 model in the tilted (20°) geometrical configuration.The models depicted use Q 3 elements with element width h = 25, 50, and 100 m, and δ = 1.43h.We impose a volumetric yielding criterion within the fault zone.

Figure D5 .
Figure D5.Horizontal configuration of TPV3 model, simulation using Q 3 elements with element width of h = 25 m, and δ = h, through 7 s of simulation time.Using volumetric yielding criteria.The profiles include additional receiver locations along the fault geometry, every 2 km from the nucleation up to 14 km.

Figure D6 .
Figure D6.Tilted (20°) configuration of TPV3 model, simulation using Q 3 elements with element width of h = 25, 50, and 100 m, and δ = 1.43h, through 4 s of simulation time.Using volumetric yielding criteria.The profiles include additional receiver locations along the fault geometry, every 2 km from the nucleation up to 10 km.

Figure D7 .
Figure D7.Filtered (Butterworth with f c = 10 Hz) slip rate profiles for results using Q 2 elements and alternative blending parameters A = 18/δ, φ c = 0.65δ, following the parameter choices for the lower left inset of Figure 9.

Table 1 Parameters
Describing Our Kostrov-Like Self-Similarly Propagating Kinematic Shear Crack Model