A comparative analysis of continuum plasticity, viscoplasticity and phase-field models for earthquake sequence modeling

This paper discusses continuum models for simulating earthquake sequences on faults governed by rate-and-state dependent friction. Through detailed numerical analysis of a conventional strike-slip fault, new observations regarding the use of various continuum earthquake models are presented. We update a recently proposed plasticity-based model using a consistently linearized formulation, show its agreement with discrete fault models for fault thicknesses of hundreds of meters, and demonstrate mesh objectivity for slip-related variables. To obtain a fully regularized fault width description with an internal length scale, we study the performance and mesh convergence of a plasticity-based model complemented by a Kelvin viscosity term and the phase-field approach to cohesive fracture. The Kelvin viscoplasticity-based model can introduce an internal length scale and a mesh-objective response. However, on grid sizes down to meters, this only holds for very high Kelvin viscosities that inhibit seismic slip rates, which renders this approach impractical for simulating earthquake sequences. On the other hand, our phase-field implementation for earthquake sequences provides a numerically robust framework that agrees with a discrete reference solution, is mesh objective, and reaches seismic slip rates. The model, unsurprisingly, requires highly refined grids around the fault zones to reproduce results close to a discrete model. Following this line, the effect of an internal length scale parameter on the phase-field predictions and mesh convergence are discussed.

rate regime with earthquake rupture and seismic wave propagation (in the order of 10 0 m/s), and a postseismic period accommodating rapidly induced changes before returning to tectonic loading in the interseismic period. Considering the highly nonlinear and multi-scale temporal nature of these phenomena, accurate numerical modeling of the dynamic earthquake rupture is a crucial step not only in earthquake applications [42], but also in broader instances that two bodies under dynamic friction contact are investigated, e.g. the stick-slip behavior observed in pile foundation and monopile installation under impact loading [75]. Rocks sliding across a contact are notably described by a frictional behavior that evolves throughout an earthquake sequence [e.g., 24,46]. A dynamic frictional contact governs the state of stress transfer between rocks and controls rupture nucleation, its subsequent propagation involving seismic waves propagating through the bulk media, and rupture arrest. Traditionally, friction coefficients during single earthquakes are described by a slip-dependent behavior (e.g., linear slip weakening friction) [39], which is not able to simulate recurrent stick-slip behavior nor capture the slip rate and contact history dependent friction behavior observed in laboratory experiments [e.g., 24,67]. Such variations are captured in rate-and-state dependent frictional formulations [46,63,76], which are a class of empirical, laboratory-based relations widely accepted within the earthquake community due to their success in capturing and predicting observations in different periods throughout earthquake sequences [40].
The numerical simulation of earthquake sequences is computationally challenging due to the ten to twenty orders of magnitude variations in spatial and temporal time scales involved [6]. Therefore it has conventionally been done by adopting simplifying physical assumptions and numerical tools, including efficient modeling using a spectral boundary integral formulation and anti-plane fault formulations [1,41,42]. Other contributions focus on enhancements of the physical problem, by improving the rheological description in the medium [e.g., 1,2,28] and fluid flow within the fault [62,74] or studying deformation from first physical principles on laboratory scales [e.g., 13,56,60,83]. The richness of physical mechanisms in the latter branch of models allows it to accurately study physical interactions, yet its comparatively higher computational cost inhibits extension to larger scales.
The common assumption in the numerical studies of earthquake sequence on continuum scales of kilometers and larger is that faults are infinitely thin, single-planar structures within highly simplified elastic media. In fact, faults are complex zones that are heterogeneous in terms of lithology and strain accommodation leading to highly variable mechanical and hydrological properties across their width [5,29]. Continental fault zones consist of several principal slip planes of 0.1-20 mm within a centimeters wide ultra-cataclasite embedded in a meters-thick foliated cataclasite, which is surrounded by highly damaged host rock of a few hundred meters [10]. In subduction zones, deformation localizes on finite thicknesses of less than 5 meters at the shallow portions [11], while being part of an up to a few kilometer wide subduction channel at a few to tens of kilometers depth [e.g., 12,77]. Moreover, fault zones intersect and are geometrically rough, while their geometry and location are regularly not well known. Representing this natural complexity as a single, infinitely thin fault at a predefined unknown location limits our ability to understand and forecast slip along faults, mechanisms governing slow slip events, and diffusion of fluids perpendicular to faults. As a result, these models can not reproduce certain features such as water lubrication within the fault [23], or the thermal pressurization which dynamically increases pore pressures during earthquakes [37], both causing a change in physical parameters controlling the fault. Some of these models are restricted as the fault's geometrical properties conform to the underlying grid layers, whereas this limitation has been relaxed by advancements in use of generalized finite element method for faults with simple fric-tion laws [14,44,59,65]. These enhancements are seldom applied in earthquake studies as the complexities associated with their implementations are rather extensive. Other recent contributions to facilitate fault geometrical description is the application of a phase-field approach to fracture propagation in geological materials [8,30]. These studies are however done for a simplistic frictional behavior, and similar to existing studies with the generalized finite element method, can not be used for modeling earthquake sequences.
To ensure flexibility of the numerical model and extension to realistic fault zone and regional architectures, a continuum representation of a fault is favorable in order to understand fault zones themselves and tectonic settings with complex fault geometries and/or in settings with no or little constraints on fault geometries. A continuum approach has gained recent traction in studies within the earthquake and tectonic communities [71,73,78]. An invariant plasticity formulation has been proposed to simulate complex and evolving fault zones with earthquake sequences using strongly slip-rate dependent friction [79]. This has been extended with adaptive time stepping to model earthquake sequences using rate-andstate dependent friction [36,57,58] or solid-fluid coupling [55]. These methods are spatially and temporally discretized using a fully-staggered, conservative finite difference method combined with marker-in-cell technique [32,34]. The finite difference method, although not as widely used as the finite element method, has distinct features that make it an appealing method, particularly in the earth sciences [32]. It removes the need for remeshing when combined with a marker-incell technique, which allows for an efficient combination of Lagrangian and Eulerian grids [34]. Other advantages include relying on a strong form, easier implementation, and some desired properties in the resultant stiffness matrices (such as a lower number of non-zero entries).
These models assume a fault is approximated by 1-2 layers of elements, which introduces mesh sensitivity to the fault geometrical distribution (e.g, its width). In these models, earthquake sequence characteristics are convergent upon grid size reduction for predefined, fixed width faults [36,79] due to an implicit regularization in which a rate-dependent formulation makes the slip rates objective [51]. However, it is expected that the mesh-dependent fault width description may introduce numerical instabilities upon applying it to complex fault geometries that are not aligned with the grid [36] and numerical inaccuracies may arise when considering physics in which an accurate description of fault width and local gradient fields is desired. Also for a growing fault problem slip rates are not yet objective during the initial strain localization stage, which introduces a grid size dependence for fault angles [58]. This is remedied in an ad-hoc manner using a heuristic fix in which fault width is dependent on slip rate [57], as inspired by findings of increasing localization in laboratory experiments [e.g., 21,70]. Nonetheless, grid size convergence and an accurate representation of fault widths can not be obtained without an internal length scale. Absence of an internal length scale and local loss of ellipticity not only introduces inaccuracies in fields and mesh-dependency, it also leads to numerical instability of the boundary value problem.
Mesh dependency and the absence of an internal length scale have been thoroughly investigated in various studies, and possible treatments have been discussed. Generally, a localized fracture is simulated in a diffuse and regularized manner to retain the mesh objectivity of the numerical solutions. Among different approaches, gradient enhanced damage model [54], Cosserat theory [49], Kelvin viscosity treatments [17,18,27,84] and more recently, the phase field approach to fracture [7,31,47,53,82] are proposed for modelling the localized fracture in a diffuse and regularized manner. Notably, phase field fracture methods have been extended to accommodate Drucker-Prager failure surfaces for both initiation [19] and growth [50].
In this paper, we investigate the applicability of two of these computational mechanics solutions and compare them to updated continuum and existing discrete approaches for describing earthquake sequences. We start by reformulating the plasticity based model of Herrendorfer et al. [36] to a consistently linearized classical plasticity formulation, highlighting its performance as well as contrasting its predictions against a conformal finite element model with zero-thickness interface elements [68]. This comparison demonstrates the validity of continuum models for earthquake sequences and suggests refinements up to a few hundred meters are appropriate to simulate earthquake sequences. To incorporate an inherent length scale, we extend our consistently linearized plasticity model with a Kelvin viscosity regularization to find that fault widths are regularized, yet seismic slip rates are inhibited. Finally for the first time we include rate-and-state friction formulation into a phase field approach and demonstrate objective earthquake characteristics as well as fault widths are obtained. These detailed comparisons will benefit the accurate simulation of multi-scale fault zones [e.g., 15,43] and earthquake sequences in regional and plate tectonic contexts [e.g., 16,33,80].

