Casper Pranger
email: casper.pranger@geophysik.uni-muenchen.de

Patrick Sanan

Dave A May

Laetitia Le Pourhiet

Alice-Agnes Gabriel

Rate and State Friction as a Spatially Regularized Transient Viscous Flow Law

We reformulate the empirical rate and state friction law as a bulk viscous flow law 11 in terms of anelastic shear strain rate.

12

• We show how mesh independence is achieved by including a gradient-like non-local 13 anelastic shear strain rate equivalent. 14 • We show analytically and numerically that the proposed continuum model closely 15 reproduces existing results of rate and state friction.

show that the expected tendency of the continuum rheology for runaway localization and mesh dependence can be halted by including an artificial diffusion-type regularization of anelastic strain rate in the softening law. We demonstrate analytically and numerically using a simplified fault transect that important aspects of the frictional behavior are not significantly affected by the introduced regularization. Any discrepancies with respect to the interfacial description of fault behavior are critically evaluated using 1D numerical velocity stepping and spring-slider experiments. Since no new physical parameters are introduced, our model may be straightforwardly used to extend the existing modeling techniques. The model predicts the emergence of complex patterns of shear localization and delocalization that may inform the interpretation of complex damage distributions observed around faults in nature.

Plain Language Summary

How, where, and when earthquakes nucleate is one of the great questions in science and society that, despite steady progress, has hardly been answered to any practical degree. Based on field observations, laboratory experiments, and theoretical work it is believed that a cocktail of escalating mechanical, chemical, and thermal grain-scale processes cause the sudden and rapid onset of earthquakes. The net effect of these processes are characterized by an immediate strengthening and a gradual weakening response to deformation and are unified in simplified form in the theory of 'rate and state friction'. This theory is commonly used in computer simulations of earthquake sequences.

We point out that rate and state friction, unlike some physical theories of earthquake

Introduction

Over the second half of the last century the study of the sliding behavior of frictional surfaces, such as those believed to occur in the earth's crust, has led to a general understanding that it is governed by competition of stabilizing viscous-like effects and potentially destabilizing processes that affect the texture of the interface in a time-dependent manner and self-organize into periods during which the interface is arbitrarily close to elastic stick punctuated with periods of anelastic slip [START_REF] Rabinowicz | The intrinsic variables affecting the stick-slip process[END_REF][START_REF] Brace | Stick-slip as a mechanism for earthquakes[END_REF][START_REF] Bowden | Friction, lubrication and wear: a survey of work during the last decade[END_REF][START_REF] Scholz | Detailed studies of frictional sliding of granite and implications for the earthquake mechanism[END_REF][START_REF] Ida | Cohesive force across the tip of a longitudinal-shear crack and griffith's specific surface energy[END_REF][START_REF] Dieterich | Time-dependent friction and the mechanics of stick-slip[END_REF][START_REF] Ruina | Friction laws and instabilities: A quasistatic analysis of some dry friction behaviour[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF][START_REF] Dieterich | Direct observation of frictional contacts: New insights for state-dependent properties[END_REF].

Rate-and State-dependent Friction

This has led to the proposition of the phenomenological rate-and state-dependent friction law by [START_REF] Dieterich | Time-dependent friction and the mechanics of stick-slip[END_REF]Dieterich ( , 1979a)); [START_REF] Ruina | Friction laws and instabilities: A quasistatic analysis of some dry friction behaviour[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF], which reads in its general form as

f = f 0 + a log V + b log Θ. (1)
In this law, the friction coefficient f is given by the linear combination of a reference friction coefficient f 0 , a logarithmic contribution from dimensionless slip rate V multiplied by a coefficient a, and another logarithmic contribution from a dimensionless state variable Θ multiplied by a coefficient b. All quantities are taken to be positive.

The rate effect or direct effect is to a variable degree strengthening with increasing slip rate and thus exerts a stabilizing influence on the frictional interface for any strictly positive a. The state effect or evolution effect is to be governed by an evolution law that has the evolving steady state Θ ss V (t) = 1. Then, if b > a, this allows an externally loaded system with a frictional interface to move to a configuration with lower elastic strain energy and thus a frictional-mechanical instability may occur.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

One of the most widely used state evolution laws is the aging law of [START_REF] Ruina | Friction laws and instabilities: A quasistatic analysis of some dry friction behaviour[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF], given by

