Qualitative study of a geodynamical rate-and-state model for elastoplastic shear flows in crustal faults

The Dieterich-Ruina rate-and-state friction model is transferred to a bulk variant and the state variable (aging) influencing the dissipation mechanism is here combined also with a damage influencing standardly the elastic response. As the aging has a separate dynamics, the overall model does not have a standard variational structure. A one-dimensional model is investigated as far as the steady-state existence, localization of the ataclastic core, and its time response, too. Computational experiments with a damage-free variant show stick-slip behavior (i.e. seismic cycles of tectonic faults) as well as stable slip under very large velocities.


Introduction
In the last decades the mathematical interest in geophysical problems was steadily growing.While there is already a large body of work in atmospheric and oceanographic fluid flows, the mathematics for geophysical models for solid earth is much less developed.The latter concerns in particular the deformation and motion of lithospheric plates in the upper crust, in particular earthquakes.The difficulties in these models is the complex behavior of rock that behaves elastically like a solid in the case of seismic waves on short time scales but behaves like a viscoplastic fluid when considered over centuries.However, very slow motion of long periods are crucial for building up internal stresses that are then released in short rupture events triggering earthquakes.Only recently, a new class of periodic motions in the Earth crust was detected by evaluating GPS measurements, namely the so-called "episodic tremor and slip" (cf.[3,28]): Here all motions are so slow that no seismic waves are emitted, but there exist two distinct regimes, one involving inelastic motions and one involving slow smooth slip.These events are observed in so-called subduction zones and have periods in the range of a few years while the overall shear velocity rate is in the range of millimeter per year.
In addition to these temporal time scales there are also several spatial scales involved.For instance, between tectonic plates there form weak regions called faults that are relatively narrow but may accumulate relatively large deformations, in particular in rapid shearing events.We refer to [26,34,38,40,42] for some recent efforts in geodynamical modeling towards a better understanding of these phenomena.On the mathematical side the work started less than a decade ago and is still comparably small, see [22,23,25,37,39,46].Moreover, there is a dichotomy with respect to bulk interface models, where most of the nonlinear effects are localized in the interface (e.g. by a so-called rate-and-state dependent friction law), and pure bulk models where typically only existence results for solutions are obtained but no qualitative behavior of the solutions can be deduced.
With this work we want to initiate a mathematical study where pure bulk models are considered but still interesting qualitative features can be deduced.In this first study we will confine ourselves to a simplified "stratified" setting where only shear deformations are considered that depend on a one-dimensional variable x ∈ (−H, H) representing the transverse direction to a straight fault or damage zone between two compact rocks representing two plates that move with respect to each other, see Figure 2.1.The continuum model is given in terms of • the shear velocity v = v(t, x) ∈ R, • the elastic strain ε = ε(t, x), • the plastic strain p = p(t, x), • the internal damage variable α = α(t, x), and • the internal aging variable θ = θ(t, x).
The model to be studied in its simplest form is the following system of five partial differential equations posed for (t, x) ∈ (0, ∞) × (−H, H) (see (2.12) for the more general case treated below): ∂ .
p as well as on the aging variable θ.Thus, we are able to mimic the commonly used Dieterich-Ruina rate-and-state friction law [21,49] where now the aging variable can be interpreted as the "state" while the dependence on π = .
p gives the rate dependence.
p = 0 and Sign( . Our paper is organized as follows: In Section 2 we provide the background from geodynamics introducing the rate-and-state friction models with a given interface and our distributed-parameter model which is slightly more general than (1.1).In particular, Section 2.2 discusses the steady-state equation where .v = .
p is independent of time.The full evolutionary model is then introduced in Section 2.3.
The analysis of steady states is the content of Section 3. In Theorem 3.1 we provide an existence theorem for steady states under quite natural assumptions and arbitrary shear velocities v(±H) = ±v ∞ .The proof relies on a Schauder fix-point argument and we cannot infer uniqueness, which is probably false in this general setting.In Proposition 3.4 we show that for steady states the limit η → 0 + in (1.1b) can be performed in such a way that accumulation points are still steady states.
In Section 4 we discuss the full dynamic model, show its thermodynamic consistency, and derive the natural a priori estimates.For our main existence result we restrict to the case without damage, i.e.C is independent of α and α ≡ 1 solves (1.1b).The result of Theorem 4.1 is obtained by time discretization and a staggered incremental scheme mimicking the solution of the static problem in Theorem 3.1.The analytical aspects are nontrivial because of the non-variational character of the problem, the non-polynomial friction law (2.4) leading to usage of Orlicz spaces, and the lack of compactness for the elastoplastic wave equation.
The final Section 5 is devoted to a numerical exploration of some simplified models that show the typical behavior expected also for the full model.The simplified model is obtained from (1.1) by neglecting α as in Section 4 and by further ignoring inertia (i.e.setting = 0 and choosing η = 0), see Section 5.1: In Section 5.2 we discuss the steady states (θ stst , π stst ) where π stst = Π(σ stst , θ stst ).We do a parameter study for varying κ and v ∞ and obtain a monotone behavior with respect to v ∞ , namely θ stst is decreasing and π stst is increasing.We always observe spatial localization in the sense that π stst is supported on The pure existence of steady states does not say anything about stability in the dynamic model (1.2).In Section 5.3 we provide a two-dimensional ODE model where there is a unique steady state that is unstable for small positive v ∞ and convergence of general solutions to periodic motions.Similarly, Section 5.4 shows simulations for system (1.2) which shows convergence towards (θ stst , π stst ) if v ∞ is large but predicts convergence towards time-periodic solutions that also have a clearly defined plastic zone smaller than (−H, H), see Figures 5.6 and 5.7.
A surprising effect is that the width 2h of the core of the fault (the active cataclastic zone) does not tend to be 0 if the plasticity gradient is ignored by setting η = 0, and even not if the aging gradient is ignored by setting κ = 0.In Proposition 3.5 we show that under natural assumptions on the rate-and-state friction law one obtains a linear dependence h = h * (v ∞ , 0) = |v ∞ |/π * for shear velocities with |v ∞ | < Hπ * , where π * is uniquely determined by the friction law and the aging law.
Another noteworthy effect is that the length scale of the aging qualitatively influences the character of response, varying in between the stick-slip and the sliding regimes.In particular, for very large shear velocities v ∞ (which are not relevant in usual geophysical faults in the lithosphere) the fault goes into a continuous sliding mode and no earthquakes occur.Actually, this is a recognized attribute of this friction model which in [4] has been compared to the observation of our "everyday life when one often manages to get rid of door-squeaking by a fast opening".In contrast under very slow shear velocities, the friction threshold is not reached for large time spans after a relaxation.Only when enough shear stress has build up, the threshold can be overcome.But then not only stresses are released but also the aging variable is reduced which leads to a much larger stress release than needed.Hence, another long waiting time is needed until next "earthquake" will start.
2 Setup of the geodynamical model

Geodynamical background
Earth's crust (together with lithosphere) is a rather solid rock bulk surrounding the lower, more viscous parts of the planet.It is subjected by damage typically along thin, usually flat weak surfaces, called faults, which exist within millions of years.The faults may exhibit slow sliding (so-called aseismic slip) or fast rupture (causing tectonic earthquakes and emitting seismic waves) followed by long period or reconstruction (healing) in between particular earthquakes.The former phenomenon needs some extra creep-type rheology modeled using a plastic strain variable or some smoothing of the activated character of the frictional resistance at very small rates (cf.Remark 3.3) and will not be scrutinized in this article, while the latter phenomenon needs some friction-type rheology.Thus faults can be modeled as frictional contact surfaces or as flat narrow stripes.
As for the frictional contact, the original Dieterich-Ruina rate-and-state friction model [21,49] prescribes the tangential stress σ t on the frictional interface as where the normal stress σ n is considered to be given (= a so-called Tresca friction model) and v is (the norm of) the tangential velocity jump along interface.The (given) parameters a and b are the direct-effect and the evolution friction parameters, respectively, d c is the characteristic slip memory length, and v ref reference velocity.If a−b > 0, we speak about velocity strengthening while, if a−b < 0, we speak about velocity weakening -the latter case may lead to instabilities and is used for earthquake modeling.The friction coefficient µ = µ(v, θ) depends in this model on the velocity magnitude v and an internal variable θ being interpreted as an aging variable, sometimes also as damage.The evolution of θ is governed by a specific flow rule typically of the form of an ordinary differential equation at each spot of the fault, say: .
with some continuous nonnegative functions f 0 and f 1 More specifically, f 0 (θ) = 1 and f 1 (θ) = θ/d c with d c > 0 is most common, considered e.g. in [6-8, 14, 16, 27, 41, 50]; then for the static case v = 0, the aging variable θ grows linearly in time and has indeed the meaning of an "age" as a time elapsed from the time when the fault ruptured in the past.
The steady state cf.[36], and then θ stays bounded and asymptotically approaches θ ∞ in the steady state . This suggests to interpret θ rather as a certain hardening or "gradual locking" of the fault in the "calm" steady state v = 0.An obvious undesired attribute of (2.1) is, as already noted in [21, p.108], that, "as v or θ approach zero, eqn.(2.1) yields unacceptably small (or negative) values of sliding resistance" µ.Therefore, (2.1) obviously violates the Clausius-Duhem entropy inequality, although being used in dozens of geophysical articles relying that in specific applications the solutions might not slide into these physically wrong regimes.Nevertheless, a regularization leading to µ > 0 and thus to a physically correct non-negative dissipation is used, too, typically as [20], cf.e.g. also [36]: In what follows, we will therefore have in mind rather (2.4) than (2.1).For an analysis and numerics of the rate-and-state friction in the multidimensional visco-elastic context we refer to [35,[37][38][39].
Since the velocity occurs in the aging flow rule (2.2), this nonisothermal friction model however does not seem consistent with standard thermodynamics as pointed out in [44] in the sense that the evolution (2.2) does not come from any free energy.On top of it, it has been known from the beginning of this rate-and-state model that it does not fit well some experiments [48] and (rather speculative) modifications e.g. by using several aging variables (which naturally opens a space for fitting more experiments) have been devised, cf.[49].
A rather formal attempt to overcome the mentioned thermodynamical inconsistency has been done in [39] by introducing two energy potentials.Thermodynamically consistent models have been devised either by using isothermal damage with healing [46] or by nonisothermal damage when temperature variation was interpreted approximately as a sliding velocity magnitude v.The latter option uses the idea that the slip of the lithospheric fault generates heat which increases temperature on the fault.In geophysical literature, the heat produced during frictional sliding is believed "to produce significant changes in temperature, thus the change of strength of faults during seismic slip will be a function of ... also temperature", cf.[10, p.7260].The usage of an (effective) interfacial temperature discussed in [14,16] following ideas from [31].In [5,10,11,51] the classical rate-and-state friction law is also made temperature dependent.Experimentally, even melting of rocks due to frictional heating is sometimes observed.
A simplified friction model µ(v) = µ 0 + aln(b|v|+1) or µ(v) = µ 0 + (a−b)ln|v/v ref | is sometimes also considered under the name rate-dependent friction [19,32,49,52] and was analyzed in [33] as far as its stability.In contrast, the above mentioned variant of temperature dependent friction can be called purely state dependent.
The friction model is sometimes "translated" into a bulk model involving a plastic-like strain and the sliding-friction coefficient µ then occurs as a threshold (a so-called yield stress) in the plastic flow rule, cf.[44,Sect. 6], or [15,16,26,32,52], known also under the name a shear-transformation-zone (STZ) concept referring to a (usually narrow) region in an amorphous solid that undergoes plastification when the material is under a big mechanical load.Instead of velocity dependence (2.4), one should play with dependence on the strain rate, cf.(2.7) below.These options can be "translated" into the bulk model by making the yield stress µ dependent, beside the strain rate, also on an aging variable θ, or on an temperature, or on a damage, or on various combination of those.Altogether, one thus get a wide menagerie of friction-type models.
Here we consider, as rather standard in geophysical modeling as (2.4), an isothermal variant and make µ dependent on strain rate and on aging.We consider also damage (or phase-field) as usual in fracture mechanics to illustrate its a different position in the model.The main phenomena are that aging evolution does not directly contribute to energetics when influencing only dissipative "friction" µ.This is similar to a cam-clay model [12,13] where the dissipative response is controlled through an internal variable whose rate, however, does not explicitly contribute to energetics.On the other hand, damage (or phase-field) influences the elastic response through the elastic response in the stored energy and is also driven by the resulting driving force from it.Also, we adopt the (realistic) assumption that the elastic strain (as well as its rate) is small, which makes possible to let µ dependent on the plastic strain rate rather than elastic strain rate and to put it into the standard framework of rate-dependent plasticity.The plasticity is consider without any hardening which otherwise might dominate with big slips on long time scales and would unacceptably corrupt the autonomy of the model.In principle, damage may also influence friction µ like in [46,47] but we will not consider it.

The one-dimensional steady-state model
It is generally understood that fracture mechanics and in particular fault mechanics is very complex and difficult to analyze.Therefore, we focus to a very simplified situation: a flat fault which is perfectly homogeneous in its tangential direction.Thus all variables depend only on the position in the normal direction and the problem reduces to be one dimensional, cf.We ask a question about existence of a steady state in the situations where the sides of the fault move with a constant speed in opposite directions.The model is thus expressed in rates rather than displacements and plastic strains.Such steady states are also called aseismic slips (sliding), in contrast to seismic slips which are dynamical phenomena related with a stick-slip motion and earthquakes.For the relation of the aseismic slip (fault growth) and orientation of faults see [40].The aseismic slip can be also understood as creep, within which the Maxwellian viscoelastic rheology is manifested.
The variables of our steady-state model will thus be: • π plastic strain rate (in 1/s), • ε elastic strain (dimensionless), • α damage (dimensionless, ranging over [0, 1]), and • θ aging (in seconds), and later also • σ a stress (or, in one-dimensional case, rather a force in J/m=N).These first five variables are to satisfy the following system of five equations (inclusions): where (•) x denotes the derivative (later also partial derivative) in x.Actually, (2.5c) contains a set-valued term ∂ π R(π, θ) = µ(π, θ)Sign(π) and is thus an inclusion rather than an equation.There, we have denoted by " Sign" in set-valued sign function, i.e.
−1 for π < 0. (2.6) This system arises as a steady state from an evolution model (2.12) below.In particular, the equation (2.5b) arises from the additive (Green-Naghdi's) decomposition of the total strain into the elastic strain and the plastic strain, cf.(2.12b) below.Written in terms of rates and taking into account that the rate of the elastic strain is zero in the steady state, we arrive at (2.5b).In fact, the velocity v here enters the rest of the system only through the boundary condition (2.8) below, in contrast to the full evolutionary model later in Section 4 where velocity acts through the inertial force.The data (or constitutive relations) in the model (2.5) are: µ = (π, α) a yield stress (in the one-dimensional model in N=J/m)), C = C(α) elastic modulus (smooth, nondecreasing, in N=J/m), f 0 aging rate (dimensionless), f 1 "contra-aging" coefficient (in seconds), G c fracture toughness (in a one-dimensional model in N=J/m), η > 0 a length scale coefficient for π (i.e. for the cataclastic zone, in W/m), > 0 a length scale coefficient for the damage (in meters), κ > 0 a length scale coefficient for the aging (in m 2 /s), while f 0 and f 1 are essentially borrowed from (2.3).Actually, v in (2.4) has the meaning rather of a difference of velocities across the contact interface than a velocity itself which would not be Galilean invariant.In a variant of the bulk model, µ should depend rather on a shear rate and, instead of the coefficient 1/v ref , one should consider a h/v ref with h a certain characteristic width of the active slip area, likely to be identified with the width of the cataclastic core zone, cf. Figure 2.1.Thus, we consider In comparison with (2.2), the steady-state equation (2.5e) contains the length-scale term κθ xx .Also damage equation (2.5e) contains a length-scale term 2 α xx competing with the driving force 1  2 C (α)ε 2 coming from the α-dependence in (2.5a).Note that the gradient term in (2.5c) applies to plastic rate and no gradient term involves directly the plastic strain, similarly as in [17,45].This eliminates spurious hardening-like effects by large slips accumulated on faults in large time scales, which would otherwise start dominating and corrupt the autonomous character of the model.
We have to complete the system (2.5) by suitable boundary condition.Specifically, we choose the boundary conditions with θ ∞ from (2.3).Let us mention that we use the mathematical convention that α = 1 means undamaged material while α = 0 means maximally damaged material.From (2.5a), we can see that C(α)ε is constant on the damage domain D = [−H, H], say = σ.From this, we can express for all x ∈ D . (2.9) If C(•) is increasing, one can conversely express α as a function of ε, but we will eliminate ε rather than α.Also the equation (2.5b) can be eliminated because the velocity v occurs only in the first boundary condition in (2.8).This condition then turns into an integral side constraint We can thus reduce (2.5) to the system of three elliptic ordinary-differential equations with the integral and the boundary conditions (2.11c) It is noteworthy that (2.10b) decouples from (2.10a,c) which arises not from necessity but rather from our desire for simplicity and for consistency with the standard rate-andstate friction as in Section 1: we assumed that µ, f 0 , and f 1 are independent of α.The system (2.10a,c)-(2.11a,c)thus represents a nonstandard non-local two-point boundaryvalue problem for the functions (π, θ) on D and one scalar variable σ.When solved, the two-point boundary-value problem (2.10b)-(2.11b)can be solved for α.Then ε is obtained from (2.9).Eventually, the velocity v can be calculated from (2.5b) when using also (2.11a).

The evolutionary model
We will now investigate an evolution version of the steady-state model (2.5), which in particular explains how (2.5) have arisen.In addition to the variables needed in Section 2.2, we now will exploit also: • p plastic strain (dimensionless) and • mass density (in one-dimensional model kg/m).An additional ingredient will be a dissipation potential ζ for damage, which is convex with subdifferential ∂ζ and has physical dimension J/m.The evolution variant of (2.5) then looks as: It is to be completed with boundary conditions as (2.8) with possibly time dependent boundary velocity with θ ∞ constant in time.The (Green-Naghdi's) additive decomposition is written in rates, which just gives (2.12b).Obviously, the steady-state variant of (2.12) where all time derivatives vanish yield just (2.5).
The system (2.12a-d) has a rational physical background while (2.12e) expresses some extra phenomenology controlling the nonconservative part in (2.12c).For = 0, the system (2.12a-d) represents the so-called Biot equation ∂ .
The underlying specific stored energy and the dissipation potential (in terms of the rates of plastic strain p and damage α) behind this model are x and (2.14a) where often C(α) = ( 2 / 2 0 + α 2 )C 0 with some 0 .The constants and 0 are in meters while the fracture toughness G c is in J/m 2 , cf. [30,Eqn. (7.5.35)], or rather in J/m in our 1-dimensional model.This is known as the Ambrosio-Tortorelli functional [2].

Analysis of the steady state model
Further on, we will use the standard notation for the function space.In particular, C(D) will be the space of continuous functions on D and L p (D) will denote the Lebesgue space of measurable functions on the domain D = [−H, H] whose p-power is integrable (or, when p = ∞, which are bounded), and W k,p (D) the Sobolev space of functions in L p (D) whose k-th distributional derivative belongs to L p (D).We abbreviate Besides, H 1 0 (D) will denote a subspace of H 1 (D) of functions with zero values at x = ±H.In Section 4, for the time interval I = [0, T ] and a Banach space X, we will also use the Bochner spaces L p (I; X) of Bochner-measurable functions I → X whose norm in in L p (I), and the Sobolev-Bochner space H 1 (I; X) which belong, together with their distributional time derivative, into L p (I; X).

Existence of steady states
Let us recall the standard definition of a weak solution to the inclusion (2.5c) as a variational inequality to be satisfied for any π ∈ L 1 (D), where ∂ π R(π, θ) = µ(π, θ)Sign(π).We will prove existence of solutions due to even a stronger concept of a classical (also called Carathéodory or strong) solution, namely that |π xx | is integrable (actually in our case even bounded) and holds a.e. on D. As mentioned in Section 1, the rate-and-state friction model lacks standard thermodynamical consistency, which is reflected in the steady-state case by a lack of joint variational structure.Nevertheless, the two equations (2.5c) and (2.5e) for π and θ, respectively, have an individual variational structure governed by the functionals where ϕ 0 and ϕ 1 are primitive functions to f 0 and f 1 , respectively.Then, the pair (θ, π) is a desired solution if and only if Since both functionals A π (•) and B θ (•) are strictly convex, the solutions operators θ = S A (π) = argminA π and π = S B (θ) = argminB θ (•) are well-defined.The existence of steady states will be proved by a Schauder fixed-point theorem applied to S A • S B .Theorem 3.1 (Existence of steady states).Let the following assumptions hold: and non-increasing on (−∞, 0], inf R µ(0, θ) > 0, (3.4a) (ii) Moreover, any solution satisfied Proof.For a given θ, equation (2.10a) with the nonlocal condition in (2.11) is equivalent to π = S B ( θ) = argminB θ (•).The monotonicity of µ(•, θ) assumed in (3.4a) ensures the uniform convexity of the functional B θ (•).Therefore the minimizer π = S B ( θ), which clearly exists by the direct method in the calculus of variations, is uniquely determined.Moreover, it depends depends continuously on θ with respect to the weak topology on H 1 (D).Thanks to (3.4a), for v ∞ given, B θ (•) is coercive uniformly with respect to θ, and therefore the minimizer π = S B (θ) can be a priori bounded in H 1 (D) independently on θ.
For a given π, equation (2.5e) is equivalent to θ = S A (π) = argminA π (•).As f 1 is nondecreasing and f 0 is nonincreasing, the functional A π (•) is convex, and it is to be minimized on the affine manifold {θ ∈ H 1 (D); θ(±H) = θ ∞ }, cf. the boundary conditions (2.11).Therefore this boundary-value problem has a unique weak solution θ ∈ H 1 (D), which depends continuously on π and can be bounded independently of θ when taking into account the mentioned a priori bound for π.
Altogether, we obtain a mapping θ → θ = S A S B (θ) which is continuous with respect to the weak topology on H 1 (D) and valued in some bounded set (depending possibly on a given v ∞ ).By the Schauder fixed-point theorem, this mapping has a fixed point θ.This thus determines also π = S B (θ) and σ.
We discuss further qualitative properties of solution pairs (θ, π) that arise from the specific form of the steady state equations (2.5)-(2.8).As our above result does not imply uniqueness of solutions, our next results states that there are solutions with symmetry and, under a weak additional condition, these solutions are also monotone on [0, H].For the latter we use the technique of rearrangements, which strongly relies on the fact that we have no explicit x-dependence in our material laws.For general function f ∈ L 1 (D) we define its even decreasing and even increasing rearrangements f dr and f ir via then there exists an even, monotone pair (θ, π), i.e. it is an even pair such that additionally [0, H] x → θ(x) is nondecreasing and [0, H] x → π(x) is nonincreasing.
To obtain the evenness we simply restrict the existence theory developed in the proof of Theorem 3.1 to the closed subspaces of even functions.By the uniqueness of the minimizers of A π and B θ it is clear that S A and S B map even functions to even functions.Hence, Schauder's fixed-point theorem produces an even solution.
For showing the existence of monotone pairs we rely on classical results for rearrangements, see e.g.[29], namely the Polya-Szegö inequality To exploit the theory of rearrangements we define the closed convex sets For the last identity we use ϕ 1 (θ) dr = ϕ 1 (θ dr ) which holds because of Summing the three relations gives A π (θ dr ) ≤ A π (θ).
Similarly, we derive S B : For this we use assumption (3.6), which gives R(π, θ) = R(π, 0) + B(θ)|π|, and the three relations  1, real faults may go into so-called aseismic slip (also called aseismic creep), where one observes pure sliding like predicted by our steady state solutions constructed above.However, for our simplified evolutionary model introduced in Section 5 (cf.(5.1)) numerical simulations predict instability of the steady state and the development of stick-slip oscillations, see Section 5.4.In the former case, stresses remain low and never challenge the plastic yield stress µ(0, θ ∞ ) at the core of the faults, a fact which is unfortunately not covered by our model.One possible modification for modeling this effect would be to replace the setvalued Sign(•) in (2.5c) by some monotone smooth approximation, e.g.π → tanh(π/δ) with 0 < δ 1.
3.2 Asymptotics of the plastic zone for η → 0 and κ → 0 The gradient term in (2.5c) and in (2.10a) controls in a certain way the width of the cataclastic zone where the slip is concentrated.There is an expectation that, when suppressing it by η → 0, the slip zone will get narrower.It is however a rather contra-intuitive effect that the zone eventually does not degenerate to a completely flat interface like it would be in so-called perfect plasticity where the plastic strain rate π would be a measure on D. Here, in the limit, π only looses its W 2,∞ -regularity as stated in Theorem 3.1 for η > 0 but remains in L 1 (D).
In this section, let us denote the solution obtained as a Schauder fixed point in the proof of Theorem 3.1 by (ε η , v η , π η , α η , θ η , σ η ).There is a subsequence such that, for some Proof.From the proof of Theorem 3.1, we can see that the a priori bounds for are independent of η > 0 and π η H 1 (D) = O(1/ √ η).Moreover, from π η = S B (θ η ), we can easily see that even R(π η , θ η ) is bounded in L 1 (D).Using (3.10a) we can apply the criterion of de la Valleé Poussin [18] and obtain that {π η } η>0 is weakly compact in L 1 (D).
Then the limit passage in the weak solution to (2.5)-(2.8)for η → 0 is quite easy.The only nontrivial point is the limit passage in the variational inequality (3.1).We first use The liminf estimate follows because R(•, θ) is convex and continuous such that D R(•, θ) dx is weakly lower semicontinuous on L 1 (D).The penultimate integral in (3.12) converges to 0 because θ η → θ uniformly on D due to the compact embedding The variational inequality (3.12) does not contain any x-derivatives any more and hence is equivalent to the pointwise inequality R( π, θ(x)) − σ( π−π(x)) ≥ R(π(x), θ(x)) a.e. in D. But this is equivalent to σ ∈ ∂ π R( π(x), θ(x)) and hence (3.11g) holds.
We are now ready to study the limit κ → 0 as well, which is really surprising because we are losing all control over spatial derivatives and all the modeling length scales induced by η and κ tend to 0. In such a situation the usual compactness arguments fail and fast spatial oscillations, i.e. microstructures, may appear.Indeed we will see in Remark 3.6 that there are indeed many complicated solutions without any length scale.However, it is surprising that it is possible to show that natural solutions exist, namely even, monotone pairs (θ, π).The idea is to use for κ > 0 and η = 0 the even, monotone pairs (θ κ , π κ ) obtained from Proposition 3.2 and the subsequent limit η → 0 in Proposition 3.4.The monotonicity of the pairs (θ κ , π κ ) allows us to deduce pointwise convergence, which is good enough to pass to the limit κ → 0 even in nonlinear functions.
Here the monotonicities are kept, i.e. θ 0 = θ ir and π 0 = π 0 dr , but the continuity of the limits might be lost.Moreover, π 0 (0) = ∞ might be possible.Step 3. Limit passage in the equations: Since Π is continuous, the pointwise convergence yields the limit relation π 0 (x) = Π(σ 0 , θ 0 (x)) for all x ∈ D. (3.17) For the equation determining θ we can use the a priori estimate κ θ κ 2 L 2 ≤ C * and pass to the limit in the weak form of (κ θ κ x ) x + f 0 (θ κ ) = π κ f 1 (θ κ ), i.e. in the integral identity This provides the pointwise relation From (3.17) and (3.18) we immediately see that (3.13) holds.We next observe that θ = Θ f (π) is well-defined by the implicit function theorem using (3.4c).Thus, the solutions satisfy θ 0 (x) = Θ f (π 0 (x)) for a.a.x ∈ D. Henceforth, recalling µ(π) = µ(π, Θ f (π)), the minimization problem (3.13) is equivalent to σ ∈ µ(π)Sign(π) and D π dx = 2v ∞ .Defining the function R(π) = π 0 µ(s) ds, this is equivalent to the following problem: However, this minimization problem is well understood via the convex hull R * * , see [9,Ch. 2].By our assumption (3.14) we know that R * * has the form and satisfies R * * (π) R(π) for π ∈ (0, π * ) and R (π) > 0 for π ≥ π * , see Figure 3.2.As our R is superlinear, a minimizer always exists.Moreover, recalling that v ∞ /H > 0 is the average value of π : D → R, the minimizer is unique if and only if the tangent at π = v ∞ /H is not in the interior of an interval on which R * * is affine.In the open interval (0, v ∞ /H) the minimizers π attain only the values 0 and π * on sets with the corresponding measures to fit the average.However, by constructing the even, nonincreasing rearrangement, we find a unique minimizer, where only the value at the two jump points x = ±h = v ∞ /π are free.
From these uniqueness results we also obtain the convergence of the full family by the standard contradiction via compactness.With this, Proposition 3.5 is established.
Remark 3.6 (Nonuniqueness of solutions).We want to emphasize that the uniqueness result for κ = η = 0 at the end of Proposition 3.5 concerns only even, monotone solutions.
Because of κ = η = 0 there are indeed infinitely many solutions, as we can "rearrange" the function values of (θ, π) freely.In the case v ∞ < π * H, we can choose any open set P ⊂ D with D 1 P dx = 2v ∞ /π * and the function is a solution of (3.13) as well.

Analysis of the evolutionary model
We now consider the evolutionary model (2.12).The energetics (2.14) behind this model can be revealed by testing momentum balance (2.12a) by x/H, the plastic flow rule (2.12b) by .p, and the damage rule (2.12c) by .
α. Using the Dirichlet boundary condition for the velocity at x = ±H, we have v(±H) = 0, as needed.The first test gives, in particular, the term where also (2.12b) has been used.This test of the inertial form gives Combining it with the tests of (2.12b) by .p and of (2.12c) by .
α which give we altogether obtain the energy balance x kinetic and stored energies dx where τ ∈ R 2 is the traction on the boundary (i.e.here two forces at x = ±H) defined as a functional τ, (z(H), z(−H)) = D .vz + C(α)εz x dx for any z ∈ H 1 (D), cf.e.g.[30,Sect.6.2].
Further on, we will be interested in an initial-value problem.For this, we prescribe some initial conditions, i.e.
A definition of the weak solutions of particular equations/inclusions in (2.12) can be cast by standard way, using convexity of the involved functionals.Let us specify, rather for illustration, the weak formulation for the inclusion (2.12c) exploiting that µ( .p, θ)Sign( .
p.This leads to the variational inequality to be valid for any π ∈ L ∞ (I×Ω).
Beside the previous assumptions, we now also assume The definition of weak solutions to (2.12) with (2.13) and (4.4) is standard and we will not write it explicitly; the variational inequality (3.1) is to hold integrated over I. Furthermore, we also exploit the superlinear growth of R(•, θ) from (3.10a), namely which is a standard estimate for µ ∈ ∂ψ(π), Note that the standard model (2.4) complies with assumption (3.10a).
Relying formally on the tests leading to (4.3), after integration in time on the interval [0, t] when using also the by-part integration, we obtain Moreover, the aging equation (2.12e) has to be tested separately by using the test function θ−θ ∞ , which has zero traces for x = ±H.Integrating the result over [0, t] leads to When summing (4.8) and (4.9), we can use the Hölder and a (generalized) Young inequality to estimate the resulting right-hand side.Actually, the only nontrivial term is 4.9) and it can be estimated as where Φ * is the Fenchel-Legendre conjugate of Φ, i.e.Φ * (s) = sup π∈R πs − Φ(π) .The term 1 2 µ( .p, θ)| .
p| in (4.10) can then be absorbed in the left-hand side of (4.8) while Eventually, the last term in (4.8) can be estimated as (1+|v| I) and using Gronwall's inequality, from the left-hand sides of (4.8) and (4.9) we can read the a priori estimates ) By comparison, we will get also an information about .
, and also about The rigorous existence proof of weak solutions is however very nontrivial and seems even impossible for the full dynamical model (2.12) with damage.Some modifications by involving some additional dissipative terms or some higher-order conservative terms seem necessary, cf.[30,Sect.7.5] or also [46] for the model without aging.Consistently also with the computational experiments in Section 5 below, we thus present the rigorous proof only for a model without damage, i.e. for C > 0 constant.Theorem 4.1 (Damage-free case -existence and regularity of solutions).Let (3.4a,c,d) with µ smooth, (4.6), and (4.7) hold, and > 0 be a constant and v ∞ ∈ W 1,1 (I).Then: (i) There is a weak solution to the initial-boundary-value problem for the system (2.12a-c,e) with the boundary conditions (2.13) and the initial conditions (4.4).
(ii) If sup 0≤θ≤θ∞ µ(•, θ) does not have a growth more than O(|π| s ), then these solutions are, in fact, regular in the sense that p ∈ W 1,s (I; and also each such weak solution satisfies the energy balance (4.3) without α-terms integrated over a time interval [0, t] with any t ∈ I.
Let us note that the O(|π| s )-growth condition in the point (ii) surely covers the model (2.7) for any 1 ≤ s < ∞.
Sketch of the proof.Actually, the above formal procedure is to be made first for a suitable approximation whose solutions exist by some specific arguments, and then to pass to the limit.Imitating the split for the static problem used in the proof of Theorem 3.1, we choose a staggered time discretization.We take an equidistant partition of the time interval I by using the time step τ > 0, assuming T /τ integer and considering a sequence of such τ 's converging to 0. Then, recalling ∂ π R(π, θ) = µ(π, θ)Sign(π), we consider a recursive boundary-value problem for the system to be solved for k = 1, 2, ..., T /τ starting for k = 1 from the initial conditions v 0 τ = v 0 , ε 0 τ = ε 0 , and θ 0 τ = θ 0 .The boundary conditions for (4.12) are like in (2.8) but now with time-varying velocity v ∞ , i.e.
The system (4.12a-c) has a variational structure with a convex coercive potential For a sufficiently small τ > 0, this potential is convex and coercive on L 2 (D) 2 × H 1 (D).Minimization of this functional on an affine manifold respecting the boundary conditions v(±H) = ±v k ∞ , π(±H) = 0, and θ(±H) = θ ∞ gives by the standard direct-method argument existence of an (even unique) minimizer, let us denote it by . This minimizer satisfies (4.12a,b) in the weak sense and also the inclusion xx in the weak sense.Then we can solve (4.12d) by minimization of the convex functional where ϕ i are the primitive functions to f i , i = 0, 1.This functional is coercive on a linear manifold of the space H 1 (D) respecting the boundary condition (2.8).Let us denote its unique minimizer by θ k τ .We introduce the piecewise affine continuous and the piecewise constant interpolants.Having {v k τ } T /τ k=0 , we define for (k−1)τ < t ≤ kτ with k = 0, 1, ..., T /τ .Analogously, we define also ε τ , or θ τ , etc.This allows us to write the system (4.12) in a "compact" form: . . .
By modifying appropriately the procedure which led to the a priori estimates (4.11ac,e), we obtain here ≤ C , and here also (4.18d) All these estimates hold also for the piecewise affine interpolants, and (4.18d) holds also for θ τ .The last estimate is obtained by comparison from ξ τ = (Cε τ +η(π τ ) xx )/µ(π τ , θ τ ) when testing it by functions bounded in L 2 (I; H 1 (D)) and using the smoothness of 1/µ(π τ , θ τ ).Then, by the Banach selection principle, we obtain subsequences indexed, for simplicity, again by τ ) weakly* converging in the topologies indicated in (4.18), and we pass to a limit for τ → 0 and are to show that such limit (let us denote it by (v, ε, π, θ, ξ)) solve the continuous problem with π = .
p.For this, one uses the Aubin-Lions compactness theorem adapted for the time-discretization method as in [43,Sect. 8.2].Thus we can rely on that The limit passage in the linear hyperbolic equation (2.12a) is due to a weak convergence of both v and ε and also the limit passage in the linear equation (2.12b) is easy via weak convergence.Yet, there is one peculiarity in the limit passage in the nonlinearity in (2.12c) for which a strong convergence of ε is needed, but we do not have any information about space gradient of ε.The other peculiarity is a need of the strong convergence of .
p which is needed for (2.12e), but we do not have any information about .
For the equality in (4.20), we used (4.17b) together with its limit obtained by the weak convergence, i.e. .ε = v x − π, and also (4.17a,c) for the identity It is important, that (4.20) holds for any ξ ∈ Sign(π) and, at this moment, we do not assume that ξ comes as a limit from the (sub)sequence {ξ τ } τ >0 .
p k and the fixed boundary conditions, so that we do not need to rely on the monotonicity of ∂ π R(•, θ) which may not be strong.
Having the strong convergence (4.21) at disposal, the limit passage is then easy, showing that the previously obtained weak limit (v, ε, π, θ) is a weak solution to the system (2.12).In particular, from the inclusion in (4.17c) one obtains ξ ∈ Sign(π) by using maximal monotonicity of the graph of the set-valued mapping Sign : L 2 (I; H 1 (D)) ⇒ L 2 (I; H 1 (D) * ) and the strong convergence (4.21c).Thus (i) is proved.
If s ≥ 2, the procedure which led to the energy balance (4.3) considered here without α-terms but integrated over a time interval [0, t] was indeed rigorous.This is because v ∈ L 2 (I; H 1 (D)), as can be seen by comparison from (2.12b), is in duality with .
v ∈ L 2 (I; H 1 (D) * ) and with (Cε) x ∈ L 2 (I; H 1 (D) * ), so that testing the momentum equation (2.12a) and the related by-part integration is legitimate.Similar arguments concern also the aging rule (2.12e).Since ηπ xx ∈ L 2 (I×D) if s ≥ 2, also the test of the plastic rate equation (2.12c) by π ∈ L s (I×D) is legitimate together with the related by-part integrations.
In this case when s ≥ 2, also (4.17d) can be tested by .
Remark 4.2 (Stability and time-periodic solutions).In geodynamics the phenomenon called episodic tremor and slip describes time-periodic motions in subduction zones where shorter periods of plastic slips alternate with longer periods with slow slip events.Hence, it would be interesting to complement our existence result for "transient events" governed by the above initial-value problem by a theory for time-periodic solutions.The aim would be show that there is a period t * > 0 and a solution of the system (2.12) with the boundary conditions (2.8) satisfying ( .v, .ε, .
θ) ≡ 0 and instead of (4.4).Of course, a general question is that of stability of the steady state solutions (π, θ) obtained in Section 3 or potentially of such time-periodic solutions as described here.As we will see in the following section, one indication of the existence of time-periodic solutions is the loss of stability of the steady state solution.But because of the complexity of the model, these questions are beyond the scope of this paper.
Remark 4.3 (Asymptotics for η → 0 and κ → 0).Unlike to the case for steady solutions for (2.5) as in Section 3.2, it is not possible in the evolutionary model (2.12) to pass to the limit for η → 0. In particular, a limit passage in the term Cε η ( π η − . p η ) occurring in (4.5) seems to be out of reach.The substitution (4.2a) by a convex term in .
p could not help, being not weakly upper-semicontinuous.If also (3.10a) holds, then like in Propositions 3.4 and 3.5, we can at least obtain some uniform bounds, in particular for the plastic strain rate π = .p in the Orlicz space L Φ (I×D) with Φ from (3.10a), i.e.I D Φ(π(t, x)) dx dt < ∞.Yet, the limit passage for η → 0, even while keeping κ > 0 fixed, remains intractable.

Illustrative numerical simulations
We illustrate the response of the evolutionary model in Section 4 by a simplified model derived in Section 5.1.This model still has exactly the same steady states as the full model, such that all the theory of Section 3 applies to it, when ignoring statements about the damage variable α.We expect that the simplified model is still relevant as far as usually observed dynamical features concern.Moreover, it also displays the effect of the free boundary occurring between the elastic zone and the plastic zone.In Section 5.2 we show by numerical simulations that the steady states localize for v ∞ → 0 in such a way that π stst has support (i.e. the so-called cataclastic zone) in . Moreover, we show that, when keeping v ∞ = 0 fixed but sufficiently small, we obtain a support with h * (v ∞ , κ) → v ∞ /π * for κ → 0 + .In Section 5.3 we study an ODE model for scalars θ(t) and σ(t) which displays the effect of oscillatory behavior for |v ∞ | < v crit while solutions converge to the unique steady state for |v ∞ | > v crit .Finally Section 5.4 presents simulations for the simplified evolutionary model.In particular, we observe again that for small nontrivial values of |v ∞ | we have oscillatory behavior, where the plastic zone is spatially and temporarily localized in the sense that the support of π(t, •) is compactly contained in D = [−H, H] for all t ∈ [0, T per ] and that π(t, x) = 0 for all x ∈ D and all t ∈ [t 1 , t 2 ] for a nontrivial interval [t 1 , t 2 ] ⊂ [0, T per ].For |v ∞ | large, we find convergence into a steady state with a nontrivial plastic (cataclastic) zone.All the following results are derived from numerical experiments only.

The simplified model without damage
To display the main features of our rate-and-state friction model we reduce the full evolutionary model (2.12) by making the following simplifications: • we neglect inertial effects (i.e.we set = 0 in (2.12a)), thus making the system quasistatic but still keeping a rate-and-state dependent plasticity; • we choose η = 0 for the length-scale parameter in (2.12c) as analyzed in Section 3.2 for the steady-state solutions; • we neglect all damage effects through α and omit (2.12d) as we did in Theorem 4.1.Because of = 0, the momentum balance leads to a spatially constant stress σ(t) = Cε.As now C is constant, also ε(t) is spatially constant.Integrating (2.12b) over x ∈ D = [−H, H] and using the boundary condition for v from (2.13) gives the following coupled system for σ, π = .
Here the nonsmoothness due to the plastic behavior is realized by the nonsmooth function π = Π(σ, θ) defined in (5.2).For all the following simulation we choose the following parameters and functions: Subsequently, we will only vary the coefficient κ > 0 and the shear velocity v ∞ .

Steady states
We first discuss the steady states for (5.3), which are indeed a special case of the steady states obtained in Proposition 3.4.Numerically, we always found exactly one steady state θ stst = Θ(v ∞ , κ), but were unable to prove its uniqueness rigorously.When varying the parameters v ∞ and κ we can easily observe clear trends for (θ stst , π stst ), where the associated plastic flow rate is given by π stst = P (σ stst , θ stst ), see Figure 5.1.We first observe that for fixed κ the functions θ stst and π stst depend monotonically on v ∞ in the expected way, namely θ stst decreases with the shear velocity v ∞ , while π stst increases, which fits to the relation 2v ∞ = D π stst (v ∞ , κ; x) dx.
Stationary profiles θstst of the aging variable Moreover, for v ∞ → 0 + the scaled plastic rate π stst /v ∞ converges to a nontrivial limit with localized support, while θ stst converges uniformly to θ ∞ .For larger and larger v ∞ the plastic zone occupies more and more of the domain D = [−1, 1] and θ stst is very small in most of the plastic zone, namely θ When reducing the size of κ we also see that the size of the plastic zone shrinks.For small v ∞ it can be seen that the support of Finally, we want to study the case corresponding to Proposition 3.5, where v ∞ is kept fixed and the limit κ → 0 is performed.In Figure 5.3 we show plots of the steady states (θ κ stst , π κ stst ) for three different values of v ∞ for a sequence of decreasing κ.We clearly see the predicted development of convergence against towards the limit (θ 0 stst , π 0 stst ) taking only two different values.Moreover, the values are roughly independent of v ∞ , where the active plastic zone (−h, h) behaves like h = v ∞ /π * , as proved in Proposition 3.5.

An ODE model showing oscillations in time
Oscillatory behavior is most easily seen in a simple finite dimensional model, consisting only of σ(t) and θ(t), where we may consider θ(t) as the average of θ(t, x) over the critical plasticity region where π(t) = P (σ(t), θ(t)) is positive.We also refer to the analysis of a spring-slider model in [33] as well as the geophysical paper [1].
Rescaled stationary profiles πstst/v∞ of the plastic strain rate (5.5) Here h ∈ ]0, H[ represents the width of the plastic zone, which has to be adapted accordingly.We may consider (5.5) as an evolutionary lumped-parameter system, which in geophysical literature is often referred to as a 1-degree-of-freedom slider and is considered as a basic test of every new friction model.The nice feature of this ODE model is that the steady states can be calculated explicitly, and even a stability analysis can be performed.Indeed there is exactly one steady state, namely Instead of performing a rigorous analysis, we simply display the solution behavior of this ODE by a few numerical results.We find that for small positive v ∞ we obtain oscillatory behavior, while for larger v ∞ the solutions converge to the steady state, see Figure 5.4.Indeed, the oscillations can be interpreted physically in terms of geophysical processes as seismic cycles.
During the oscillatory behavior there is a large part of the interval where there is no plastic slip (i.e.π(t) = 0).In these intervals the stress is growing linearly with a slope that is proportional to v ∞ , and the aging variable θ is relaxing exponentially back to its equilibrium value θ ∞ .However, if the stress reaches a critical value, then the plastic strain rate is triggered, which leads to reduction of the aging variable.This leads to a simultaneous weakening of the plastic yields stress µ(π, θ) such that π can grow even more.As a result the stress is drastically reduced in a rather short time interval, and θ is reduced almost down to 0 (refreshing).If the inertial term would be included, then this fast rupture-like processes could emit elastic waves, i.e. earthquakes.Because of the stress release the plastic strain rate reduces to 0, and the process starts again by a slow aging and building up the stress.In fact, choosing h = 0.3 a closer analysis of the system shows that the steady states are stable if and only if v > v (1) ∞ ≈ 0.17462.However, stable oscillations are already seen for v < v (2) ∞ ≈ 0.175452.A careful analysis of the trajectories in the phase plane for (θ, σ) reveals that for v ∞ ∈ (v ∞ ) there are two periodic solutions, as smaller unstable one that encircles the stable fixed point and a larger stable one that encircles the unstable one, see Figure 5.5.Thus, in the small parameter interval (v ∞ ) we have coexistence of a stable fixed point and a stable periodic orbit.

Convergence to steady states versus oscillations for (5.3)
The behavior of the evolutionary coupled system (5.3)coupling the parabolic PDE for the aging variable θ(t, x) to the ODE for the stress σ(t) displays roughly a similar behavior as the lumped ODE system (5.5).For large |v ∞ | one observes convergence into the steady states analyzed in Section 3 and displayed numerically in Section 5.2.For small nontrivial values of v ∞ one observes oscillatory behavior.Of course, the new feature is the spatial distribution of the plastic rate π(t, x) = Π(σ(t), θ(t, x)) and the aging variable θ(t, x).In most cases one observes that π(t, x) has a nontrivial support in the sense that the support of π(t, •) is compactly contained in (−H, H).Moreover, in the oscillatory case, we also observe that there are large parts of the periodicity interval, in which there is no plastic flow at all (i.e.π = .p = 0), but there is aging and slow building up of stress.The outer one is stable and is approached by the blue trajectories from inside and outside.The unstable periodic orbit lies between the orange and the brown trajectory.
Then, in sudden plastic bursts there is a strong plastic flow that leads to stress release and refreshing, i.e. reduction of θ almost down to 0 inside the cataclastic zone.In the case κ = 0.04 and the smaller shear rate v ∞ = 0.15 one observes oscillatory behavior.In fact, we start the solution very close to the steady state and the solution needs some time to develop the instability but then it switches quickly into a periodically looking regime, see Figure 5

Figure 3 . 1 :
Figure 3.1: Two examples of functions f and their decreasing and increasing rearrangements f dr and f ir .

7 )
and the Hardy-Littlewood inequality (cf.[24,Ch. 10]) D f dr g ir dx = D f ir g dr dx ≤ D f g dx ≤ D f dr g dr dx = D f ir g ir dx.(3.8)While the upper estimate is classical and works for integration overD = B R (0) ⊂ R d or D = R d ,the lower estimate is special to D ⊂ R 1 , see [24, Eqn.(10.2.1)].

π
dr B(θ) dr dx = D |π| B θ dr dx, where we used that B is nondecreasing.This finishes the proof of existence of even, monotone pairs.Remark 3.3 (Aseismic-slip regime).Under very low shear velocities |v ∞ |

18 Figure 5 . 4 :
Figure 5.4: Solutions (θ(t), σ(t)) together with π(t) = P (σ(t), θ(t)) for h = 0.3 and three different values of v ∞ .In the first two cases the solutions start very close to the unstable steady state.In the third case the solution starts far away but soon returns to the stable fixed point.

Figure 5 . 5 :
Figure 5.5: The (σ, θ) phase plane for h = 0.3 and v ∞ = 0.175, where all trajectories rotate clockwise around the fixed point (σ stst , θ stst ) ≈ (1.973, 0.168).There are two periodic solutions.The outer one is stable and is approached by the blue trajectories from inside and outside.The unstable periodic orbit lies between the orange and the brown trajectory.