Numerical modelling of earthquake sequences
Let ∈ R d be a finite domain as shown in Fig. 1a, which describes a solid host rock of size L × L with the external boundary ∂ . The rock is homogeneous for simplicity in this study, and contains a pre-existing fault zone with an internal interface surface denoted by I .
To describe the mechanics of the problem, the momentum equation is considered: where σ i j represents the total stress tensor (σ i j = τ i j + δ i j P, with deviatoric (τ i j ) and hydrostatic (P) contributions distinguished), andü i is the second time derivative of a displacement component (i, j = x, y for the present 2D case). In this study, we ignore gravity and consider a compressible material with a finite bulk modulus. Additionally, heat conservation and shear heat production (shear heating), which is related to the dissipation of mechanical energy during frictional contact, is not considered. These uncoupled assumptions are frequently made in modelling of earthquake sequences [3], although the rate-dependence in rate-and-state friction empirically parameterizes some of the bulk weakening effects of temperature [22]. An elastic material behavior is assumed for the bulk, while a non-linear fault behavior is introduced using the rate-andstate dependent friction formulation as will be detailed.
Interfacial shear tractions t c are described through a frictional contact formulation: where μ is the friction coefficient which relates the shear and normal tractions. Fault frictional behavior is generalized following observations on rock contact from laboratory experiments, including so-called "slide-hold slide" and "velocity stepping" experiments. These form the basis for an empirical friction formulation in the framework of rate-and-state friction (RSF), which relates the friction coefficient to the magnitude of the slip rates V and a so-called state variable θ : with the parameters a and b the proportionality constants, μ 0 being the reference friction coefficient defined at a steady state slip rate V 0 , and δ being the characteristic slip distance. The second term in Eq. (3) introduces the instantaneous "viscosity-like" direct effect, because it represents the immediate response of μ, and hence τ s , to a change in V, which is proportional to a. The last term is referred to as the evolution effect, which is proportional to b and is described by the evolving state variable θ . θ is described as the state variable within the earthquake community, as it describes the state of the interface rather than being a solution variable and has amongst others been interpreted to describe the time that both fault sides have been in contact [25]. Different evolution laws have been proposed to parameterize the change of θ as a func-  [4], including the aging law [67] based on the observation of time-dependent healing at stationary contact, the slip law where state evolution only occurs when the fault slips, and the Nagata law [48]. Earthquake sequence simulations commonly apply one of the first two laws and we adopt the aging (or Ruina or slowness) law as in Herrendörfer et al. [36], which is defined as Analyzing a steady-state solution (i.e., dθ dt = 0) of Eq. (3) readily shows that upon an increase in slip rate, the friction coefficient decreases (i.e. rate-weakening) if a − b < 0, which allows for earthquake nucleation and facilitates earthquake propagation. Friction and fault strength increase with slip rate if a − b > 0 (i.e., rate-strengthening), which generally does not allow for earthquake nucleation and eventually leads to rupture arrest.
Equation (3) is singular for a limit case of V = 0 and for V V 0 can lead to negative friction coefficients. A regu-larized, commonly used version, which is largely equivalent for slip rates approaching V 0 , circumvents these limitations and is adopted here [42,63]: where all the parameters are identical to what described for the original form in 3.