Θ = r 0 (1 -V Θ) , (2)
with r 0 a rate constant that is commonly expressed by dividing the reference velocity v 0 by a critical slip distance d c . It should be noted that several other evolution laws that govern the state variable have been proposed. The most commonly cited contender is the slip law (Dieterich, 1979a;[START_REF] Ruina | Friction laws and instabilities: A quasistatic analysis of some dry friction behaviour[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF]. In this work we do not consider it because its functional form ( Θ ∝ Θ V log Θ V ) is unfortunately not amenable to the type of analytical treatment that is performed. We note that rate and state friction captures only the first-order behavior of the stick-slip cycle and that laboratory experiments have brought to light many secondary effects (e.g. [START_REF] Ruina | Slip instability and state variable friction laws[END_REF][START_REF] Marone | The effect of loading rate on static friction and the rate of fault healing during the earthquake cycle[END_REF][START_REF] Mair | Friction of simulated fault gouge for a wide range of velocities and normal stresses[END_REF][START_REF] Chester | Effects of temperature on friction: Constitutive equations and experiments with quartz gouge[END_REF][START_REF] Passelègue | Initial effective stress controls the nature of earthquakes[END_REF].

Following Amonton's law the friction coefficient is expressed as the ratio of shear stress τ to normal stress σ, implying a cohesionless fault. It is further assumed that the fault is always critically loaded, i.e. shear stress is equal to the shear strength and V > 0.

The assumption of criticality is necessary because (1) degenerates at V = 0, but has also been suggested to be realistic for faults in nature [START_REF] Bak | Earthquakes as a self-organized critical phenomenon[END_REF].

The rate and state friction law is frequently applied in numerical studies of fault slip, whether that be a study over the course of a single earthquake or a complex sequence of slip transients (e.g. [START_REF] Dieterich | Modeling of rock friction: 2. simulation of preseismic slip[END_REF][START_REF] Okubo | Dynamic rupture modeling with laboratory-derived consti-tutive relations[END_REF][START_REF] Rice | Spatio-temporal complexity of slip on a fault[END_REF][START_REF] Rice | Slip complexity in earthquake fault models[END_REF][START_REF] Ben-Zion | Dynamic simulations of slip on a smooth fault in an elastic solid[END_REF]Lapusta et al., 2000;Lapusta & Liu, 2009;[START_REF] Rubin | Earthquake nucleation on (aging) rate and state faults[END_REF][START_REF] Gabriel | The transition of dynamic rupture styles in elastic media under velocity-weakening friction[END_REF][START_REF] Jiang | Deeper penetration of large earthquakes on seismically quiescent faults[END_REF]D. Li & Liu, 2017;[START_REF] Erickson | The Community Code Verification Exercise for Simulating Sequences of Earthquakes and Aseismic Slip (SEAS)[END_REF]. The use of interfacial friction laws in these models necessitates the treatment of faults as mesh features on which internal boundary conditions can be applied.

It can be laborious to construct such meshes for non-trivial fault geometries, and the procedure does not scale well with increasing fault network complexity. There is currently no established method to construct meshes with time-dependent geometry and topology reflecting an evolving fault system, but we refer to K. [START_REF] Okubo | Dynamics, radiation, and overall energy budget of earthquake rupture with coseismic off-fault damage[END_REF] for work in this direction. With some exceptions (e.g. D. Li & Liu, 2017;[START_REF] Sathiakumar | Earthquake cycles in fault-bend folds[END_REF][START_REF] Perez-Silva | 3d modeling of longterm slow slip events along the flat-slab segment in the guerrero seismic gap, mexico[END_REF] many studies focus on the already complex behavior of long-term sequences of slip transients on linear or planar faults in simple domains.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Continuum Fault Rheology

There is strong field and experimental evidence that faults in nature are not infinitely thin planes but consist of complex evolving networks of strongly localized shear zones within a wider region of damaged host rock [START_REF] Tchalenko | Similarities between shear zones of different magnitudes[END_REF][START_REF] Granier | Origin, damping, and pattern of development of faults in granite[END_REF][START_REF] Chester | Ultracataclasite structure and friction processes of the punchbowl fault, san andreas system, california[END_REF][START_REF] Katz | Geometry and kinematic evolution of riedel shear structures, capitol reef national park, utah[END_REF][START_REF] Savage | Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones[END_REF]Faulkner et al., 2011;[START_REF] Barth | Slip localization on the southern alpine fault, new zealand[END_REF][START_REF] Passelègue | Frictional evolution, acoustic emissions activity, and off-fault damage in simulated faults sheared at seismic slip rates[END_REF][START_REF] Perrin | Off-fault tip splay networks: A genetic and generic property of faults indicative of their long-term propagation[END_REF][START_REF] Locatelli | Intermediate-depth brecciation along the subduction plate interface (monviso eclogite, w. alps)[END_REF][START_REF] Locatelli | Fluid pulses during stepwise brecciation at intermediate subduction depths (monviso eclogites, w. alps): First internally then externally sourced[END_REF][START_REF] Pozzi | A new interpretation for the nature and significance of mirror-like surfaces in experimental carbonate-hosted seismic faults[END_REF][START_REF] Pozzi | Coseismic ultramylonites: An investigation of nanoscale viscous flow and fault weakening during seismic slip[END_REF][START_REF] Pozzi | Coseismic fault lubrication by viscous deformation[END_REF]Ritter, Rosenau, & Oncken, 2018;[START_REF] Ritter | Sandbox rheometry: Co-evolution of stress and strain in riedel and critical wedgeexperiments[END_REF][START_REF] Ross | Hierarchical interlocked orthogonal faulting in the 2019 ridgecrest earthquake sequence[END_REF]. This knowledge has spurred a search for intermediate-scale continuum homogenizations of distributed micro-scale processes that obey the rate and state framework in a somewhat generalized sense. Currently proposed models focus on rheological feedback mechanisms involving temperature [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] Braeck | Spontaneous thermal runaway as an ultimate failure mechanism of materials[END_REF][START_REF] John | Generation of intermediate-depth earthquakes by selflocalizing thermal runaway[END_REF][START_REF] Garagash | Seismic and aseismic slip pulses driven by thermal pressurization of pore fluid[END_REF][START_REF] Roubek | A note about the rate-and-state-dependent friction model in a thermodynamic framework of the Biot-type equation[END_REF][START_REF] Rice | Stability and localization of rapid shear in fluid-saturated fault gouge: 1. linearized stability analysis[END_REF][START_REF] Platt | Stability and localization of rapid shear in fluid-saturated fault gouge: 2. localized zone width and strength evolution[END_REF][START_REF] Thielmann | Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating[END_REF][START_REF] Thielmann | Grain size assisted thermal runaway as a nucleation mechanism for continental mantle earthquakes: Impact of complex rheologies[END_REF][START_REF] Pozzi | Coseismic fault lubrication by viscous deformation[END_REF], grain size [START_REF] Rozel | A thermodynamically self-consistent damage equation for grain size evolution during dynamic recrystallization[END_REF][START_REF] Thielmann | Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating[END_REF][START_REF] Thielmann | Grain size assisted thermal runaway as a nucleation mechanism for continental mantle earthquakes: Impact of complex rheologies[END_REF][START_REF] Barbot | Modulation of fault strength during the seismic cycle by grainsize evolution around contact junctions[END_REF][START_REF] Pozzi | Coseismic fault lubrication by viscous deformation[END_REF], porosity [START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF][START_REF] Niemeijer | A microphysical model for strong velocity weakening in phyllosilicate-bearing fault gouges[END_REF][START_REF] Chen | Rate and state frictional and healing behavior of carbonate fault gouge explained using microphysical model[END_REF]Van den Ende et al., 2018), fluid pressure [START_REF] Barbot | A unified continuum representation of post-seismic relaxation mechanisms: semi-analytic models of afterslip, poroelastic rebound and viscoelastic flow[END_REF][START_REF] Garagash | Seismic and aseismic slip pulses driven by thermal pressurization of pore fluid[END_REF][START_REF] Rice | Stability and localization of rapid shear in fluid-saturated fault gouge: 1. linearized stability analysis[END_REF][START_REF] Platt | Stability and localization of rapid shear in fluid-saturated fault gouge: 2. localized zone width and strength evolution[END_REF][START_REF] Poulet | Thermo-poromechanics of chemically active creeping faults: 3. the role of serpentinite in episodic tremor and slip sequences, and transition to chaos[END_REF]Veveakis et al., 2014;[START_REF] Alevizos | Thermo-poro-mechanics of chemically active creeping faults. 1: Theory and steady state considerations[END_REF]Rattez, Stefanou, & Sulem, 2018;[START_REF] Rattez | The importance of thermo-hydro-mechanical couplings and microstructure to strain Resnyansky[END_REF], damage [START_REF] Lyakhovsky | A non-local visco-elastic damage model and dynamic fracturing[END_REF][START_REF] Lyakhovsky | Damage-breakage rheology model and solid-granular transition near brittle instability[END_REF], 2014a;[START_REF] Lyakhovsky | Dynamic rupture in a damage-breakage rheology model[END_REF]Kurzon et al., 2019Kurzon et al., , 2020)), granular physics (Daub et al., 2008;Daub & Carlson, 2008, 2009;[START_REF] Daub | Pulse-like, crack-like, and supershear earthquake ruptures with shear strain localization[END_REF][START_REF] Hermundstad | Energetics of strain localization in a model of seismic slip[END_REF][START_REF] Elbanna | A two-scale model for sheared fault gouge: Competition between macroscopic disorder and local viscoplasticity[END_REF] or combinations thereof. Continuum models employing a rate and state formulation with the slip rate and tractions replaced by strain rate and stress invariants have been developed by [START_REF] Herrendörfer | An invariant rate-and state-dependent friction formulation for viscoeastoplastic earthquake cycle[END_REF] and applied to self-organizing fault systems in [START_REF] Preuss | Seismic and aseismic fault growth lead to different fault orientations[END_REF][START_REF] Preuss | Characteristics of earthquake ruptures and dynamic off-fault deformation on propagating faults[END_REF][START_REF] Dalzilio | Subduction earthquake sequences in a non-linear visco-elasto-plastic megathrust[END_REF]; [START_REF] Behr | Transient slow slip characteristics of frictional-viscous subduction megathrust shear zones[END_REF]. Since the state of the subsurface is typically difficult to accurately characterize on a regional basis, these models are important to develop a generic understanding of the long-term evolution of seismogenic fault systems.

A problem commonly associated with strain softening rheologies without internal length scale is unconstrained localization (Hobbs et al., 1990). As a shear zone localizes to the scale of a mesh element or grid cell they become poorly resolved and cause mesh dependence of the simulation. Models that possess an internal length scale perpendicular to the direction of shear may not suffer from this issue, provided that this length scale is resolved by the discretization. For example, as thermo-rheological feedback mech-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth anism causes localization the heat diffusion term may gain prominence until localizing and delocalizing (e.g. diffusing) influences find a balance. Unfortunately some of the proposed fault rheologies don't possess an internal length scale, and some possess one that has physical meaning and can't be changed without changing the behavior of interest.

If such a physical length scale is very small in nature, then that may place insurmountable constraints on the mesh resolution in two-or three-dimensional regional simulations (e.g. [START_REF] Rice | Stability and localization of rapid shear in fluid-saturated fault gouge: 1. linearized stability analysis[END_REF][START_REF] Platt | Stability and localization of rapid shear in fluid-saturated fault gouge: 2. localized zone width and strength evolution[END_REF]. A few models do contain a controllable internal length scale; these are the unified rate and state friction theory of [START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF] and the damage-breakage rheology of Lyakhovsky andBen-Zion (2014b, 2014a); [START_REF] Lyakhovsky | Dynamic rupture in a damage-breakage rheology model[END_REF]; Kurzon et al. (2019Kurzon et al. ( , 2020)).

Objective

In summary, we state the need for a rigorous treatment of a continuum framework in which Dieterich-Ruina-type rate and state friction can be embedded based on the following observations:

1. Rate and state friction is a useful and powerful homogenization of the complex micro-scale processes that give rise to slip transients in nature.

2. Numerical modeling exploration of the complex and evolving relation that fault networks have with their tectonic environment is hampered by the long-standing challenges with mesh adaptive generation as well as lack of detailed knowledge of the structure of the subsurface.

3. The description of fault friction as a continuum process can be an elegant and practical way to avoid these problems provided they possess a controllable internal length scale in order to produce sensible results given reasonable computational resources.

4. A continuum faulting rheology may readily model important transient effects that might be missed in an equivalent interfacial description.

5. Much effort has been dedicated to understanding rate and state friction as an ingredient in numerical models -not all proposed continuum rheologies can equally benefit from this.

In response, we develop a Dieterich-Ruina-type continuum reformulation of rate and state friction that 1. possesses a controllable internal length scale λ 0 that enforces a constraint on localization, 2. yields mesh-independent results once λ 0 is sufficiently well resolved, 3. retains the parameters to the interfacial friction (1) and evolution (2) laws and yields similar trajectories of V and f for the same choice of parameter values, but that also 4. has interesting implications for fault zone behavior during and between earthquakes.

The model we propose here bears close relation to the earlier works of [START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF] and [START_REF] Herrendörfer | An invariant rate-and state-dependent friction formulation for viscoeastoplastic earthquake cycle[END_REF]. Contrary to [START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF] we refrain from in-depth discussion of fault physics but concentrate on the introduced spatial regularization, which is also a necessary extension of [START_REF] Herrendörfer | An invariant rate-and state-dependent friction formulation for viscoeastoplastic earthquake cycle[END_REF].

In Section 2 we will motivate our choice of constitutive model, focusing on the formal links between interfacial and continuum models. In Section 3 we will present our equivalent continuum rate and state rheology. In Section 4 we build our rheology into a numerical model of a simplified one-dimensional fault transect to demonstrate both its localization and delocalization behavior as well as its similarities and differences with respect to the original laws of Dieterich and Ruina. We close with a discussion -in which we will comment on possible links to existing theories of the physics of faulting -and a summary in Sections 5 and 6.

Material Model and Continuum Mechanics

We will first describe how the existing interfacial description of rate and state friction can be embedded in a continuum mechanical framework with internal discontinuities, and then generalize this to internal friction and distributed anelastic deformation in an interface-free model.

A Fault Reference Frame

We consider the immediate neighborhood Ω ⊂ R 3 of a section of fault Γ ⊂ Ω (Figure 1) that is described as the set of points belonging to a differentiable 2D manifold. We invoke the existence of a surjective function x f ( x) that maps a coordinate x ∈ Ω to its nearest Euclidian neighbor x f ∈ Γ. Next, we define the fault transect T ( x f ) ⊂ Ω as the 1D curve that includes all the points x ∈ Ω that map to a particular fault co-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth Ω Τ ⊂ Ω Γ ⊂ Ω φ → n → t Figure 1.
A region Ω in the enclosing a portion of a fault interface Γ, with normal and tangent vectors n and t indicated. A transect T is drawn perpendicular to the fault. The fault interface may be substituted by a shear band on which anelastic shear strain is distributed according to ϕ( x), as indicated by the colored contours.

ordinate

x f ∈ Γ, in short T = x -1 f ({ x}) ⊂ Ω. Finally, we let n( x) ∈ R 3 , x ∈ Ω be a
unit vector field locally tangent to the transecting curve T ( x f ( x)) and therefore normal to the fault.

The state of the body Ω is described by a symmetric Cauchy stress tensor field σ(t, x) ∈ R 3 ⊗R 3 and a velocity vector field v 0 v(t, x) ∈ R 3 . Their governing equations are stated after the introduction of some further scaffolding. Let t denote time.

We define the normal traction σ = -n•σ • n (positive in compression), the shear traction vector τ = σ•n+σ n, and the shear traction magnitude τ = √ τ • τ . The rupture process is most efficient when the slip rate vector V ∈ R 3 is parallel to the in-plane shear traction vector τ , and so alongside the full velocity field v we define a purely slipinduced velocity double couple field v around the coordinate x as

v (t, x) := 1 2 V (t, x f ) t( x f ) sgn(n( x f ) • [ x -x f ]), (3)
with t ∈ R 3 := τ /τ the unit tangent vector and sgn the sign function, which has sgn(0) = 0.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Assumptions

In the following, we will assume that continuum processes that occur around the fault are predominantly confined to the small neighborhood Ω, which itself includes only a small section of a whole fault. The neighborhood Ω is assumed to be sufficiently small that in-plane variations of fault properties such as state, slip rate, and curvature are negligible. This assumption comes at a loss of generality, especially around fault kinks, branches, and tips, but does not neccesarily limit the applicability of our formulation. We will revisit this assumption in the Section 5.5 of the Discussion.

We also assume that the material has uniform static properties throughout the domain. This means that the gouge layer is effectively assumed to be infinitely wide, or alternatively that the only quality that distinguishes the fault gouge from the host rock is it's 'state'. Undoubtedly there is a collection of state variables that govern the longterm structure of fault zones in nature that are not modeled here. We will revisit this assumption too in Section 5.5, but in the mean time note that its impact is largely confined to inferences made about the de-localization behavior of rate strengthening faults in Section 3

Generalization of a Fault to a Shear Band

We now generalize the fault interface to a distributed shear band by substituting a smoothed step function Φ = Φ(t, x; λ 0 ) ∈ [-1, +1], λ 0 > 0 (e.g. Φ(t, x; λ 0 ) := tanh(x/λ 0 ))

for the sign function in (3). We require that the chosen function converges in a pointwise manner to the sign function as λ 0 → 0 + , such that (3) can be seen as the result of taking the limit

v (t, x; λ 0 ) := 1 2 V (t, x f ) t( x) Φ t, n • x -x f λ 0 v (t, x) := lim λ0→0 + v (t, x; λ 0 ) . (4)
Recognizing that fault slip is an anelastic process, the anelastic strain rate tensor ε¬e is defined to be the symmetric gradient ∇ s of the (continuous) slip-induced velocity field v (t, x; λ 0 ):

ε¬e (t, x) := v 0 ∇ s v (t, x; λ 0 ) := 1 2 v 0 (∇ v ) T + (∇ v ) (t, x; λ 0 ) ≈ 1 2 v 0 λ 0 V (t, x f )ϕ t, n • x -x f λ 0 t ⊗ n + n ⊗ t ( x) =: γ 0 γ(t, x)ŝ( x). (5)

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Here the equivalent anelastic shear strain rate (a scalar value) has been introduced as

γ 0 γ(t, x) = v 0 V (t)ϕ(t, x), (6)
with reference value γ 0 = v 0 /λ 0 and strain rate distribution ϕ(t, x; λ 0 ) = ∂ x Φ(t, x; λ 0 ) (illustrated in blue shading in Figure 1). The symbol ŝ denotes the Schmidt tensor ŝ = 1 2 t ⊗ n + n ⊗ t , and its Frobenius norm ŝ = tr ŝT ŝ = 1/ √ 2, which is an important property to maintain when generalizing ŝ later on because it is compatible with the interpretation of γ 0 γ as the fault-perpendicular derivative of fault-parallel velocity, i.e. as an accurate measure of simple shear. In the derivation of (5) the gradients of V and the orthonormal bases t and n have been dropped under the assumption set out in Section 2.2 that they are small compared to the gradient of Φ. For this assumption to be met, the across-fault length scale λ 0 must be small compared to the along-fault length scales that exists in the interfacial rate and state formulation.

Since the derivative of half a smoothed step function

1 2 Φ(x) ∈ [-1 2 , + 1 2 ]
has the properties of a distribution -is a non-negative and importantly integrates to unity over its domain -we may express the integral relation

V(γ) := γ 0 T γ(t, x)dx = v 0 V (t, x f ) (7)
over the fault transect T , which is again compatible with the notion that γ 0 γ is a measure of simple shear in the fault reference frame. We note that the above integral strictly requires that all anelastic strain occurs over the length of the fault transect T , implying that the anelastic strain rate distribution has compact support in space. In practice we can relax that requirement if said distribution is narrow and thin-tailed.

Subtracting the anelastic shear strain rate tensor ε¬e from the full strain rate tensor ε = ∇ s v results in the elastic strain rate tensor. This additive decomposition of elastic and anelastic strains is known as the Maxwell model. Application of Hooke's law of linear elasticity to the elastic strain rate tensor leads to

σ = S [ ε -γ 0 γŝ] (8)
with fourth-order elastic stiffness tensor S and the dot over a symbol denoting the time derivative. The relation above is standard in the context of elasto-plasticity (e.g. [START_REF] Mühlhaus | A variational principle for gradient plasticity[END_REF], in which γ 0 γ is known as the plastic multiplier and ŝ coincides with the derivative of the magnitude of shear traction τ with respect to the full stress ten-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth sor σ. In this case τ is to be interpreted as the non-associated plastic potential of a material that does not undergo anelastic volume change.

For completeness we list ordinary differential equation that governs v, the momentum balance equation:

v 0 ˙ v -g = ρ -1 ∇ • σ. ( 9
)
Here g is the acceleration due to gravity and ρ the mass density. Gravity may be ignored in the remainder of this work at no loss of generality and the inertial term ˙ v may be ignored at some loss of generality under the assumption of static momentum balance.

Plasticity and Coordinate Invariance

For a spontaneously developing fault zone the shear and normal stress can not be defined in relation to a known plane. Plasticity models avoid this problem through the use of stress tensor invariants and scalar parameters. The Mohr-Coulomb plasticity model generates anelastic shear strain parallel to the Schmidt tensor ŝ = 1 2 (n ⊗ t + t ⊗ n), in which t makes an angle φ = tan -1 f with respect to the largest principal stress axis σ 1 , and t and n form an orthogonal triad of vectors together with the intermediate principal stress axis σ 2 . In other words, anelastic shear strain occurs in a plane perpendicular to σ 2 . The corresponding cohesionless yield criterion is written as

1 2 (σ 1 -σ 3 ) = -sin(φ) 1 2 (σ 1 + σ 3 ) , (10)
with σ 1 and σ 3 the magnitudes of the largest and smallest principal stresses σ 1 and σ 3 .

The quantity -1 2 (σ 1 + σ 3 ) = σ m is called the mean stress and 1 2 (σ 1 -σ 3 ) = σ d the effective deviatoric stress. We note that the angle of a shear band that is generated does not need to be parallel to the direction of shear strain (Vardoulakis, 1980;Vermeer, 1990;[START_REF] Marone | Coulomb constitutive laws for friction: Contrasts in frictional behavior for distributed and localized shear[END_REF]Kaus, 2010;Le Pourhiet, 2013), and moreover that a local change in friction is not expected to immediately alter the macroscopic fault angle (e.g. [START_REF] Preuss | Seismic and aseismic fault growth lead to different fault orientations[END_REF].

In the Mohr-Coulomb model, the friction coefficient f = tan(φ) of a yielding or even of a non-yielding material may be expressed from (10) as a function of the stress tensor σ as

f (σ) = σ d (σ) σ m (σ) 2 -σ d (σ) 2 -1 2 . ( 11
)

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

τ e := τ √ 2 = 1 2 tr (τ T • τ ). ( 12
)
In this model the existing definition of the Schmidt tensor ŝ is no longer applicable because the unit vectors t and n are not known a priori, and are not even uniquely constrained a posteriori. Keeping in line with plasticity theory, ŝ is redefined as the stress derivative of the non-associated plastic potential τ e :

ŝ := ∂τ e ∂σ = τ τ e .
This definition still satisfies the criterion set out in the preceding section that ŝ = 1/ √ 2.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

We calibrate the Drucker-Prager model to the Mohr-Coulomb model around a reference stress state σ 0 of simple shear combined with isotropic compression, which encompasses all stress states possible in 2D plane strain and may be considered the most relevant stress state even in 3D tectonic settings. Let σ 0 be given in ordered principal stress space by

σ 0 := -p 0       1 1 1       + τ e,0       +1 0 -1      
.

At σ = σ 0 , Mohr-Coulomb and Drucker-Prager measures of effective shear and normal stress coincide (up to a sign) and thus the Drucker-Prager model may be written as

τ e (σ) = sin(φ) p(σ). ( 13
)
Away from σ 0 (and towards a uniaxial stress state) the Drucker-Prager yield surface becomes an increasingly worse approximation to the Mohr-Coulomb yield surface (Figure 2). Higher-order approximations like the Willam-Warnke yield envelope (Ulm et al., 1999) may be considered too. The yield surfaces discussed in this section serve as part of a general model proposition -in the one-dimensional numerical experiments considered in this work only the two effective stresses of Amonton's friction law are defined.

Similar to the Mohr-Coulomb model ( 11), the friction coefficient f = tan φ of the Drucker-Prager model may also be expressed from (13) as a function of the stress tensor σ as

f (σ) = τ e (σ) p(σ) 2 -τ e (σ) 2 -1 2 . ( 14
)
In both the Mohr-Coulomb model and the Drucker-Prager model, the friction coefficient becomes ill-posed whenever the effective shear stresses σ d or τ e exceed the effective normal stresses σ m or p. We have found a practical solution to be to add a constant value to the effective normal stresses, which for steady-state friction behaves as a cohesive strength.

Despite adopting a plasticity framework to enforce the frictional yield constraint we emphasise that the resulting model is still best regarded as 'frictional viscous' rather than frictional plastic because it lacks a distinct boundary between yielding and nonyielding states in space and time -the yield constraint (an equality, not an inequality) is enforced everywhere and anytime, and at any stress.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Continuum Evolution Law

Our goal is to find a continuum state evolution law that produces nearly the same history of loading and sliding as the original formulation of Dieterich and Ruina on a discrete fault. Like V (t) is a 'global' measure of γ(t, x) over a fictitious fault transect T , we introduce a local variable θ(t, x) of which Θ(t) is a global measure. We target continuum friction and evolution laws expressed in terms of γ and θ of the form

f = f 0 + a log γ + b log θ, (15)
θ = r 0 (c 2 -c 1 M(γ) θ) , (16)
with M a non-local interaction operator that will be further elaborated in Section 3.3 and c 1 and c 2 additional coefficients that are required to calibrate ( 15) and ( 16) to their respective interfacial equivalent. The local friction and evolution laws ( 15) and ( 16) retain the structure of their interfacial counterparts ( 1) and ( 2) and generalize the continuum formulations of [START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF] and [START_REF] Herrendörfer | An invariant rate-and state-dependent friction formulation for viscoeastoplastic earthquake cycle[END_REF]. We will revisit this connection in the Section 5.4 of the Discussion. The targeted form of the continuum equations does not come out of the blue. In coming sections we hope to make clear how it arises.

Analytical Framework

By 1. reorganizing the rate and state friction law (1) into a definition of Θ(V, f ), 2. differentiating (1) with respect to time and reorganizing it into a definition of V (t), and 3. substituting the definition of Θ(V, f ) and the definition of Θ that follows from the aging law ( 2) into the definition of V (t), we arrive at the ordinary differential equation

a V = b r 0 V 2 -V Θ(V, f ) -1 + V ḟ , (17a) Θ(V, f ) -1 := V a b exp (-[f -f 0 ]/b) (17b)
in which the friction coefficient f = f (σ; t, n) is treated as a property of -rather than a constraint on -the stress tensor σ.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

By applying the same procedure to the proposed continuum friction and aging laws ( 15) and ( 16) we obtain the analogous ODE for γ:

a γ = b r 0 c 1 γM(γ) -c 2 γ θ(γ, f ) -1 + γ ḟ , (18a) θ(γ, f ) -1 := γ a b exp (-[f -f 0 ]/b), (18b)
in which f = f (σ) is given by the yield criterion ( 14) and ḟ = ḟ (σ, σ) is expressible using the same yield criterion and the elastic constitutive equation ( 8).

Equations (17a,b) and (18a,b) are useful in mathematical and numerical analysis of the problem because they eliminate a variable and an algebraic constraint, and at the same time provide an ideal reference frame for establishing, and if necessary influencing, the ability of (18a,b) to reproduce the predictions made by (17a,b) in the sense of the integral relation ( 7).

In the following derivation we will assume that the field of effective friction f ( x), x ∈ Ω in the neighborhood of a point x f on the fault is the same whether it is generated by a hard or a soft discontinuity. This assumption ties into the principal assumption set out in Section 2.2, which is revisited in Discussion section 5.5. We will also assume that M( (t)♦(x, . . .)) = (t)M(♦(x, . . .)). This assumption is satisfied in Section 3.3. Combining (17a,b), (18a,b), (6), and (7) leads to the following useful relation that describes the evolution of the anelastic shear strain rate distribution over time:

φ ∝ V λ 0 c 1 ϕ M(ϕ) -ϕ -Θ(V, f ) -1 λ a b 0 c 2 ϕ 1+ a b -ϕ . ( 19
)
This differential equation tells us several things. Most importantly, for γ to remain in accordance with V through the integral relation ( 7) over the transect T , the same integral of the left hand side of ( 19) must be zero. This requirement can only be met at abitrary (V, f ) if both the term multiplied by V and the term multiplied by Θ(V, f ) -1 integrate to zero over T . This in turn can not in general be guaranteed unless

c 1 := C 1 (ϕ) := λ -1 0 ϕ(x)M(ϕ)(x) dx -1 , (20a)
c 2 := C 2 (ϕ) := λ -a b 0 ϕ(x) 1+ a b dx -1 . ( 20b
)
We consider it undesirable for coefficients to depend in a time-dependent way on integrals of the modeled quantities and therefore will be restricting our attention to specific regimes of interest during which the values of the coefficients c 1,2 can be predicted analytically. We then employ those predicted values as model-specific constants in time and

Accepted Article

This article is protected by copyright. All rights reserved. space. We accept that this incurs a potential error whenever the state of the model is outside the selected regime, and will critically evaluate this error using numerical models in Section 4.

Assuming the coefficients c 1,2 have been chosen appropriately and for simplicity that M(ϕ) = ϕ, two additional observations can be made on the basis of ( 19):

1. The term that is multiplied by

V promotes localization if M(ϕ) = ϕ because ϕ 2
is a narrower distribution than ϕ (Figure 3).

2. The opposite (delocalization; Figure 3) holds true for the term that is multiplied by Θ(V, f ) -1 on account of the corresponding minus sign and the requirement that

a, b > 0.
This localizing and delocalizing behavior is of great interest and will be treated in detail in the following sections.

In the sections that follow we will assume to be operating on the fault-perpendicular transect T and will denote with the scalar coordinate x the distance with respect to the fault core, i.e.

x = n( x) • [ x -x f ( x)]
. This setting is in line with our assumption set out in Section 2.2 that across-fault variations in anelastic strain rate are more compact than its along-fault variations.

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Runaway Sliding and Localization

Under rate-weakening conditions a < b, earthquakes occur as superexponential solutions to (17a,b). Prior to inertial damping the seismic slip rate V behaves as

V ∝ V 2 .
This ODE has the solution

V (t) = V (0) 1 - t t * -1 , (21)
in which t * ∝ V (0) -1 is the time of the (hypothetical) singularity.

Assuming M( ) = , the same behavior occurs in (18a,b) in the limit

γ ∝ c 1 γ 2 , (22)
which is similarly solved by

γ(t, x) = γ(0, x) 1 - t t * (x) -1 , in which t * (x) ∝ (c 1 γ(0, x)) -1 .
We can make the following qualitative observations:

1. small spatial variations in initial conditions γ(0, x) can be amplified to infinity due to the fact that points with larger initial conditions on γ are closer to the strain rate asymptote at t = t * (x) than points with smaller initial conditions on γ, and 2. γ(t, x) and V (t) do not generally satisfy the integral relation ( 7) for all time t ∈ [0, t * ).

Noting that this regime of runaway sliding corresponds in (19) to the limit

φ ∝ V λ 0 c 1 ϕ M(ϕ) -ϕ ,

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth we can further illustrate the observations made above by choosing the class of solutions

(Appendix A) V (t) ∝ 1 - t t * m -1 -1 , (23a) ϕ(t, x; m) := l(t)c 1 (m) 1 + 2x l(t) m -1 , m ∈ R ≥ 2 (23b) l(t) ∝ V (t) (1-m) -1 , ( 23c
)
γ 0 γ(t, x; m) := v 0 V (t)ϕ(t, x; m), (23d)
c 1 = c 1 (m) := π m csc π m , (23e)
in which the family of strain rate distributions ϕ(t, x; m) generated by the exponent m ∈ R ≥ 2 (Figure 4c) are of generalized Cauchy type and evolve in accordance with a timevariable characteristic width l(t) (Figure 4d). We emphasize that, though illustrative, this class of solutions is by no means unique because (18a,b) contains no constraint on the spatial distribution of strain rate. Nevertheless, we note the following:

1. In the limit m → ∞, the distribution ϕ(t, x; m → ∞) reduces to a uniform distribution on the constant interval x ∈ [-l(0)/2, +l(0)/2]. The parameter c 1 (m → ∞) = 1 (Figure 4a) and the solution (23a) reduces to the interfacial solution (21).

However this limit also requires infinite mesh resolution or alternatively the definition of two mesh discontinuities, which defeats our intent.

2. In the other end-member case where m = 2, the distribution ϕ(t, x; m = 2) is known as the Cauchy distribution, and its characteristic width l(t) ∝ V (t) -1 for a given initial condition l(0) = l 0 . The parameter c 1 is chosen in such a way that the distribution ϕ(t, x; m) integrates to one and the time of the singularity matches the prediction of the interfacial model. However, the trajectory of V towards the asymptote is unavoidably different (Figure 4b).

We therefore come to the qualitative conclusion that for initially smooth and numerically resolvable strain rate distributions, orders-of-magnitude increases of slip rate readily correspond to orders-of-magnitude increase of strain rate locality, placing insurmountable demands on mesh resolution. If left unconstrained, this process culminates in finite

-19-

Accepted Article

This article is protected by copyright. All rights reserved. time to an infinite slip rate concentrated in an infinitesimally thin region of space -a plane.

Non-local Equivalent Strain Rate

In order to counter this tendency towards complete localization to a plane we introduce the non-local interaction operator M, which mollifies its operand by means of convolution with the symmetric distribution m(x):

M(γ)(t, x) := supp(m) m(ξ)γ(t, x -ξ)dξ. (24)
We note that M(γ)(t, x) still satisfies the relation (7) given that the distribution m integrates to unity over its domain.

Evaluation of M(γ) requires the evaluation of the integral of γ over the support of m(ξ) and therefore lacks the mathematical (as opposed to physical) locality that is desirable for continuum models. We follow [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF] in constructing and truncating a Taylor series linearization of (24). We find that a mollifier m(x) = 1 2λ0 exp -x λ0 (Figure 5) leads to series coefficients that are even powers of the length scale λ 0 :

M(γ)(t, x) = ∞ n=0 λ 2n 0 ∂ 2n x γ
This infinite sum can be rolled up into a recursive or implicit definition of

M = M im M -1 im (γ) := γ -λ 2 0 ∂ 2 x γ, ( 25
)
which is to be solved for the auxiliary variable γ = M im (γ). Alternatively, the terms beyond n = 1 can be dropped, leading to the explicit definition of

M = M ex M ex (γ) := γ + λ 2 0 ∂ 2 x γ. ( 26
)

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

We note that any truncation of the Taylor series into a linear combination of even derivatives of γ adheres to the integral property expressed in (7), even if these truncations lose correspondence to a concrete mollifier. The implicit definition (25) has numerically advantageous properties [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF], but requires the potentially costly solution of an elliptic PDE. For this reason we will use the explicit definition (26) in this work and discuss the implicit definition (25) as an option for future implementations.

We generalize the second partial derivatives ∂ 2

x with respect to the x coordinate that appear in ( 25) and ( 26) to Laplacians ∇ 2 in 3D, again assuming that variations of strain rate across the fault are much bigger than variations along it.

A convolution over slip history d(t) with a kernel w(d) ∝ exp -d/d c , which is identical to the kernel discussed above, was also employed by Ruina (1980, eq. 15, 16a, 17a) and Ruina (1983, eq. 13a,b) as a 'generator' of evolution laws. There, it represented the 'fading memory' of slip history on the fault. Here, we have consistently extended that notion with a non-local interaction of deformation that fades with distance away from a point in the shear band. Although not primarily intended to have physical meaning in this work, we note that the Laplacian of effective strain rate plays a role in granular physics (e.g. [START_REF] Bouzid | Nonlocal rheology of granular flows across yield conditions[END_REF][START_REF] Bouzid | Microrheology to probe non-local effects in dense granular flows[END_REF]. There, the length scale λ 0 is proportional to the grain size [START_REF] Mühlhaus | The thickness of shear bands in granular materials[END_REF][START_REF] Bouzid | Nonlocal rheology of granular flows across yield conditions[END_REF][START_REF] Bouzid | Microrheology to probe non-local effects in dense granular flows[END_REF]Rice, 2006, and ref-erences in the latter).

Non-locality as a Localization Limiter

We repeat the analysis of Section 3.2, but now use the explicit non-local operator

M ex in γ ∝ c 1 γ M ex (γ) and φ ∝ V λ 0 c 1 ϕ M ex (ϕ) -ϕ . ( 27
)

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

As shown in Appendix C, we find that it is solved by

V (t) ∝ 1 - t t * -1 , (28a) ϕ(x) :=          1 π λ0 cos 2 1 2 x λ0 ∀ x ∈ [-πλ 0 , +πλ 0 ] 0 ∀ x / ∈ [-πλ 0 , +πλ 0 ] , (28b)
γ 0 γ(t, x) := v 0 V (t)ϕ(x), (28c)
c 1 := 2π, (28d)
and observe that 1. The trajectory of V towards its asymptote can be made to exactly match the prediction ( 21) of the interfacial rate and state friction formulation (17a,b) by choos-

ing c 1 = 2π.
2. The cosine-squared distribution with prescribed size λ 0 is an attractive steady state during the process of earthquake nucleation. Further localization will not spontaneously occur. We interpret this distribution as striking a balance between the localizing tendency of "ϕ 2 -ϕ" and the delocalizing tendency of "ϕ∇ 2 ϕ -ϕ"

(Equation ( 19) and Figure 3).

The implicit non-local operator M im in the same context does not have analytical solutions that simultaneously satisfy V ∝ V 2 and φ = 0 with non-degenerate distribution ϕ = ϕ(x/λ 0 ).

Steady-state Friction and Deformation

As mentioned in Section 1.1, an important property of the state evolution law is the steady state Θ V = 1. We wish to retain this steady state in the continuum equivalent (16) of the state evolution law, ideally at a steady strain rate distribution ϕ(x), cf.

(19). We are thus interested in a steady-state solution to

φ ∝ λ 0 c 1 ϕ M ex (ϕ) -λ a b 0 c 2 ϕ 1+ a b , (29)

Accepted Article

This article is protected by copyright. All rights reserved. which exists as long as a, b > 0 (standard requirements) and a < b (a net weakening 572 fault). As detailed in Appendix B, this solution is given by

573 ϕ ss (x) =          ϕ 0 cos 1 2 x λ1 2[1-a b ] -1 ∀ x ∈ [-πλ 1 , +πλ 1 ] 0 ∀ x / ∈ [-πλ 1 , +πλ 1 ] (30a) 574 λ 1 := λ 0 1 - a b -1 . ( 30b
)
575 576

Faults that are net strengthening or neutral (a > b) only possess a uniform steady state 577 strain rate distribution that is bounded by the (fictitious) walls of the gouge layer (Sec-578 tions 2.2 and 5.5).

579

-23-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

The coefficients c 1 and c 2 , and the newly introduced ϕ 0 , are given by

c 1 = 1 2 π -1 2 Γ 1 + 2 1 -a b -1 Γ 1 2 + 2 1 -a b -1 [λ 1 ϕ 0 ] -2 (31a) c 2 = 1 2 π -1 2 Γ 2 1 -a b -1 Γ 1 2 3 + a b 1 -a b -1 [λ 1 ϕ 0 ] -1-a b (31b) ϕ 0 = 1 2 π -1 2 Γ 1 + 1 -a b -1 Γ 1 2 + 1 -a b -1 λ -1 1 , (31c)
with Γ the gamma function that has the property Γ(n) = (n -1)!∀n ∈ N + . The coefficients c 1 and c 2 are plotted as functions of a/b in Figure 7. The solutions (30a) are plotted in Figure 6, in which we can see that steady-state anelastic strain rate distributions have finite width and are therefore numerically resolvable as long as the aforementioned requirements on a, b, and a/b are met.

Damping

Finite time blow-up of the solution (Section 3.5) under quasi-static loading conditions is understood to be the consequence of the model's inability to radiate or dissipate liberated potential strain energy away from the fault. It appears impossible to directly counteract this problem by implementing inertia with absorbing boundary conditions in the 1D models we present here, due to the resulting model's tendency to quickly establish a global steady state at an unphysically high strain rate. Therefore, inspired by the radiation damping approximation of [START_REF] Rice | Spatio-temporal complexity of slip on a fault[END_REF] (see also [START_REF] Ben-Zion | Slip patterns and earthquake populations along different classes of faults in elastic solids[END_REF][START_REF] Rice | Slip complexity in earthquake fault models[END_REF]Lapusta et al., 2000;[START_REF] Thomas | Quasi-dynamic versus fully dynamic simulations of earthquakes and aseismic slip with and without enhanced coseismic weakening[END_REF] and the Kelvin-Voigt visco-plastic regularization of [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF] (see also F. [START_REF] Wu | Deformation trapping due to thermoplastic instability in one-dimensional wave propagation[END_REF][START_REF] Peirce | Material rate dependence and localized deformation in crystalline solids[END_REF]Duretz et al., 2019[START_REF] Duretz | Toward robust and predictive geodynamic modeling: The way forward in frictional plasticity[END_REF][START_REF] Duretz | Modeling lithospheric deformation using a compressible visco-elasto-viscoplastic rheology and the effective viscosity approach[END_REF]de Borst & Duretz, 2020;[START_REF] Stathas | The role of viscous regularization in dynamical problems, strain localization and mesh dependency[END_REF], we add a linear viscous damping term ηγ in a parallel configuration to the continuum friction law:

f = f 0 + a log γ + b log θ + ηγ.
This article is protected by copyright. All rights reserved.

a ζ(t, x) = b r 0 c 1 M ex • γ • ζ(t, x) -c 2 θ(ζ, f ) -1 + ḟ (t), (33a) θ(ζ, f ) -1 := (γ • ζ(t, x)) a b exp -1 b [f -f 0 -η γ • ζ(t, x)] , (33b)
γ(ζ) := a η W 0 η a exp ζ , (33c)
with the principal branch W 0 (♦) of the Lambert W function, which is an implicit transcendental function that is defined to be the solution to the equation W 0 exp W 0 = ♦.

While we acknowledge that the use of implicit functions in physical descriptions is not ideal, from a practical standpoint there is little problem because reliable and fast algoritms for computing the common Lambert W function are included in many programming languages [START_REF] Fritsch | Algorithm 443: Solution of the transcendental equation we w = x[END_REF]Barry, Barry, & Culligan-Hensley, 1995;Barry, Culligan-Hensley, & Barry, 1995;[START_REF] Johansson | Computing the lambert w function in arbitrary-precision complex interval arithmetic[END_REF].

As indicated in Figure 8, the function γ(ζ) is approximately linear at high ζ (and correspondingly high γ). This means that in this limit, (33) simplifies to

η γ(t, x) = b r 0 c 1 M ex • γ(t, x) -c 2 θ(γ, f ) -1 + ḟ (t), (34a) θ(γ, f ) -1 := exp -1 b [f -f 0 -η γ] , (34b)
We can now compare reaction terms of ( 18) and (34). Ignoring spatial derivatives and the temporal derivative ḟ of friction in both, the reaction term of (18) that is active at

-25-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth low strain rate looks like γ ∝ γ 2 -α(f )γ 1+ a b while the reaction term of (34) that is active at high strain rate looks like γ ∝ γ-α(f ) exp( η b γ). The former is a convex function of γ, the latter is a concave function of γ. We can thus see that the same conditions (large f , small α(f )) that readily trigger runaway sliding in (18) also introduce a stable steady state at high γ in (34) and therefore also in (33). We note that this quality of having a stable steady state at non-zero γ means that (33) can be classified as a reactiondiffusion equation of generalized Fisher or KolmogorovPetrovskyPiskunov (KPP) type [START_REF] Fisher | The wave of advance of advantageous genes[END_REF]Kolmogorov et al., 1937). This type of equation is known to exhibit outward propagating (strain rate) fronts, and consequently does not admit a steady strain rate distribution. This markedly limits the validity of the analytical values derived for c 1,2 during the coseismic phase. We will revisit this issue in detail throughout the results section.

In the interfacial representation of a fault in a homogeneous half-space, the shear traction τ (t, z) at a time t and location z along the fault may be decomposed into

f (t, z)σ = τ (t, z) = τ 0 (t, z) + τ d (t, z) -η * V (t, z),
where τ 0 (t, z) is the loading stress, τ d (t, z) is the convolution that expresses the dynamic stress transferred to the point z at time t by slip at all points within the causality cone around (t, z) [START_REF] Rice | Spatio-temporal complexity of slip on a fault[END_REF][START_REF] Cochard | Dynamic faulting under rate-dependent friction[END_REF], 1996;Lapusta et al., 2000). The term η * V exists to compensate for the removal of a singularity in aforementioned convolution, and this sets the viscosity η * = 1 2 µ c equal to half the material's shear impedance, with c s = µ/ρ the shear wave speed and in turn µ the shear modulus and ρ the mass density. In the quasi-dynamic or radiation damping approximation [START_REF] Rice | Spatio-temporal complexity of slip on a fault[END_REF], the dynamic stress transfer term τ d (t, z) is replaced by an easier to compute static stress transfer term τ s (t, z), but the damping term η * V is retained. Here we take the traction τ 0 + τ s to be the one produced by our static momentum balance equation ( 9) and elastic constitutive equation ( 8), and implicitly move over the damping term to the right-hand side of the equation. Calibration of (32) to the result gives our non-dimensional damping viscosity as

η = c 3 1 2 µ σ v 0 c s ,
with c 3 ∼ 1 an additional dimensionless calibration constant that is chosen somewhat arbitrarily to be equal to c 1 , which seems to give maximum macroscopic slip rates that are similar to those produced by the interfacial friction law with radiation damping.

Analytical predictions made in preceding sections are complemented here with numerical simulation results that better illustrate the complex time-dependent behavior of the system of equations, and allow us to establish the consequences of approximations and other model choices made in the process of formulating a local equivalent continuum rate and state friction law. As before we restrict the scope to 1D models across the shear band.

Equations

On the one-dimensional line x ∈ [-L/2, +L/2] that deforms under simple shear, analogous to the transect T defined in Figure 1, we model the compact equation ( 33a,b)

together with the one-dimensional static momentum balance equation ∂τ /∂x = 0 cf.

(9), elastic constitutive equation τ (t, x) = µ [∂v(t, x)/∂x -2γ 0 γ(t, x)] cf. (8) (with µ the shear modulus as before), yield equality τ (t) = f (t)σ, and finally the boundary con-

ditions v(t, x = ∓L/2) = ∓V p /2 combined into the single ODE ḟ (t) = µ σ V p L - V(γ • ζ) L , (33d)
with V p the 'plate' driving velocity that is applied antisymmetrically on both ends of the domain. This equation uses (7) that relates slip velocity to the integral (from -L/2 to +L/2) of anelastic shear strain rate γ, as well as (33c) for the relation between γ and the auxiliary field ζ.

The system is closed with natural boundary conditions ∂γ/∂x = 0 on γ(t, x = ∓L/2), and initial conditions f (t = 0) = f 0 + ∆f and γ(t = 0, x) ∝ ϕ ss (x) + , with ϕ ss (x) deriving from the analytical prediction (30a) and a necessary but small homogeneous background value of anelastic shear strain rate. As long as this value is small enough ( V p /L v 0 /L) and the 'tails' of the strain rate distribution can be ignored, the required computational domain size may be based on the length scale λ 0 and therefore cover only a potentially small fraction of the line [-L/2, +L/2], greatly reducing the computational cost.

The solution procedure is outlined in Appendix D.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved. -29-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Continuum Velocity Stepping Friction Experiments

Velocity stepping laboratory experiments form the foundation on which rate and state friction was proposed [START_REF] Dieterich | Time-dependent friction and the mechanics of stick-slip[END_REF]. We compare numerical velocity stepping experiments of a traditional rate and state friction-governed interface with our proposed continuum equivalent under the same conditions. We impose a driving velocity that follows a smooth square wave function with respect to sliding distance d:

log 10 (V p (d)/v 0 ) = 2 π tan -1 - 1 ν cos π d n d c , ( 36
)
with n = 20 the amount of slip weakening distances d c between velocity steps, and ν = 10 -3 a parameter that controls the smoothness of the smooth square wave, with small numbers giving the squarest result (see Figure 9a). Some smoothness is important to maintain numerically stable results. The low and high driving velocities are 10 -1 v 0 and 10 +1 v 0 respectively. Other model parameters are given in Table 4.2. Inertial effects may be ignored (both in the original lab experiments and in our simulations) due to the low driving velocity. Thus, the quasi-radiation damping viscosity η is set to zero and its effects not discussed here. In the limit η → 0 + , (33c Because the domain size L is very small and the system (33a-33d) therefore stiff, the friction tends to a stable steady state after being perturbed by a velocity step rather than developing a limit cycle. This also means that V(γ) ≈ V p (d). In fact, V(γ) and V p (d) are so close that their difference would not register in Figure 9a.

)
In the following, we use the terms time dependence and slip dependence interchangeably but note that the rate of change of either is not constant from the perspective of the other. A slip-centric presentation of velocity-stepping results is commonplace in the earlier cited literature.

In Figure 9b (and its detailed view 9e) we explore the similarities and differences in slip dependence of the friction coefficient f between the interfacial and continuum velocity stepping experiments. Based on (19), we expect that a velocity step perturbs not just the magnitude of anelastic shear strain, but also its distribution (Figure 9d) and therefore the theoretical value of the dimensionless coefficients c 1,2 (Section 3.3; Equations (20a) and (20b); Figure 9c, solid lines) away from steady state. Since these coefficients are actually set to their constant steady-state predictions (Section 3.5; Equations (31a)

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth the stress drop is already about one third of the total. The same limit cycle shown in Figure 12f in log-linear phase space (U , log 10 V ) allows us to see that there are subtle differences between the limit cycles of the interfacial and continuum models over the whole range of slip velocities.

A further result of practical interest is the difference in adaptively chosen time step between the interfacial and continuum models. 12d,g,k). From these same plots it can be seen that, just like stress begins to drop in significant advance of an event, so too does the time step of the continuum model.

Finally, we remark that spatial resolution tests indicate that the quality of the solution is only influenced by the ratio of cell size h x to λ 0 , with acceptable results achieved when h x /λ 0 ≥ 10. For the results presented here we have used h x /λ 0 = 20. Changing λ 0 in proportion to h x does nothing except to cause a wider or narrower but equally well resolved strain distribution and a virtually identical limit cycle.

Discussion

So far we have argued for a visco-elastic continuum rheology that resembles both the mathematical presentation as well as the resulting behavior of traditional interfacial rate and state friction as described by Dieterich and Ruina in their respective seminal puplications. Before addressing this work in the context of a much broader body of existing research and reflecting on critical assumptions made in this work and its extension towards the future, we briefly touch upon some alternate branches of our proposed theory that have been left out so far for the sake of clarity.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Figure 12. Panels a-g: behavior of the interfacial (bright yellow lines) and continuum (dark blue lines) models for a full limit cycle. Panels h-k: focus on the 100 seconds around an event.

From left to right, horizontal axes are: time t in years, slip deficit U in meters (linearly related to shear stress τ and friction coefficient f ), and time t in seconds. Vertical axes from top to bottom are slip deficit U in meters, slip velocity V in meters per second on a linear scale, the same on a logarithmic scale, and finally the time step size ht in seconds on a logarithmic scale.

Where limit cycles are shown (middle column, panels e-g), the cycle sense is indicated with small arrows.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

A Primitive Reformulation

In contrast to traditional models of plasticity, the continuum rate and state friction laws presented here possess a continuous time dependence that makes the anelastic shear strain rate γ a predictable quantity rather than one that is to be solved by a constrained optimization algorithm (e.g. [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF][START_REF] Duretz | The benefits of using a consistent tangent operator for viscoelastoplastic computations in geodynamics[END_REF]. This is apparent from the ODE (18a,b) for γ, but is hidden in the Dieterich-Ruina form ( 15) and ( 16) of the continuum rate and state friction rheology. However, (18a,b) lacks the simplicity and elegance of the Dieterich-Ruina system and is difficult to interpret as the resultant of physical phenomena. By defining a set of primitive variables κ and ψ as κ = r κ γ (37)

ψ = r ψ γ a/b exp (-[f (σ) -f 0 ]/b),
with the respective reference rates r κ = r 0 c 1 and r ψ = r 0 c 2 , we are able to write a more expressive formulation of (18a,b) as

κ = r κ γ = r κ ψ/r ψ b a exp ([f (σ) -f 0 ]/a) (38) ψ = r ψ exp (M ex (κ) -ψ) , (39)
that still closes the momentum balance equation ( 9) and elasto-plastic constitutive equation (8) without algebraic constraints and thus without requiring a constrained optimization algorithm. Equation (39) can optionally accommodate damping using the Lambert W function analogous to the way it was discussed in Section 3.6.

We recognize κ(t, x) as a dimensionless measure of accumulated anelastic shear strain and therefore assign to ψ the meaning of a 'virtual' or 'unrealized' strain. In spite of the conceptual appeal of ( 38) and ( 39), it can be difficult to choose appropriate initial conditions for ψ.

Quadratic Diffusion

Instead of the mollified term γ M(γ) that was introduced in ( 16) and used throughout the remainder of this work, we could have chosen to use M(γ 2 ):

a γ = b r 0 c 1 M(γ 2 ) -c 2 γ θ(γ, f ) -1 + γ ḟ . ( 40
)
and the definition of c 1 is adapted to the new way of mollifying according to the same principles as laid out in Section 3.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

This choice has advantages and disadvantages. A major disadvantage is that it does not seem to allow a return from (40) to the more elegant Dieterich-Ruina form ( 15) and ( 16), nor to the compact form ( 38)-( 39) proposed in the preceding section. A clear advantage is that, taken together with the explicit non-local operator M ex , (40) produces a degenerate reaction-diffusion equation of generalized Fisher or KolmogorovPetrovskyPiskunov (KPP) type [START_REF] Fisher | The wave of advance of advantageous genes[END_REF]Kolmogorov et al., 1937) that remains well-posed at and around γ = 0. This leads to shear zone solutions that can propagate into perfectly intact rock without taking recourse to an arbitrary small initial condition on γ. Furthermore, under quadratic diffusion, the function γ(ζ) used in the ODE form (33a-d) in Section 3.6 becomes algebraic instead of transcendentally implicit.

We have reproduced all results presented in Section 4 with this quadratic diffusion term and zero background anelastic strain rate too, and observe numerically smoother and better resolved shear fronts during the coseismic delocalization phase. These fronts also travel faster and further than in the linear diffusion case leading to a larger discrepancy between measured and used values of the coefficients c 1,2 and consequently larger deviations of the limit cycle with respect to that of the interfacial rate and state friction model. The generalized Fisher-KPP type equation has merited a large body of theoretical research (e.g. [START_REF] Sánchez-Garduño | An approximation to a sharp type solution of a density-dependent reaction-diffusion equation[END_REF], 1995[START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF][START_REF] Sherratt | Nonsharp travelling wave fronts in the fisher equation with degenerate nonlinear diffusion[END_REF][START_REF] Petrovskii | An exactly solvable model of population dynamics with density-dependent migrations and the allee effect[END_REF][START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] Gilding | A fisher/kpp-type equation with densitydependent diffusion and convection: travelling-wave solutions[END_REF]; Y. [START_REF] Wu | Stability of travelling waves with algebraic decay for n-degree fisher-type equations[END_REF]Y. Li & Wu, 2008;[START_REF] Malaguti | Asymptotic speed of propagation for fishertype degenerate reaction-diffusion-convection equations[END_REF][START_REF] Broadbridge | Exact solutions for logistic reaction-diffusion equations in biology[END_REF]) and therefore some of the above mentioned empirical findings may be given a theoretical underpinning in future work.

Relation to Regularized Damage or Plasticity Models

As shown in Section 4.3 there are circumstances under which the transient viscous rheology proposed in this work promotes a spontaneous organization of periods and regions of negligible anelastic strain rate and those of significant anelastic strain rate. Even though a critical yield stress seems to be an emergent rather than an inherent property of the system of equations and initial and boundary conditions, it is tempting to think of this rheology as a smooth plasticity model.

We have shown (Section 3.5) that our model requires regularization by a non-local strain rate measure to avoid spurious mesh dependence. Non-locality in one form or an-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth other has been applied to combat ill-posedness in e.g. [START_REF] Bažant | Continuum theory for strain-softening[END_REF][START_REF] Triantafyllidis | A gradient approach to localization of Research[END_REF]; [START_REF] Schreyer | One-Dimensional Softening With Localization[END_REF]; Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF];

de [START_REF] De Borst | Gradient-dependent plasticity: Formulation and algorithmic aspects[END_REF]; [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF]; [START_REF] Jirásek | Nonlocal models for damage and fracture: comparison of approaches[END_REF]; [START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: survey of progress[END_REF]; [START_REF] Engelen | Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour[END_REF]; [START_REF] Jirásek | Nonlocal models for damage and fracture: comparison of approaches[END_REF]Rolshoven (2009b, 2009a) The comprehensive non-local damage-breakage rheology of [START_REF] Lyakhovsky | A non-local visco-elastic damage model and dynamic fracturing[END_REF]; Lyakhovsky andBen-Zion (2014b, 2014a); [START_REF] Lyakhovsky | Dynamic rupture in a damage-breakage rheology model[END_REF], firmly rooted in thermodynamic theory and well-calibrated to match observations, is one of few damage theories developed to cover the full process of earthquake generation and healing [START_REF] Lyakhovsky | Dynamic rupture in a damage-breakage rheology model[END_REF]Kurzon et al., 2019Kurzon et al., , 2020)). As in our model, the nonlocality in theirs is of Gradient type, but differs subtly in that it acts on the damage parameter rather than on the anelastic shear strain rate. Another example is the Godunov-Peshkov-Romenski model (S. [START_REF] Godunov | Nonstationary equations of nonlinear elasticity theory in eulerian coordinates[END_REF][START_REF] Godunov | Elements of continuum mechanics and conservation laws[END_REF]Resnyansky et al., 2003;[START_REF] Romenskii | Deformation model for brittle materials and the structure of failure waves[END_REF][START_REF] Romenski | Modeling wavefields in saturated elastic porous media based on thermodynamically compatible system theory for two-phase solid-fluid mixtures[END_REF], which was first used in [START_REF] Gabriel | A unified first-order hyperbolic model for nonlinear dynamic rupture processes in diffuse fracture zones[END_REF] to simulate dynamic rupture and off-fault damage generation. This model differs from our model and the aforementioned damage-breakage rheology in a fundamental way in that shear bands are produced by damage waves described by hyperbolic equations rather than parabolic equations of reaction-diffusion type.

Our model differs from phase-field models of fracture in some respects. We have used the distribution ϕ of anelastic shear strain rate as an analytical tool during the derivation of the continuum rheology proposed in this work. While this distribution could be renormalized to form a phase field d ∈ (0, 1), the necessary inclusion of a damping viscosity in Section 3.6 causes temporal variations in the anelastic shear strain rate distribution that transcend the modeled distributions, and make renormalization to a phase field possible at best as a post-processing step in a simulation. We recognize that the inclusion of a damping viscosity constitutes the use of a 'double-well potential' in phase field terminology, which multiple authors have commented on critically in this context (e.g. Kuhn et al., 2015;[START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF].

Our linear viscous damping appears as a regularization technique in some plasticity models [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] Peirce | Material rate dependence and localized deformation in crystalline solids[END_REF][START_REF] Wu | Deformation trapping due to thermoplastic instability in one-dimensional wave propagation[END_REF]Duretz et al., 2019[START_REF] Duretz | Toward robust and predictive geodynamic modeling: The way forward in frictional plasticity[END_REF][START_REF] Duretz | Modeling lithospheric deformation using a compressible visco-elasto-viscoplastic rheology and the effective viscosity approach[END_REF]de Borst & Duretz, 2020;[START_REF] Stathas | The role of viscous regularization in dynamical problems, strain localization and mesh dependency[END_REF]. In these mod-

-39-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth els the Kelvin-Voigt arrangement of yield strength and Newtonian viscosity truncates the steady-state anelastic strain rate that may be achieved, forcing a shear zone to have a finite width in order to slide at a certain macroscopic rate. In our model, the same viscosity also introduces an effective upper bound on anelastic strain rate attained during the limit cycle and causes subsequent delocalization, but we must still rely on the gradient regularization to combat the unbounded localization that would otherwise happen even at low strain rate. A detailed analysis of the efficacy of this viscous regularization in dynamical problems is made in [START_REF] Stathas | The role of viscous regularization in dynamical problems, strain localization and mesh dependency[END_REF].

Relation to Other Transient Continuum Rheologies

The framework here proposed can be seen as a generalization of the work of Herrendörfer ically and effectively yielding a numerical method analogous to the stress glut method of [START_REF] Andrews | Test of two methods for faulting in finite-difference calculations[END_REF]. Herrendörfer's model was subsequently applied in an unconstrained evolving continuum model in [START_REF] Preuss | Seismic and aseismic fault growth lead to different fault orientations[END_REF][START_REF] Preuss | Characteristics of earthquake ruptures and dynamic off-fault deformation on propagating faults[END_REF], but notwithstanding measures put in place that acknowledge the changing distribution of shear strain rate within a shear zone, their model ultimately lacks regularizations that remove mesh dependence. In the chapters titled 'Localization of Deformation' and 'Relationship of Localization to Instability' of his PhD thesis, [START_REF] Ruina | Friction laws and instabilities: A quasistatic analysis of some dry friction behaviour[END_REF] gives a thoughtful take on aspects of the localization behavior of a strain rate formulation of rate and state friction (without spatial regularization), which is in some aspects in line with findings reported in this work, and complementary in others.

As noted in Section 1.2 of the Introduction, a variety of continuum theories have been developed to explain the general rate-and state-like behavior of deformation in faults and shear zones (e.g. [START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF][START_REF] Braeck | Spontaneous thermal runaway as an ultimate failure mechanism of materials[END_REF]Daub et al., 2008;Daub & Carlson, 2008, 2009;[START_REF] Daub | Pulse-like, crack-like, and supershear earthquake ruptures with shear strain localization[END_REF][START_REF] Elbanna | A two-scale model for sheared fault gouge: Competition between macroscopic disorder and local viscoplasticity[END_REF][START_REF] John | Generation of intermediate-depth earthquakes by selflocalizing thermal runaway[END_REF][START_REF] Thielmann | Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating[END_REF][START_REF] Thielmann | Grain size assisted thermal runaway as a nucleation mechanism for continental mantle earthquakes: Impact of complex rheologies[END_REF][START_REF] Rozel | A thermodynamically self-consistent damage equation for grain size evolution during dynamic recrystallization[END_REF][START_REF] Barbot | Modulation of fault strength during the seismic cycle by grainsize evolution around contact junctions[END_REF][START_REF] Niemeijer | A microphysical model for strong velocity weakening in phyllosilicate-bearing fault gouges[END_REF]

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth [START_REF] Romenskii | Deformation model for brittle materials and the structure of failure waves[END_REF][START_REF] Chen | Rate and state frictional and healing behavior of carbonate fault gouge explained using microphysical model[END_REF]Van den Ende et al., 2018;[START_REF] Roubek | A note about the rate-and-state-dependent friction model in a thermodynamic framework of the Biot-type equation[END_REF][START_REF] Lyakhovsky | A non-local visco-elastic damage model and dynamic fracturing[END_REF][START_REF] Lyakhovsky | Damage-breakage rheology model and solid-granular transition near brittle instability[END_REF], 2014a;[START_REF] Lyakhovsky | Dynamic rupture in a damage-breakage rheology model[END_REF]Kurzon et al., 2019Kurzon et al., , 2020;;[START_REF] Poulet | Thermo-poromechanics of chemically active creeping faults: 3. the role of serpentinite in episodic tremor and slip sequences, and transition to chaos[END_REF]Veveakis et al., 2014;[START_REF] Alevizos | Thermo-poro-mechanics of chemically active creeping faults. 1: Theory and steady state considerations[END_REF]Rattez, Stefanou, & Sulem, 2018;[START_REF] Rattez | The importance of thermo-hydro-mechanical couplings and microstructure to strain Resnyansky[END_REF][START_REF] Pozzi | Coseismic fault lubrication by viscous deformation[END_REF]. We note again that an internal length scale in these models typically arises from the inclusion of a diffusion process (of e.g. temperature, pore pressure) but that may not always be adjusted to meet the constraints imposed by scale and computational power without changing the outcome of the model. In our proposed formulation an artificial diffusion process acts directly on the anelastic shear strain rate, resulting in a robust and controllable internal length scale.

Spiers

In this last respect our work is preceded by nearly 25 years by [START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF]. In this work Sleep combined and extended earlier works [START_REF] Sleep | Ductile creep, compaction, and rate and state dependent friction within major fault zones[END_REF][START_REF] Segall | Dilatancy, compaction, and slip instability of a fluid-infiltrated fault[END_REF][START_REF] Linker | Effects of variable normal stress on rock friction: Observations and constitutive equations[END_REF][START_REF] Chester | Effects of temperature on friction: Constitutive equations and experiments with quartz gouge[END_REF][START_REF] Chester | A rheologic model for wet crust applied to strike-slip faults[END_REF] in which rate and state friction was interpreted as the product of crack generation and healing, associated rheological weakening, and dissipative heating. This physical reasoning resulted in a model that contains only quantities that are either directly measurable or can be modeled by independent methods.

This contrasts with our purely mathematical argumentation that serves to retain close correspondence to the original phenomenological description of rate and state dependent sliding on a frictional interface. [START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF] neglects the fluxes and associated spatial gradients of the pore fluid but does include heat diffusion, which they note does however not play a significant role at the scale of their numerical experiments. Instead, resembling our approach, they impose an artificial length scale and forced strain distribution by explicit mollification of the anelastic shear strain rate with a Gaussian kernel.

This leads to an aging law (Equation 53in [START_REF] Sleep | Application of a unified rate and state friction theory to the mechanics of fault zones with strain localization[END_REF]) that is structurally identical to our result (16). Like us, they find that strain localization can only occur when a < b, and that a rate-strengthening effect that activates at high strain rate leads to strain delocalization. With respect to Sleep's valuable contribution, in this work we provide a more complete argumentation for this type of spatial regularization and analysis of the resulting patterns of strain localization and delocalization over the seismic cycle.

Assumptions and Future Work

We proposed our model in a very general three-dimensional continuum mechanics framework, but for simplicity have considered only a small fault neighborhood in which

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth in-plane variations of fault properties can be neglected (Section 2.2) so that the model becomes effectively one-dimensional. This same assumption was also taken in the numerical models that we have used in our analysis. The assumption is clearly violated around fault branches, at fault kinks or on rough faults, and near the fault or rupture tip. We note however that rate and state friction was proposed based on laboratory studies that also neglect these geometrical complexities. Although rate and state faults with branches and kinks are still largely non-standard in present-day numerical modeling studies, plenty of attention has been given to the critical nucleation patch and the structure of the cohesive zone near the rupture tip (e.g. [START_REF] Rice | Spatio-temporal complexity of slip on a fault[END_REF][START_REF] Cocco | On the slip-weakening behavior of rate-and state dependent constitutive laws[END_REF][START_REF] Cocco | Physical interpretation of the breakdown process using a rate-and state-dependent friction law[END_REF][START_REF] Rubin | Earthquake nucleation on (aging) rate and state faults[END_REF][START_REF] Day | Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture[END_REF]Lapusta & Liu, 2009;[START_REF] Viesca | Stable and unstable development of an interfacial sliding instability[END_REF]Viesca, , 2016a;;[START_REF] Putelat | A phase-plane analysis of localized frictional waves[END_REF]. Given our main assumption, these features can only be accurately reproduced with our continuum formulation in higher-dimensional numerical models if the regularization length scale λ 0 is significantly smaller than the length scales associated with the critical nucleation patch and the cohesive zone. In turn the cell size must be sufficient to resolve λ 0 , and so we expect to need a grid resolution that is significantly higher than that of existing methods to simulate rate and state frictional interfaces. We have also seen that the diffusion process is associated with a smaller time scale, and thus, stricter time step constraints than the interfacial model. We note that both spatial and temporal resolution requirements already place challenging constraints on simulations of seismic and aseismic slip sequences [START_REF] Erickson | The Community Code Verification Exercise for Simulating Sequences of Earthquakes and Aseismic Slip (SEAS)[END_REF][START_REF] Jiang | Community-driven code comparisons for threedimensional dynamic modeling of sequences of earthquakes and aseismic slip[END_REF].

It seems prudent to first make a more detailed assessment of the computational demands and the techniques that may be required to meet the resolution requirements (e.g. adaptive mesh refinement, local adaptive time stepping). In the process of constructing higher-dimensional models one may first concentrate on the friendliest regions of parameter space, for example those that promote stable sliding, or undamped runaway localization as in [START_REF] Viesca | Stable and unstable development of an interfacial sliding instability[END_REF]Viesca ( , 2016a[START_REF] Viesca | On the existence of a nucleation length for dynamic shear rupture[END_REF].

The temporal patterns of localization and delocalization that occur in our model yield testable predictions that tie in to a recent surge in interest in similar patterns observed in the lab and in nature (e.g. [START_REF] Ben-Zion | Localization and coalescence of seismicity before large earthquakes[END_REF][START_REF] Mcbeck | Investigating the onset of strain localization within anisotropic shale using digital volume correlation of time-resolved x-ray microtomography images[END_REF][START_REF] Mcbeck | Deformation precursors to catastrophic failure in rocks[END_REF]McBeck, Ben-Zion, & Renard, 2020;[START_REF] Mcbeck | The competition between fracture nucleation, propagation, and coalescence in dry and water-saturated crystalline rock[END_REF]. It could be the scope of future research to reinterpret anelastic strain rate in our model as a measure of the activity of a statistical distribution of cracks of various properties and com--42-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth pare to aforementioned lab and field observations. In this context, the delocalization that is in our models induced by a rate-limiting Kelvin viscosity is reminiscent of the growth of off-fault fracture networks during dynamic rupture (e.g. [START_REF] Templeton | Off-fault plasticity and earthquake rupture dynamics: 1. dry materials or neglect of fluid pressure changes[END_REF][START_REF] Okubo | Dynamics, radiation, and overall energy budget of earthquake rupture with coseismic off-fault damage[END_REF][START_REF] Gabriel | A unified first-order hyperbolic model for nonlinear dynamic rupture processes in diffuse fracture zones[END_REF].

We have assumed an infinite gauge layer and observed that materials that are ratestrengthening at steady state (a > b) can only feature delocalization of strain. This appears inconsistent with the widespread observation of faults and fault gauges that are more or less stable (e.g. [START_REF] Ikari | On the relation between fault strength and frictional stability[END_REF][START_REF] Carpenter | Frictional properties and sliding stability of the san andreas fault from deep drill core[END_REF][START_REF] Coble | In situ frictional properties of San Andreas Fault gouge at SAFOD[END_REF][START_REF] Bedford | The stabilizing effect of high pore-fluid pressure along subduction megathrust faults: Evidence from friction experiments on accretionary sediments from the nankai trough[END_REF]. A possible explanation is that a gauge layer is initially formed by a set of strain-softening processes, subsequently chemically, petrologically, texturally or geometrically matured over time to become rate-strengthening, but generally weaker than host rock. The gauge layer walls may then act as a barrier to further delocalization. While this situation could be simulated in our approach by varying the frictional properties accordingly, the consistent way to model the long term evolution of faults likely requires at least two internal state variables. The damage-breakage rheology of Lyakhovsky andBen-Zion (2014b, 2014a) has this feature.

In this work, we have restricted ourselves to classical Dieterich-Ruina rate and state friction with aging law. In Section 1.2, we have hinted at a sizable number of physical interpretations of the phenomenological rate and state friction problem. This concerns notably the behavior of weakening and the interpretation given to the 'state' variable.

At the same time, a number of alternative phenomenological evolution laws have been proposed, among which the slip law (Dieterich, 1979a). It would be worthwhile to investigate the possibility of applying the mathematical techniques developed here to a wider range of friction and evolution laws, especially those involving multiple states, such as for example temperature, pore pressure, or grain size. However, this might not be straightforward because we have relied heavily on analytical solutions to calibrate the interfacial and continuum descriptions of friction, and finding these solutions is frequently a time-consuming task with uncertain prospects for success.

Conclusions

In this work we have carefully constructed a coordinate-invariant and mesh-independent transient visco-elastic continuum rheology that behaves in a way that is consistent with fusion that ultimately spurred us to pursue this work. We also gratefully acknowledge fruitful discussions with Robert Herrendörfer, Robert Viesca, Vladimir Lyakhovsky, Dmitry Garagash, Sylvain Barbot, and René de Borst. We greatly appreciate the thoughtful re- ing data are available in the supplementary materials. We depend on the MIT-licenced 'Scientific Colour Maps' package [START_REF] Crameri | Scientific colour maps[END_REF] for distortion-free representation of the model results -also for readers with color vision deficiencies [START_REF] Crameri | The misuse of colour in science communication[END_REF].

We then take the time derivative of the Ansatz (A2) and equate the result to (A1), yielding γ(t, x) = γ(t, x) (A3)

At x = 0, we have:

Γ(t) = Γ(t) 2 . (A4)
The solution to (A4) is

Γ(t) = Γ 0 (1 -Γ 0 t) -1 , Γ 0 = Γ(t = 0). ( A5
0 c 2 ϕ 1+ a b , (B1)
subject to a, b > 0 and a < b, with M ex given by ( 26) and the coefficients c 1 and c 2 measured from the sought solution ϕ(∞, x) by means of (20a) and (20b). We summarize the PDE that is to be solved as

ϕ 2 + ϕλ 2 0 ∂ 2 x ϕ -c ϕ 1+ a b = 0, ( B2

Accepted Article

This article is protected by copyright. All rights reserved. (B6)

We arbitrarily take the distribution ϕ(∞, x) to be symmetric around x = 0, where it reaches its maximum value ϕ 0 = ϕ(∞, 0). At this point, as well as at the extremities of the distribution where ϕ = 0, we set the gradient, i.e. f (ϕ), to zero. This eliminates the term f (0) 2 in (B6) and constrains the free parameter c to

c = 1 2 1 + a b ϕ 1- a b 0 . ( B7
)
We now multiply both sides of (B6) by 2 and take the square root to obtain

f (ϕ) = ϕ 1- a b 0 ϕ 1+ a b -ϕ 2 . (B8)
Applying the inverse function theorem to f (ϕ

) -1 = 1 -a b λ -1 1 ϕ (x) -1 = 1 -a b λ 1 -1 x (ϕ),
we obtain from (B8)

x (ϕ)

λ 1 = 1 -a b ϕ 1- a b 0 ϕ 1+ a b -ϕ 2 - 1 2
.

(B9)

We seek to integrate this relation once more over the region ε ∈ [0, ϕ) to obtain a solution for x(ϕ). The integral of the right hand side of (B8) can be reverse-engineerded from the known derivatives

cos -1 (f ) = -(1 -f 2 ) -1 2 .
and d f f α = αf α-1 , noting also that ϕ/ϕ 0 ≤ 1, to give

1 -a b ϕ 1- a b 0 ϕ 1+ a b -ϕ 2 - 1 2 = ∂ ϕ   -2 cos -1   ϕ ϕ 0 1 2 [1-a b ]     . (B10)
This yields

x(ϕ) λ 1 - x(0) λ 1 = -2 cos -1   ϕ ϕ 0 1 2 [1-a b ]   + π. (B11)
Setting x(0) λ1 = -π and inverting for ϕ(∞, x) gives ϕ(∞, x) = ϕ 0 cos 1 2

x λ 1 2[1- a b ] -1 . (B12) ϕ(∞, x) =          ϕ 0 cos 1 2 x λ1 2[1- a b ] -1 ∀ x ∈ [-πλ 1 , +πλ 1 ] 0 ∀ x / ∈ [-πλ 1 , +πλ 1 ]. (B13)
Finally, we determine the coefficient ϕ 0 that ensures that the distribution ϕ integrates to one over its domain, and the coefficients c 1 and c 2 following the expressions (20a) and (20b). Definite integrals of (B13) are evaluated using Wolfram Mathematica (Wolfram Research, Inc., 2017), and we find

ϕ 0 = 1 2 π -1 2 Γ 1 + 1 -a b -1 Γ 1 2 + 1 -a b -1 λ -1 1 (B14a) c 1 = 1 2 π -1 2 Γ 1 + 2 1 -a b -1 Γ 1 2 + 2 1 -a b -1 [λ 1 ϕ 0 ] -2 (B14b) c 2 = 1 2 π -1 2 Γ 2 1 -a b -1 Γ 1 2 3 + a b 1 -a b -1 [λ 1 ϕ 0 ] -1-a b . (B14c)
We have verified that (B13) is a solution to (B2) given (B7) using automated symbolic manipulation in Wolfram Mathematica.

Accepted Article

This article is protected by copyright. All rights reserved. where tildes indicate approximate space-discrete quantities and numerical arrays are indicated in bold face. We form a symbolic rules for computing the Jacobian matrix J of this system of non-linear space-discrete equations,

J( ζ, f ) =    ∇ ζ g( ζ, f ) ∇ f g( ζ, f ) ∇ ζ h( ζ, f ) ∇ f h( ζ, f )    ,
   + ∆ t    g( ζk , f k ) h( ζk , f k )    (D1) f =    ζk+1 f k+1    -    ζk f k    -∆ t    g( ζk+1 , f k+1 ) h( ζk+1 , f k+1 )    = 0. (D2)
Our algorithm makes an explicit-in-time prediction using (D1) and evaluates the l 2 norm of the implicit-in-time residual f . It includes the possibility to perform Newton-Raphson iterations using the Jacobian J to keep the residual bounded, although we find it to be more efficient to experimentally set the dimensionless time step ∆ t /S t sufficiently small (e.g. 0.1) to never cause tolerances to be violated. This value is then taken to correspond to a stable time step size of the explicit problem.

Figure 2 .

 2 Figure 2. Comparison of Mohr-Coulomb (thick yellow lines) and Drucker-Prager yield envelopes (thin blue lines) at increments of pressure. The field of view is the octahedral plane of principal stress space. Projections of principal stress axes are indicated with arrows, and projections of planes of simple shear are indicated with dashed lines. A friction coefficient of 0.3 was used to generate this figure.

  [START_REF] Herrendörfer | An invariant rate-and state-dependent friction formulation for viscoeastoplastic earthquake cycle[END_REF] used the Drucker-Prager model[START_REF] Drucker | Soil mechanics and plastic analysis or limit design[END_REF] as a simple and smooth approximation to Mohr-Coulomb plasticity. The model is defined in terms of the straightforwardly computable invariants 1. pressure p = -1 3 tr σ, and 2. effective shear stress τ e proportional to the Frobenius norm of the deviatoric stress tensor τ := σ + p δ, with δ the Kronecker delta:

  Figure 3. Influence of the different terms in Eq. (19) on the time derivative φ under the assumption that M(ϕ) = ϕ. The first term (blue) can be seen to promote localization, while the second term (orange) favors the opposite. Also shown in this figure in green is the delocalizing influence of the nonlinear Laplacian term ϕ(x)ϕ (x), which is introduced in Sections 3.3 and 3.4. The input function ϕ follows a cosine-squared distribution.

Figure 4 .

 4 Figure 4. Asymptotic solutions during earthquake nucleation. a) Value of the dimensionless coefficient c1(m) versus distribution exponent m (see main text) represented on a reciprocal axis that ranges between 2 and ∞. The thick line acts also as a color bar approximately indicating the value of m at which the individual lines in the remaining figures are plotted. b) Dimensionless slip rate V versus dimensionless time t/t * to the slip rate asymptote for different values m ∈ 2 1.0 , 2 1.5 , 2 2.0 . . . . c) Nondimensionalized distributions ϕ(x/l(t); m)/ϕ(0; m) versus dimensionless coordinate x/l(t). d) Dynamic dimensionless length scale l(t/t * )/l0 versus dimensionless time t/t * .

Figure 5 .

 5 Figure 5. Plot of the mollifier m(x) = 1 2λ exp -x λ .

  Figure 6. Steady-state distributions ϕ [cf. Eq. (30a)] of strain rate at various a/b versus distance x with respect to the fault core. The curve plotted for a/b = 0 also corresponds to the most extreme strain rate localization at any value of a/b (see main text).

Figure 7 .

 7 Figure 7. The coefficients c1 a b and c2 a b , given by Eq. (31a) and (31b), are plotted versus the acceptable range of a b .

Figure 8 .

 8 Figure 8. The function γ(ζ) defined in Eq. (33c) plotted on a logarithmic scale in the top left diagonal portion of the figure and on a linear scale in the opposing portion, highlighting the exponential behavior of γ(ζ) at low ζ and the linear behavior at high ζ.

Figure 9 .

 9 Figure 9. Results of the numerical velocity stepping experiments: a) dimensionless sliding velocity V closely tracking a smooth square wave signal (Eq. (36)); b) frictional response of the interfacial reference model (thick yellow line) and continuum model (dark blue line) with for reference the theoretical steady-state friction fss(V ) (black line); c) measured (solid lines; cf. Eq. (20a,b)) vs. used (dashed lines; cf. Eq. (31a,b)) values of the dimensionless coefficients c1,2; d) distribution ϕ = γ/V of anelastic shear strain rate with darker blue colors reflecting very low values and bright yellow colors reflecting high values, and contours distributed evenly on a linear scale. Reflecting the ultimate slip dependence (rather than direct time dependence) of the rate and state friction laws, curves in panels a)-d) are plotted against sliding distance d (measured in critical slip distances dc) on the vertical axis. Panel e) shows an enhanced view of the evolution of the instantaneous friction coefficient towards the steady-state value with slip distance d, simultaneously showing a positive and negative step for both the interface and continuum models. Panel f) displays the same, but artificially increases the critical slip weakening distance of the interfacial friction law by 20%.

  Some rate and state parameter values are roughly based on Erickson et al. (2020). a Using (31a) and (31b) with stated values of a and b.

  simply reduces to γ(ζ) = exp ζ or equivalently ζ(γ) = ln γ, which reduces (33a,b) back to (18a,b).

  Figure 10. Evolution of the continuum seismic cycle simulation demonstrated by a succession of six events plotted against time step count on the horizontal axis, with a seventh event stretched by a factor four in order to better show some features of an individual event. The end of each event is punctuated by a slip velocity minimum and numbered 1 to 7 on the horizontal axes. Each revolution around the limit cycle takes 83.00 years and 107,500 time steps to complete. Panel a: time step size ht (thin dark line; left axis; logarithmic scale) and slip velocity V (thick bright line, right axis; logarithmic scale). The former serves to aid the interpretation of the time-dependent progression of each event while the latter gives a sense of magnitude to the color scale of panel b. Axes are chosen so that together these curves express the imperfect reciprocal relation between time step and slip velocity. Panel b: logarithm of dimensionless anelastic shear strain rate versus signed distance x measured in number of length scales λ0 away from the shear zone center. Bright yellow colors indicate high strain rate; dark blue colors indicate low strain rate. Contours are drawn at equal intervals in log space. No explicit color scale is given because the magnitude of anelastic shear strain rate γ is tied to the value of λ0 in order to produce a length-scale independent history of sliding velocity V (panel a). Panel c: measured (solid lines; cf. Eq. (20a,b)) vs. used (dashed lines; cf. Eq. (31a,b)) values of the dimensionless parameters c1,2.

  Figure 10a shows an approximate inverse relation between the time step of continuum model to the aggregate slip velocity v 0 V , which is unsurprising seeing how d c /[v 0 V (t)] gives a sensible local time scale for both interfacial and continuum rate and state friction models. The approximate inverse relation breaks down at large time scales where the corresponding time steps are truncated to a fixed value of 0.1 yr, and at large slip velocity and small time step, where the time scale of the continuum model becomes dominated by large spatial gradients and is generally much smaller than that of the interfacial model (Figures

  ; Burghardt et al. (2012); Lyakhovsky et al. (2011); Lyakhovsky and Ben-Zion (2014b, 2014a); Lyakhovsky et al. (2016); Kurzon et al. (2019, 2020); Kiefer et al. (2018); Abdallah et al. (2020).

  et al. (2018);[START_REF] Preuss | Seismic and aseismic fault growth lead to different fault orientations[END_REF][START_REF] Preuss | Characteristics of earthquake ruptures and dynamic off-fault deformation on propagating faults[END_REF], who made the purely local substitution v 0 V → h x γ 0 γ and used a Drucker-Prager elastoplastic model similar to the one set out in Section 2.4. Setting λ 0 , η = 0 and assuming that strain rate fully localizes into a discrete Dirac function sampled every h x , we find that the coefficients c 1 and c 2 become h x /d c and 1 respectively, substitution of which into (16) yields Herrendörfer's version of the aging law. Their model was first applied along a predefined staggered grid line in[START_REF] Herrendörfer | An invariant rate-and state-dependent friction formulation for viscoeastoplastic earthquake cycle[END_REF], obeying the discrete Dirac distribution of anelastic shear strain automat-

  views of Norman Sleep and two anonymous reviewers, which encouraged us to improve the work. This research has been supported by the European Union's Horizon 2020 Research and Innovation Programme under the ERC StG TEAR, grant no. 852992. CP acknowledges additional financial support by Swiss National Science Foundation grant 200021-169880. PS acknowledges financial support from the Swiss University Conference and the Swiss Council of Federal Institutes of Technology through the Platform for Advanced Scientific Computing (PASC) program. A.-A.G. acknowledges additional support by the European Union's Horizon 2020 Research and Innovation Programme under ChEESE, grant no. 823844 and the German Research Foundation (DFG projects grants no. GA 2465/2-1, GA 2465/3-1). None of the authors declare any competing interests or financial conflicts. All scripts used to generate the figures and their underly-

1 ,

 1 ) with c a parameter to be determined later. Dividing the whole equation by ϕ, and mak-ing the substitution λ 0 → 1 -a b λ x ϕ = 1 -a b λ 1 ∂ x f (ϕ) = 1 -a b λ 1 f (ϕ)∂ x ϕ = f (ϕ)f (ϕ). (B4)-46-

  manuscript submitted to JGR: Solid Earth by exploiting symmetry across the shear zone. Integrals are evaluated using a midpoint rule. This gives the discrete systemζ = g( ζ, f ) ḟ = h( ζ, f ),

  with the upper left block a dense N ×N matrix with a dominant sparse band structure (N being the problem size), the lower right block a 1×1 empty matrix, and the off-diagonal blocks densely populated vectors of compatible shape. The system is then linearized as d time scale S t ( ζ, f ) to be resolved is computed as the reciprocal of the largest eigenvalue of the Jacobian matrix J. The actual time step ∆ t is computed as some fixed fraction of S t , optionally bounded by a maximum value and/or maximum growth rate to prevent time step overestimation as the time scale increases. The equations are discretized in time with forward and backward Euler schemes, both first-order accurate,

Table 1 .

 1 Parameter values used in the velocity stepping experiment

	manuscript submitted to JGR: Solid Earth

Informally, the action of this term may be understood to truncate the anelastic shear strain rate, which then causes the distribution of said strain to widen under continued loading. We formalize this claim by repeating the derivation in Section 3.1 of the anelastic shear strain rate ODE (18a,b), now incorporating the damping term in (32). We arrive at the following result, written in terms of the newly introduced auxiliary field ζ:

-24-

Accepted Article

Acknowledgments

We are grateful to Yehuda Ben-Zion for his recommendation to look into damage dif-

Data Availability Statement

All codes and algorithms to generate and visualize the results discussed in this work can be found as Wolfram Mathematica (Wolfram Research, Inc., 2017) notebooks in the supplementary information.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth and (31b); dashed lines in Figure 9c), we expect some differences in time or slip dependence to occur. Notably, a positive velocity step leads to a temporary localization of anelastic deformation and a corresponding decrease in the theoretical values of c 1,2 with respect to steady state. This leads to an overestimation of corresponding terms in the righthand-side of (33a) and a faster decay towards the steady-state friction coefficient (Figure 9b,e) of the continuum model with respect to the interfacial model. A negative velocity step causes temporary delocalization (Figure 9d) and under-estimation of terms involving c 1,2 in (33a), but does not lead to any noticeable difference in friction coefficient evolution towards the steady state between the two models (Figure 9b,e). We conclude that those terms involving c 1,2 are negligible in this situation.