Discretization
The problem consists of a frictional contact between two sliding blocks. To uniformly deform the rocks, we apply non-zero shear displacements at a constant velocity (v x ) along the top and bottom edges. Zero vertical displacements are assumed over all the domain edges. The boundary value problem is discretized in time and space, where in formulating the system of equations, an implicit approach is chosen with the displacements (u) as the primary global unknowns. For tem-poral discretization, the Newmark scheme [52] is applied and requires that first and second-order time derivatives of the displacement at the current time step be written as a function of the displacements at current and previous time steps: where superscripts n + 1 and n refer to the current and previous time steps, respectively. The Newmark constants are defined as: − t n is the time increment and β and γ are constant parameters. The Newmark scheme is stable provided that β > 0.25(0.5 + γ ) 2 , and γ ≥ 0.5. Here, we use β = 2 and γ = 1.5.
The spatial discretization of the problem is done by the use of standard quadrilateral or triangular grids, where displacements, as the main unknowns, are defined at nodes. Adopting the Voigt notation (for a plane strain condition and defining strains as: = x x yy xy T , where infinitesimal strains are defined as: i j = 1 2 (u i, j + u j,i )) and irrespective of the nature of the spatial discretization technique, strains are approximated at any arbitrary domain point x as: with d the n × 1 vector of nodal displacement unknowns, attributed to the nodes located inside the support covering point x (i.e. its parent element nodes, as shown in Fig. 1). The matrix B is the well-known derivatives matrix containing derivatives of the corresponding shape functions φ i , (i = 1, . . . , n) [86]. With this definition, a discretization method of either weak (e.g. FEM) or strong (e.g. FDM) form can be applied to discretize the boundary value problem. Following standard procedures and introducing (7) into (1) (or its weak counterpart) and expressing equations at the global level, the final discretized system of equations is written as: where u x and u y are the global displacement vectors in x and y directions, which are formed following a standard assembly procedure. Stiffness contributions K xx , K x y , K yx and K yy are the stiffness matrices that can be specified according to the applied numerical discretization technique. Throughout this paper, we use both finite element and finite difference methods, and actual expressions for the stiffness matrices are not presented for brevity and can be derived following standard textbooks for each numerical discretization method [32,86].
We use the Newton-Raphson method to iteratively solve the nonlinear system of equations in (8): where i refers to the current iteration number of the nonlinear solver, and the residual vector consisting of the current and previous time step displacements, velocities and accelerations is defined as: (10) and the Jacobian matrix is accordingly expressed as:

A discrete finite element method with zero thickness interface elements for frictional faults
The most common way to simulate an earthquake sequence is by treating a fault as a discrete discontinuity, i.e., of zero thickness or infinitely thin. A standard discrete fault model is described here as a reference model to validate numerical results throughout the paper. The model is discretized using the finite element method, where mesh and fault geometries conform at the edges. The interface region is defined by conventional zero-thickness interface elements located between the two blocks (dashed line in Fig. 1a) to introduce the interfacial frictional behavior. A magnified view of the discretized region is shown in Fig. 1c, where the fault is denoted by the thick gray line. Allowing only the tangential sliding to occur between the two faces of the rock at the interface, the contribution of the interface cohesive tractions acting on l on the total energy is [61]: where w c is the slip displacement between the two sides of the interface, and t c is the shear component of the traction vector acting along the interface. Within each interface element and at a sample integration point denoted by gp in Fig. 1c, the slip displacement can be evaluated as: where interface shape functions and displacement vector are defined through: and

A discretized plasticity model with mesh-adjusted fault representation
Generation of a conformal mesh is a computationally expensive task in many applications and a limiting factor for the discrete fault model. A continuum representation of faults should bypass this limitation by allowing the fault geometry to be independent of the grid. We adapt the continuum, rate-and-state friction model based on Drucker-Prager plasticity of [36] to ensure it is consistently linearized, such that in terms of convergence behavior in a non-linear scheme, it outperforms the original model using Picard iterations [36]. Following a plasticity formulation, we assume that the total strain rate is decomposed into an elastic˙ e and a viscoplastic contribution˙ vp along the region that corresponds the faults (blue elements in Fig. 1d) [27]: where the visco-plastic strain rate is zero at the host rock and is defined along the fault region similar to a rate-independent plasticity model as: with strain and stress components defined in indicial notation, and g is the plastic flow rule, andλ is a plastic multiplier. Assuming Drucker-Prager plasticity with zero cohesion [20], the rate-dependent yield function is expressed as where following the consistency plasticity model, the yield function is complemented by the Kelvin viscosity (η vp ), J I I is the second invariant of the deviatoric stress tensor and μ is the internal friction coefficient. The plastic flow rule for the case with non-associated plasticity is defined as where ψ is the dilation angle and is taken as zero in this study. The internal friction coefficient (μ = μ(θ, V )) is defined according to rate-and-state friction (Eq. (5)). Following [36], we adapt this equation to a continuum model in which tractions defined on the interface are replaced by invariants of the stress tensor throughout the bulk, i.e., shear and normal tractions are replaced by the second invariant of the deviatoric stress tensor and the normal stress, respectively. A similar approximation is made for slip rates, which are written as a function of the second invariant of the deviatoric plastic strain rate˙ I I ( p): vp assuming slip occurs across one cell [79], such that: where W f a length parameter taken as the size of the grid cells and a representation of fault thickness. Whenever the resulting yield function (18) is violated, corresponding stress values are returned to the yield surface (i.e. f =0) by using a standard implicit backward-Euler scheme. Note that for the present model plasticity along the fault layer is always present irrespective of the time and the magnitude of applied strain.

Consistent tangent operator
To achieve quadratic convergence of the Newton-Raphson solver discussed in Sect. 2, a consistent linearization formulation requires that the relationship between stresses and strains is identified in their incremental form.
Following standard procedures detailed in [35,69], the consistent tangent operator is defined as: (22) where: with D p the consistent tangent operator that reduces to the isotropic linear elastic stiffness tensor D el in the case that plastic deformations are absent. In the above expressions, as the time increment, and I the fourth-order symmetrization tensor,

Plastic multiplier increment 1
For evaluation of the plastic multiplier increment ( λ), the yield condition (18) is satisfied during the simulations (i.e. f =0). For the present case using Drucker-Prager plasticity [20], the yield function can be written in residual form by making use of the trial stress space: where the dilation angle (ψ) is assumed to be zero and G is the shear modulus. The plastic multiplier (λ n+1 ) and the current state parameter (θ n+1 ) are expressed aṡ Considering Eqs. (20), (24), (25) and (26) becomes a nonlinear function of λ. Thus values of the plastic strain rate as a function of λ can be obtained by a local iterative solution algorithm.

Continuum models with predefined fault thickness
The model presented so far assumes that the fault thickness is equal to one grid layer (Fig. 1d), which is the size meshdependent models eventually localize to. It lacks an internal length scale for the fault and thus a geometrically objective strain field around the discontinuity can not be attained. In this section two methods are presented for the modelling of faults in a continuum setup, which each incorporate an internal length scale and assume a finite thickness for the fault that is independent of the grid. For the first method, we extend the explained model with a simple modification to represent the fault zone at predefined widths. In the second method, we developed a model that describes faults governed by rate-andstate dependent friction using a phase field representation.

A plasticity based model using Kelvin viscoplasticity
It is known from studies in earth science applications that incorporating a Kelvin viscous rheology (i.e., a parallel arrangement, see [27], Fig. 1) introduces an internal length scale in problems that involve strain localization or nonassociated plasticity. Since the original continuum models proposed are similar in form to the geodynamic models to which this extension was proposed, we study whether a relatively simple adaptation can ensure grid-independent fault thicknesses. To allow the finite thickness to be governed by the Kelvin viscosity parameter, we assign a region of weak layer with thickness W f (Fig. 1e) to represent the fault layer.
To extend the fault layer beyond one nodal layer, fault slip rates need to be integrated over the fault thickness (W f ): This simple adaptation allows the fault localization zone to be governed by the Kelvin viscosity, and not by the grid itself.