In Figure 9c and 9d it can also be seen that the distribution of anelastic shear strain rate tends more slowly to a steady state than the friction coefficient itself, and moreover that this decay is slower after a positive velocity step than after a negative one. In fact, the 40 d c wavelength we use for the intput signal ( 36) is too small to enable a somewhat complete evolution towards steady-state of the anelastic strain rate distribition during the high-velocity regime, but we have taken care that this does not meaningfully impact the results during the subsequent low-velocity regime.

Most of the difference between the interfacial and continuum velocity stepping results disappears if the critical slip weakening distance d c is seen as an additional calibration parameter. Continuing the analogy to the classical experimental setup, when we regard the continuum formulation proposed here to be the 'generator' of experimental data, and the result of the existing interfacial rate and state friction and aging laws as a data fitting curve, we observe a critical slip weakening distance d c that is about 20% larger than the value that was used to generate the data. The result is plotted in Fig- ure 9f. There it can be seen that the discrepancies remaining after the calibration of d c are small compared to the expected noise level of experimental results.

Continuum Spring-slider Experiments

Whereas velocity stepping experiments are useful for studying the way in which friction can be attracted to a steady state, permanently out-of-equilibrium 'limit cycle' behavior is more interesting in the study of earthquake dynamics. We achieve this behavior by increasing the size of the domain and the distance of the boundary conditions

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth by many orders of magnitude (L = 80 km) with respect to the velocity stepping experiments. Thereby, we greatly reduce the effective stiffness of the medium as can be seen in (33d). This is analogous to performing the so-called spring-slider experiment in which a mass is dragged over a surface by a spring that is tensed at a constant rate to give an educational example of the apparent stick-slip behavior of seismogenic faults. Our goal in performing this experiment is to determine the degree to which the here proposed continuum friction laws match the predictions made by Dieterich and Ruina's interfacial laws and to describe the transient behavior of anelastic strain rate in the added dimension.