Phase-field methods
In the phase-field method, also known as the regularized Griffiths' brittle theory, a sharp discontinuity is translated into a smooth field that is defined independently of an underlying grid. Although the phase-field approach has fundamental similarities to the plasticity-based models discussed in the previous section, its formulation relies on energetical concepts. The total potential energy is subdivided into an elastic strain energy term and a surface energy term associated with the formation and evolution of cracks. In essence, the phasefield representation of a sharp discontinuity (i.e., a fault) can be done by defining an internal length scale l s and a phase field variable φ (Fig. 1f). These two quantify the widening or thickness of the damage zone. The phase-field variable φ defines the state and evolution of damage differentiating between a fully damaged zone (φ = 1), an intact zone (φ = 0), and a partially damaged zone (0 < φ < 1). φ is an additional degree of freedom to the system and is coupled to the deformation field. When applying a phase field approach, a regularized representation of a discrete discontinuity ( l ) can be written as [53]: where the planar crack surface density function is defined as for a vanishing length scale (l → 0), this term reproduces a discrete crack line. The energy functional E of the boundary value problem that was discussed in Sect. 3.1 can be expressed as where u is the total strain energy, s refers to the surface energy due to the formation of cracks, and b sums the remaining contributions that are unchanged and due to external loading. For a linear elastic bulk material, the total strain energy can be expressed as where C 0 is the linear elastic stiffness matrix, is the standard elastic strain tensor defined in Eq. (7), and g(φ) is the degradation function due to damage evolution, here taken as where k is a small parameter and is added to ensure a wellconditioned system for the broken state as φ gets close to 1.
Defining g c as the critical Griffith-type energy release rate of the solid material per unit area, the fracture energy due to formation of crack surfaces is written as The total potential energy functional can therefore be defined as a volume integral over the bulk material After using Gauss' divergence theorem, one can derive the final form of the coupled u − φ equations After writing the equations in a weak form, a non-linear solution scheme detailed in Sect. 2.1 is adopted.
Rate-and-state friction in a phase-field method The phasefield approach was initially proposed to simulate brittle fracture, which does not account for stress transfer between the discontinuity edges. However, earthquakes occur on faults located at large depths and pressures inside the Earth's lithosphere, such that both sides are always in contact and a frictional formulation (e.g., rate and state friction) applies.
Recently, various attempts were made to formulate a phasefield approach that can be used to model cohesive fractures [9,53,85]. In a general interpretation, this requires the addition of cohesive terms to the total energy functional where the interface integrals are transformed into a phasefield volume integrals. For a stationary fault considered here, evaluating original line integrals is computationally more efficient. We will demonstrate that these approximations still lead to very good agreement with the reference solution of a discrete fault.
As shown in Sect. 2.2 for the rate and state friction formulation, exchanged tractions between the interface surfaces are expressed as a function of slip rates. Therefore the most important ingredient of the phase-field formulation for defining RSF faults remains in evaluating the local gap vectors. Smeared displacement jump approximation One of the main computational steps in using a phase field approach to cohesive fracture is to represent a sharp discontinuity by a smooth representation. This step requires the definition of discontinuity opening/sliding displacements (or velocities). In the literature, different methods have been explored, some of which rely on analytical expressions and are limited to only simplified setups. However, the most practical approach relies on a heuristic level set framework by assuming that no unique slip profiles exist for a smeared discontinuity representation, but that it can be approximated [53]. In this way, the gap vector is approximated by evaluating the displacements at specific locations around the discontinuity (Fig. 1f). This implies that the evaluated slip values are subject to change depending on the sample points' location, taken as an input to the problem. Therefore it is necessary to thoroughly study the effect of these sampling locations on the solutions. It should be mentioned that the advantage of the adopted level set approach is that displacement gaps can be conveniently evaluated at the sides of a crack/fault. However, the proposed formulation can be extended to the line integral approach [85] at some additional cost with more complications in the implementation/linearization of the arising terms.
By defining an offset distance which is taken as a function of the length scale parameter in the phase-field representation (l o = α φ l s ), one can define the slip displacements for a sample gauss point (Fig. 1f) as where w t and w b represent the value of the tangential displacement at locations denoted in Fig. 1f at the top (t) and bottom (b) sides of the discontinuity, and s represents the local coordinate vector along the interface. The location of these sample points is chosen by projecting an integration point on the line of discontinuity (or a volume integration point in case the cohesive terms are evaluated in a volume format) to the two lines representing the top and bottom edges (dashed lines in Fig. 1f). The parent element of each of the two points is then identified using an efficient octree structure [26], and the corresponding coordinates of the sampling points in a local coordinate system are evaluated using an inverse mapping technique. The local displacement vectors are then evaluated as

Model setup
We simulate a 2-D plane strain setup of 150 km×150 km representing a strike-slip fault, comparable to the setup used in [36] (Fig. 1a). Initial displacements and velocities are set to zero throughout the bulk media. A small uniform pressure of 5 MPa is predefined all over the domain, representing compressed rocks located at a few kilometers depth. The two blocks are then deformed under monotonically increasing shear displacements with a constant velocity of V 0 = 2 × 10 −9 m/s as described for the boundary conditions in Sect. 2.
The material properties for the fault and host rock are chosen within a range that is usually adopted in classical earthquake cycle simulations ( Table 1). The rate and statedependent friction parameters are adopted from [36]. The central segment of the fault is characterized as a steady-state rate-weakening zone with a − b < 0, where the length of the seismogenic zone at the center is set equal to 76 Km. This segment is bounded by two segments of 32 km and 42 km (shown as shaded gray region in Fig. 1a) described by rate strengthening friction to inhibit rupture propagation to the boundaries and thereby minimize boundary artifacts. The asymmetrical configuration ensures that the rupture initially nucleates at the right limit of the seismogenic zone.
In the continuum models, the fault (red shaded region in Fig. 1b) is identified by a very small initial state variable, which indicates recent movements along that fault and allows plastic deformations to localize there. For the bulk medium a very large initial state variable is set (θ 0 = δ V 0 exp(40) s), such that off-fault plasticity is neglected.

Results
In the simulations, the time step ( t) is adaptively adjusted to resolve all stages of an earthquake sequence, namely the interseismic (years to more than ten thousands of years), coseismic (milliseconds or less to minutes), and postseismic (minutes to years) phases. Following [36,42], the time increment t is inversely proportional to the maximum slip velocity along the fault line (V max ) and defined as where θ max is in general a function of pressure distribution at the rock [36]. Unless specified otherwise, we use a fixed value of θ max = 0.2 during the simulations [36].

Continuum slip-rate dependent plasticity model without added length scale
The standard model is discretized using a non-staggered finite difference discretization method (Fig. 1b) and is used here for the simulation of earthquake cycles. Kelvin viscosity parameter is set to zero (η vp = 0), and the rate-and-state dependent friction formulation characterizes rocks in frictional contact during episodes of stick-slip behavior. The widely different time scales are captured using adaptive time stepping (Fig. 2b) within a time-dependent non-linear system of equations, where the yield function (18) is satisfied at all times, and plastic deformation is evaluated in the form of plastic slip rates across the fault zone. As such, and unlike more traditional friction formulations, the RSF model predicts non-zero slip rates at all times (Fig. 2a). Due to loading at the far-field boundaries, slip rates increase at an exponential rate within the velocity strengthening zone, defined by the direct effect of RSF, until the maximum slip rate in the velocity strengthening zone reaches the loading rate equal to the reference rate V 0 (Fig. 2a). During this initial loading stage, the velocity weakening zone hardly experiences any slip (Fig. 2c). After this initial loading stage, quasi-characteristic earthquakes recur periodically. Below we summarize the evolution of an earthquake sequence, as we refer the reader to [36] for a detailed description of the operating physics. During the interseismic period solved using large time steps (Fig. 2b), slip propagates from the rate-strengthening segments into the central rate-weakening segment (black lines in Fig. 2b). Once slip within the rate-weakening segment occurs over the critical nucleation size [e.g., 66], a dynamic instability or earthquake nucleates. During this coseismic period, inertia controls and limits slip rates, rupture speeds, and seismic waves propagation throughout the domain (Fig. 2d).
The earthquake is characterized by a sharp, orders of magnitude increase in slip velocities over a period of seconds.
The minimum time step in the analysis is on the order of milliseconds (Fig. 2b), and maximum slip rates along the fault increases to values up to meters per second (Fig. 2a). Due to the non-symmetric configuration, rupture initiates from the right tip of the velocity weakening zone and propagates towards the left. As the rupture propagates, slips accumulates (blue lines, Fig. 2c) and the state variable decreases by orders of magnitude. Seismic waves are emitted due to variations in rupture speed (Fig. 2d) and can affect rupture propagation, particularly upon reflections [38,45,81]. In this model setup and without absorbing boundary conditions, waves reflected from the boundaries affect the rate-weakening segment at the end of the rupture and only cause very mild variations in stress unable to perturb the periodic earthquake sequence [36]. The last phase of an earthquake sequence, the postseismic period, is a phase of relaxation, where increased slip rates and stresses in the surrounding rate-strengthening segments decrease again. After this re-equilibration, the whole system starts to prepare again for the next earthquake sequence.

Comparison with discrete fault model
The numerical results of the continuum, finite-difference model are compared to a reference discrete fault, finite element model, which uses a zero thickness interface formulation (described in Sect. 2). The two simulations are run under similar initial and boundary conditions, although small discrepancies may arise from inherently different boundary condition implementations in FDM and FEM models. Results for slip and slip rates are shown at time steps of 50, 150, 200, and 300 of each simulation (Fig. 3). The fault zone distribution in the continuum model is described by a center line (blue) and an edge line (red). The center line of the continuum model shows good agreement with the slip observed at the infinitely thin fault line (Fig. 3a-d). In terms of slip rates, the continuum model adequately represents the discrete finite element model (Fig. 3e-h). Existing differences can be attributed to the kinematic differences between the models, and that in the continuum representation of the fault, nonunique slip profiles are defined, such the single line slip rate values need to be combined for an accurate representation of the fault zone. Differences are particularly minor in slip rates considering that due to slip-rate dependent time stepping (eq (40)), attributed times have not been kept identical between the different grids.

Grid refinement study
In the slip-rate dependent plasticity model, the fault thickness is inherently adjusted to the mesh, such that we also set the fault thickness equal to the size of one element (L f = y, hence, the mesh-adjusted fault representation). Considering a single earthquake, the significance of this assumption on the quality of numerical results is evaluated. To limit the number of considered time steps, the value of θ max is set equal to 1.
A mesh refinement analysis using three different grid resolutions is considered, where the grid size is set equal to L/α, with α = 150, 300, 600. Slip curves show a convergent pattern at the coseismic phase (Fig. 4a). Similar convergence is observed for the slip rates, which are shown in Fig. 4b for the high slip rate in the coseismic period. Corresponding time step values are indicated in Fig. 4c at all times of the dynamic rupture simulations. All curves confirm that the continuum, slip-rate dependent plasticity model can adequately predict the relevant quantities of interest for earthquake rupture and earthquake sequence applications (i.e., those proportional to slip rates, including recurrence interval). Although the assumption of grid-dependent fault thickness implies that the model lacks an internal length scale, the rate-and-state friction formulation [36] is cast in a format that slip rates are normalized with respect to the grid size [51,79]. This is evident in results in Fig 4a,b, which show that the slip-rate dependent plasticity model provides objective estimates for global slip and slip rates upon refining the grid.
However, from a mathematical point of view, the model is not objective, as the local plastic strains do not follow the same width and maximum upon refinement (Fig. 4d). Even though upon integration across the fault zone the amount of strain converges, the maximum strain keeps increasing sharply upon reducing the grid size. This is attributed to the fact that the standard model lacks an internal length scale that can bound strain rates.