Damping plays an essential role in these experiments to close the limit cycle that otherwise extends to V → ∞ (Section 3.6). As discussed in the same section, our implementation of damping as a linear bulk viscosity precludes a steady-state anelastic strain rate distribution at high slip velocity, a fact that directly contradicts one of our precepts (Section 3.1). Therefore we anticipate a larger discrepancy between interface and continuum models here than was observed during the velocity stepping experiments.

Figure 10 depicts the evolution of quantities of interest over the course of multiple orbits of the solution as it converges to the limit cycle. The first important observation is that the continuum model still appears to exhibit spontaneous limit cycle behavior. Due to the appropriately chosen initial condition, the wind-up period is short and there is little variability between successive revolutions. As predicted (Section 3.5), the nucleation phase is marked by runaway amplitude increase of a quasi-steady strain rate distribution (panel b). As before, we measure the dynamic value of the coefficients c 1,2 from the strain rate distribution following (20a,b), and compare to the analytically derived steady-state values (31a,b). In line with analytical predictions this phase is also accompanied by measured values of c 1,2 (panel c) that are slightly lower than the constant values that are prescribed (dashed lines in panel c). The opposite situation occurs throughout the remainder of each cycle where anelastic strain is more widely distributed.

We can attribute this directly to the effect of damping at high strain rate, which exerts a delocalizing influence. Figure 11 provides a more intuitive view of the short-lived delocalizing behavior of a strain pulse using linear rather than logarithmic scales. We emphasize that the ultimate extent of coseismic delocalization compared to the length scale λ 0 is dependent on model and material parameters, among which the effective stiffness of the medium.