Fixed fault thickness by adding Kelvin viscosity
An investigation of the model described in Sect. 2.4.1 is presented. The fault layer is predefined such that it covers multiple grid cells (Fig. 1e) and has a fixed thickness equal to L/4000 (irrespective of the grid). To ensure that this small fault thickness is discretized adequately and with multiple elements, we generated a locally refined triangular grid using the Bisection algorithm [64]. Introducing a finite thickness for the fault does not include a length scale for the fracture on its own. For regularization, the Kelvin viscoplasticity material model is used, which is shown to introduce an internal length scale in some applications in earth science [27].
As mentioned in the numerical method section, slip values have to be estimated through an integral form (27). For ease of numerical implementations in this section, we simplified the slip expression into a closed-form relation similarly defined for the standard model but with a different thickness: V ≈ W f˙ I I (vp), where W f is used in this expression instead of the grid size in the standard model, and represents the actual fault thickness (W f = L/4000).
Since Kelvin viscosity is mainly chosen based on numerical considerations, we investigate Kelvin viscosities over a wide range from η = 2.5 × 10 3 Pa s up to η = 2.5 × 10 18 Pa s (Fig. 5). Maximum slip rates are calculated during the simulations and are shown in Fig. 6. For each Kelvin viscosity, we use three different grids with minimum sizes of L/4000/h, with h the number of elements that cover the fault thickness, and is set equal to 4, 8, and 16 for the three adopted grids, respectively. We report the grid size convergence of the second invariant of the strain tensor distribution in Fig. 5, where the same quantities are plotted along a vertical cross-section in Fig. 7a-f. We observe that for smaller kelvin viscosities, strain still tends to localize into one grid upon refinement of the grid, demonstrating that fault width is still grid-size dependent. As the kelvin viscosity parameter increases, the length scale effect becomes clear. Strain distributions reach mesh objectivity for Kelvin viscosities that exceed η = 2.5 × 10 15 Pa s. If we analyze the maximum slip rate through time for this fault system, we observe that for smaller Kelvin viscosities, slip rates still show a sudden increase towards seismic slip rates beyond the applied loading rate right after the initial loading phase (Fig. 6). By gradually increasing the Kelvin viscosity values, the graphs clearly show that the slip rates plateau around the loading rate V = 10 −9 m/s. For Kelvin viscosities larger than 2.5 × 10 15 Pa s, they thus reach steady-state creep rates many orders of magnitude below seismic slip rates (i.e., 10 −2 up to 10 1 m/s). This suggests that the large Kelvin viscosities act as a diffusive element, which inhibits the occurrence of earthquakes.
These observations are in line with recent studies on applications of the Kelvin viscosities in dynamics [72], where for a rectangular region under dynamic loading it was found that viscoplasticity in the presence of inertia does not prevent mesh dependency. We extend that to a more sophisticated earthquake application, where seismic waves have a dominant effect. Regularization is only present for a steady-state situation, in which the viscosity parameter is high enough to inhibit seismic rates and earthquake occurrence.

A phase field model for earthquake sequences
We simulate loading and seismic slip up to slip rates of ≈ 0.3 m/s in the reference strike-slip model setup (Sect. 2.5) using the phase-field method described in Sect. 2.4. The phase-field form is discretized using the finite element method, along with a penalty approach used to impose the pre-defined phase-field values along the fault line. φ = 1 is enforced as an initial condition at equidistance locations along the fault line, l . As we limit this study to the case of a stationary fault, the phase-field equation is solved once and the phase-field variable is used throughout the dynamic problem. However, the formulation described in Sect. 2.4 is general and can be used for non-stationary or evolving faults as well. To reproduce results close to a discrete reference model with zero thickness, relatively small length scale parameters are used. To resolve such small length scales and adequately discretize the fault thickness, a locally refined and structured Q4 finite element grid is used.
In contrast to what is mentioned in Sect. 2.4.1, phasefield representation of the fault does not take the fault thickness as an input parameter. Instead, the geometry and thickness of the fault are proportional to the length-scale parameter l s . In the numerical simulations, the length-scale  (Fig. 8). For all length-scale parameters, we observe an increase in maximum slip rates up to seismic slip rates. Our fault zone simulated using a phase field approach is thus able to generate dynamic earthquake ruptures. The exact maximum slip rate reached after 1500 time steps averages at 0.3 m/s for the widest fault (L/500), while it reaches about 0.315 m/s for the thinnest fault (L/6000). This small difference suggests that a bit more energy is used to break a wider fault zone, such that that energy could not be used towards accelerating the rupture.
A convergence analysis with minimum grid sizes of 20 m (L/7500), 10 m (L/15,000) and 5 m (L/30,000) is performed for each length-scale parameter as well (different lines in Fig. 8). We note that the fault thickness is calculated based on the phase-field parameters and is defined over multiple grid cells irrespective of the grid, i.e., not bounded by a single grid layer as used for the standard model. The numerical results show convergence with increasing grid size for slip rates and a mesh objective solution for all length-scales. This is also unlike the behavior reported for the fixed thickness representation of faults in the previous section. Grid size convergence is attained faster for the larger length scales (e.g., Fig. 8a). Moreover, for smaller length-scales, a finer grid is required to obtain a converged solution (e.g., Fig. 8d). These numerical results confirm a convergent pattern at all time steps and for all local variables including the fault thickness.
A major challenge in the phase field representation of fault zones governed by rate-and-state dependent friction is the evaluation of slip rates along the fault zone, and specifically at the integration points of the cohesive integral contribution. As detailed in Sect. 2.4.2, slip is not uniquely defined across the fault following the level set approach. Instead it is sampled at specific locations identified by a user-defined offset distance l o . Thus the approximated slip (rate) values are used for the evaluation of rate-and-state friction coefficients in Eq. (5) and throughout the simulations. To test the impact of slip sampling point location, we vary the offset distance l o = αl s by changing α, which changes the location of the sampling points. The exact slip and slip rate values at specific times are sensitive to the location of these sampling points for length-scale parameters of L s = L/3000 or 50 m and larger, as peak slip values will differ by more than 10% (Fig. 9).
Phase-field results are also compared with the reference finite element model, which assumes a discrete fault with zero thickness (black line in Fig. 9). Phase field results for an offset length of α = 0.5 do not agree well with the discrete fault results, which can be attributed to the fact that offset is smaller than the length scale value. The agreement becomes better as the offset length α is increased, and a value larger than the length scale is chosen. Clearly, all the values larger than the length scale (even the largest ones) provide a good approx-imation of the slip and slip rate values. However, the range of slip curves considerably narrows down as the length scale parameter becomes smaller. This confirms our observation that the impact of offset distance is only pronounced when relatively large length scales, i.e., large fault thicknesses, are chosen. Overall, upon using thin enough fault zones (i.e., about a 100 m or less for the rate-and-state friction parameters) combined with small enough sampling distances (i.e., an off-set length of at least one length-scale) results of our phase field fault zone approximation agree extremely well with the discrete fault results.

Summary and conclusions
We simulate faults as predefined regions of finite thickness, which are governed by a rate-and-state-dependent friction formulation to simulate fully dynamic earthquake sequences within a continuum model. The adopted continuum approaches readily allow for the computationally efficient simulation of geologically realistic structures, including the simulation of evolving faults and off-fault damage on both present day and geological time scales, simulating known and unknown faults in regional plate tectonic settings, and simulating geologically realistic fault zones that may involve multiple principle slip locations and/or experience complex fault geometries and interactions. To ensure robust numerical simulations in this growing field, we evaluate the numerical performance and mesh objectivity of different global and local variables characterizing earthquake sequences.
Following the classical theory of plasticity, we derived a consistently linearized version of a continuum, slip-rate dependent plasticity model recently proposed by Herrendörfer et al. [36]. In this model, discretized using a fully staggered, conservative finite difference scheme, the thickness of a fault is essentially defined by one grid cell. The fault slip rate is estimated as a function of the plastic strain rates normalized by that one grid size. The model is shown to simulate periodic sequences of fully dynamic, quasi-characteristic earthquakes and accurately reproduces the results of a discrete finite element model with zero thickness interface elements. Upon reduction of the grid size, key variables of interest, such as slip rates, slip and recurrence intervals are shown to remain objective. However, it is expected that the mesh-dependent fault width description may introduce numerical instabilities upon applying it to complex fault geometries not aligned with the grid and numerical inaccuracies may arise when considering physics in which an accurate description of local gradient reference fields is critical. In addition, for predefined or shallowly dipping faults, this model provides accurate results even for remarkably coarse resolutions, i.e., for grid sizes on the order of hundreds of meters instead of (tens of) meters for the same friction and elastic parameters. This supports and promotes the usage of this continuum rate-and-state friction formulation for regional and fault zone simulations in many scenario's, including regional and fault zone settings with shallowly or steeply dipping faults, such as subduction zones and strikeslip settings.
In an attempt to incorporate a mesh insensitive and finite width fault description, such that fluid flow across heterogeneous fault zones can also be accurately simulated, this plasticity based model was extended by incorporating a simple constant Kelvin viscosity term within the consistency plasticity formulation, which is known to introduce an internal length scale in many applications. Our results show that at high Kelvin viscosities, diffusion is strong enough to diffuse away accelerating slip rates and restrict them to values many orders of magnitude below seismic slip rates, hence inhibiting the occurrence of an earthquake. For smaller viscosities, seismic slip rate can be approached, but plastic deformation still localizes into one grid cell, as the viscosity effect as an internal length scale is not realized. This may be a result of the meter-scale grid sizes utilized. However, even if finer grids would be able to regularize dynamic slip, it is not feasible to simulate regional scales using sub-meter resolutions. These results are in line with recent observations regarding the use of Kelvin viscosity as a regularization in dynamical problems [72].
Finally, a phase-field approach for a cohesive fracture is implemented and for the first time applied to the modeling of earthquake sequences on faults governed by rate-andstate friction. Unlike plasticity-based continuum models for meter-scale grid sizes, we show that the phase-field model allows for a mesh objective finite thickness for the fault, while maintaining a seismic response at corresponding seismic slip rates on the order of m/s. The seismic slip and slip rates at specific time steps also agree very well with the discrete reference solution. The effect of a length scale parameter used as an input for the phase-field approach was discussed in detail, and it was concluded that offset length values larger than the length scale parameter give good approximations of the slip (rate) values. The optimum fault representation was found for an offset length equal to 2 × l f . However, it should be mentioned that the choice of this distance is a simplified way to approximate fault slips, and more importantly, the choice of ideal value depends on the problem. The results further confirm that the phase-field approach, among various continuum models discussed in this paper, is a strong candidate for further explorations in earthquake applications.
Funding This study is part of the research programme DeepNL, financed by the Dutch Research Council (NWO) under project number DeepNL.2018.033.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.