Accepted Article

This article is protected by copyright. All rights reserved. Large peaks in the measured values of c 1,2 shown in panel c of Figure 10 occur during the post-and interseismic phase, where anelastic shear strain rate is broadly and quite uniformly distributed. The values that are actually used (dashed lines) are way too small here, however, terms that do not involve these coefficients dominate the time-dependent behavior in this regime and the discrepancy's consequences are limited. This is largely confirmed in Figure 12, which tabulates the evolving scalar properties of a limit cycle of the interfacial and continuum models of rate and state friction. In particular, the duration of the limit cycle (first column, panels a-c, time in years on the horizontal axis) is very similar even though the timing of the event in relation to its enclosing slip velocity minima is somewhat different.

Markedly different are the coseismic aspects of the limit cycle. Panel e of the middle column and panels h-j of the right column of Figure 12 clearly show this. The continuum model has a faster stress drop (panel h, U ∝ f ) and a more symmetric slip rate response with time compared to the interfacial model. The limit cycle depicted in linear phase space (U , V ) in Figure 12e shows a skewed triangular trajectory of the interfacial model and a more parabolic trajectory for the continuum model. It appears that the amount of stress drop that happens before peak slip rate in the interfacial model is small -about 10% of the total stress drop. At the same point in the continuum model

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth rate and state friction on an interface. We have shown that inclusion of a diffusion-like spatial regularization ensures a limit to strain localization and thus guarantees mesh convergence. In a simplified 1D fault transect, important metrics of the seismic cycle -such as slip rate and friction -are independent of the diffusion length scale associated with the regularization. However, throughout this work we have assumed the regularization length scale to be small compared to the length scales associated with other features of interest, such as fault curvature or along-strike variations of slip rate and stress. Therefore, going forward, high-resolution 2D or 3D numerical models are required to apply this model to the study of the seismogenic behavior of emerging and evolving fault zone networks. Our continuum rheology resembles a reaction-diffusion equation for anelastic strain rate. Processes described by such equations are ubiquitous in nature, and it is tempting to compare temporal patterns of localization and delocalization produced by our model with natural observations.

Accepted Article

This article is protected by copyright. All rights reserved. In this appendix we derive a generic family of solutions to the runaway slip regime discussed in Section 3.2. We consider the ODE γ(t, x) = γ(t, x) 2 , (A1) in which t and γ are taken to be non-dimensionalized in such a way that the rate constant equals one and is thus dropped. We emphasize that since no derivatives with respect to x are included, the ultimate collective behavior of γ(t, x) is not uniquely defined.

Even so, our Ansatz is a product-wise decomposition of the solution into a purely timedependent term Γ and a self-similar contribution g:

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth then we obtain

and

The function g(t, x/l(t)) integrates to

and thus the distribution ϕ(t, x) corresponding to g(t, x/l(t)) is given by

Finally, we use ( 7) and (A11) to write (A5), (A8), and (A9) in terms of V (t), and obtain:

Appendix B Derivation of the Steady-state Strain Rate Distribution

Here we derive the steady-state solution presented in Section 3.4 to the time-dependent

Appendix C Derivation of the Runaway Strain Rate Distribution

We simply note that the time-depedent PDE

to which a steady-state solution is sought in Section 3.5, closely relates to (B1) if in that equation the substitution a b → 0 is made. Making the same substitution in the result of Appendix B, we obtain:

Appendix D Numerical Solution Procedure

We discretize the Laplacian using a standard second-order accurate central difference stencil, with natural boundary conditions of the same accuracy implemented by staggering the fields with respect to the physical domain walls. The problem size is halved

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

We make our algorithms and scripts to produce figures available as Wolfram Mathematica (Wolfram Research, Inc., 2017) notebooks in the Supplementary Materials.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth Faulkner, D., Mitchell, T., Jensen, E., & Cembrano, J. ( 2011). Scaling of fault damage zones with displacement and the implications for fault growth processes.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Accepted Article

This article is protected by copyright. All rights reserved.