# Numerical simulation of multiscale fault systems with rate- and state-dependent friction

We consider the deformation of a geological structure with non-intersecting faults that can be represented by a layered system of viscoelastic bodies satisfying rate- and state-depending friction conditions along the common interfaces. We derive a mathematical model that contains classical Dieterich- and Ruina-type friction as special cases and accounts for possibly large tangential displacements. Semi-discretization in time by a Newmark scheme leads to a coupled system of non-smooth, convex minimization problems for rate and state to be solved in each time step. Additional spatial discretization by a mortar method and piecewise constant finite elements allows for the decoupling of rate and state by a fixed point iteration and efficient algebraic solution of the rate problem by truncated non-smooth Newton methods. Numerical experiments with a spring slider and a layered multiscale system illustrate the behavior of our model as well as the efficiency and reliability of the numerical solver.

## Authors

• 7 publications
• 3 publications
• 2 publications
10/29/2021

### An Assessment of Solvers for Algebraically Stabilized Discretizations of Convection-Diffusion-Reaction Equations

We consider flux-corrected finite element discretizations of 3D convecti...
04/04/2020

### Numerical Analysis of History-dependent Variational-hemivariational Inequalities

In this paper, numerical analysis is carried out for a class of history-...
04/21/2020

### Gradient discretization of two-phase flows coupled with mechanical deformation in fractured porous media

We consider a two-phase Darcy flow in a fractured porous medium consisti...
04/06/2020

### A BDF2-Semismooth Newton Algorithm for the Numerical Solution of the Bingham Flow with Temperature Dependent Parameters

This paper is devoted to the numerical solution of the non-isothermal in...
02/16/2021

### A positivity-preserving and convergent numerical scheme for the binary fluid-surfactant system

In this paper, we develop a first order (in time) numerical scheme for t...
04/12/2021

### Analysis of algebraic flux correction for semi-discrete advection problems

We present stability and error analysis for algebraic flux correction sc...
06/06/2019

### Mathematical modeling of a Cosserat method in finite-strain holonomic plasticity

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Stress accumulation and release in geological fault networks play a crucial role in earthquake dynamics. The phenomenology of faults is ranging from subduction zones like the Nasca plate and strike-slip faults like the San Andreas fault to multiscale fault systems like the Atacama zone. Strongly varying time scales between the occurrence and duration of slip events suggest to complement experimental studies in the field (or in the lab [35]) by numerical simulations.

In the underlying mathematical description, the Dieterich–Ruina model of rate- and state-dependent friction (RSF) [37] has become a standard for the frictional behaviour along the faults [7, 32, 34]. It can be regarded as an extension of simple Tresca friction with rate- and state dependent friction coefficient that is increasing/decreasing with increasing/decreasing sliding velocity or slip rate involving some relaxation effect as expressed by the state . The variational structure of RSF has been identified and first exploited by Pipping et al. [30]

The simulation of rupture and slip events in seismic hazard analysis has quite a history (cf., e.g., [3, 8, 20, 24] and the references cited therein). Utilizing a discontinuous Galerkin (DG) scheme in space in connection with arbitrary high-order (ADER) time integration, de la Puente et al [10] developed a numerical method for the dynamic simulation of slip events. This method was later generalized to three space dimensions [26] and cast into the software package SeisSol that was successfully utilized for the simulation, e.g., of the 2016 Kaikōura earthquake cascade [42]. More recently, a different approach based on a diffuse representation of faults was first applied to subduction zones [18, 41, 43] and later extended to strike-slip faults [9]. This approach has the potential to allow for much more complicated fault systems because the faults have to be no longer resolved exactly by the underlying finite element mesh. However, this advantage currently comes with high computational cost due to a lack of efficient algebraic solution techniques.

In this paper, we extend a variational approach to the simulation of subduction zones [29] to a layered fault system with RSF. More precisely, we consider the deformation of a geological structure with non-intersecting faults that can be represented by viscoelastic bodies undergoing small viscoelastic deformations and large tangential displacements with RSF contact conditions. Assuming existence of a sufficiently regular contact mapping, we formulate a general mathematical model that contains Dieterich–Ruina friction as a special case. Fault opening is forbidden for notational convenience, but could be included in a straightforward way.

Time discretization is performed by a classical Newmark scheme, resulting in a coupled system of non-smooth, convex minimization problems that has to be solved in each time step. Decoupling this system by a fixed point iteration leads to a problem for the velocity (and thus for the rate) with given state, and an independent state problem with given rate. Both the rate and the state problem can be rewritten as convex minimization problems admitting unique solutions. For a related coupled problem with unilateral contact, as arising from the mathematical description of subduction zones, existence and uniqueness of solutions was established by Pipping et al. [27, 30] using fixed point arguments.

Spatial discretization of the rate problem is performed by a mortar method in the spirit of Krause and Wohlmuth [23, 46, 45]. This approach has the advantage that it provides nodal block separation of the non-smooth nonlinearity which allows for direct application of globally convergent Truncated Non-smooth Newton Multi-Grid (TNNMG) methods [14, 15, 16]. The state problem is discretized by piece-wise constant finite elements. For given rate, the resulting algebraic problem decouples into independent scalar problems for each of the nodal values, which can be solved, e.g., by simple bisection or even explicitly.

In our numerical experiments, we consider a spring slider with bodies and a layered network with bodies separated by 4 faults subject to prescribed velocities at the upper boundary. We perform self-adaptive time stepping to efficiently resolve strongly varying velocities during loading, rupture, and sliding. Spatial discretization is based on triangulations as obtained by adaptive refinement concentrated at the faults. The associated hierarchy of finite element spaces is used for the algebraic TNNMG solver of the rate problems with given state as arising in the fixed point iteration mentioned above.

For the spring slider we observe the periodic occurrence of mostly unilateral slip events, similar to related simulations of subduction zones [29]. These slip events are nicely captured by adaptive time stepping, while the number of outer fixed point iterations and inner multigrid iteration remains almost the same for all time steps. Simulation of the layered network exhibits an interesting coincidence of periodic slip events along the upper fault with loading phases and oscillatory behavior on the others. We observe essentially the same efficiency of time stepping, fixed point iteration, and multigrid as for the spring slider which illustrates the robustness of our numerical solution procedure, also with respect to the number of faults.

## 2. Mathematical modelling

### 2.1. A layered fault system with rate-and state-dependent friction

We consider a geological structure containing a system of faults which is represented by a deformable body with reference domain , , that, along the faults, is decomposed into subdomains , ,

 ¯¯¯¯Ω=I⋃i=1¯¯¯¯Ωi.

We assume that these subdomains are non-empty, bounded Lipschitz domains, do not penetrate each other and are layered in the sense that at most two subdomains are in contact at any point in (see Figure 1). Then, the subdomains can be ordered such that there is a common interface , , and all other intersections of subdomains are empty. Setting for notational convenience, the boundary of is disjointly decomposed according to , into a Dirichlet, a Neumann, and a contact boundary, respectively. We set

 ΓD=I⋃i=1ΓDi,ΓN=I⋃i=1ΓNi,ΓF=I−1⋃i=1ΓFi,i+1.

For with , , we define the restrictions and of to with

 vT,i=vi+1|ΓFi,i+1,vB,i=vi|ΓFi,i+1i=1,…,I−1,

denoting the restrictions from the top and the bottom , respectively (see Figure 1). It is convenient to identify and with functions and defined on by and , . Let where stands for the outer normal to , . Note that is an inward normal to on . In particular, and are top and bottom normals on , respectively. In the following, most quantities will be defined in terms of the bottom side.

We suppose that a body force acts on all of and surface forces act on the Neumann boundary . On the Dirichlet boundary the velocity of the displacement field of the deformable body is fixed at all time instants . We set on for convenience, though all further considerations can be generalized to the inhomogeneous case in a straightforward way.

We assume that the boundary forces are compressive in the sense that no fault opening occurs. This means that neighboring bodies and , , remain in contact throughout the evolution.

We consider a deformation field

 u=(u1,…,uI)∈H1(Ω1)d×⋯×H1(ΩI)d

where is the deformation of the subdomain . Throughout the paper we assume that the deformations within each subdomain are small, while the relative displacement between different subdomains can be large. Thus we will use a (geometrically) linear elastic approach inside of the subdomains while the coupling conditions have to take care of potentially large deformations.

Large deformation coupling conditions will be defined in terms of the deformed subdomains. Given the deformation fields the associated displacements are given by leading to the deformed subdomains . The actual contact surface of the deformed subdomains is then given by . In the following, we assume that each is injective, i.e. that each is regular enough to avoid self-penetration of . Furthermore, we assume that deformations are small, such that different surfaces do not get in contact after deformation. Then, the deformed contact boundary can be pulled back to the bottom and top reference domain according to

 ΓF,uB =(Id+uB)−1(Cu)⊂ΓF, ΓF,uT =(Id+uT)−1(Cu)⊂ΓF.

In the following, we will parameterize the top reference domain over the bottom one by the bijective contact mapping

 πu:ΓF,uB→ΓF,uT,πu=(Id+uT)−1∘(Id+uB),

which maps each bottom point to the unique top point , such that the corresponding displaced points and coincide. As a consequence, the deformed contact boundary can be parametrized both over using and over using .

Now, consider any piecewise defined scalar or vector field

 v=(v1,…,vI)∈H1(Ω1)k×⋯×H1(ΩI)k

with or . Then, we define the jump of across the deformed contact boundary on the contact reference domain according to

 (1) [v]u=vB−vT∘πuon ΓF,uB.

Contact conditions and friction laws will be stated in terms of normal and tangential components on the deformed contact boundary . To this end, let and denote by an outer normal to at the point , , i.e. is the pullback of an oriented normal field of the deformed contact boundary to using the map . Then we can decompose any vector field on according to its normal and tangential components with respect to the deformed configuration as

 v=vt+vnnu,vn=v⋅nu,vt=v−vnnu.

It is important to note that the tangential and normal component are defined in terms of the -dependent normal field that is defined (piecewise) on the deformed bottom subdomains and not with respect to the normal field of the reference subdomains .

We state a closed-fault condition (no penetration and no fault opening) by prescribing that the relative motion of the deformed subdomains and is tangential to the actual contact surface , i.e.

 (2) 0=[˙u]u⋅nu.

As a consequence, the jump of the relative tangential velocity satisfies

 [˙u]ut=[˙u]u−([˙u]u⋅nu)nu=[˙u]uon ΓF,uB.

The closed-fault condition is complemented by the balance of normal forces

 (3) (σ(u)n)B=−ωu(σ(u)n)T∘πu

on , where

denotes the stress tensor on

. Note that the normal force is a force per surface area, such that the change of the area element induced by the pullback to using has to be compensated by the weighting factor , while the minus sign compensates for the change of the normal direction on opposing sides. Note that the balance of normal forces (3) can alternatively be phrased as a jump condition with a transformed weighting factor .

Utilizing , we prescribe a rate- and state-dependent friction law of the form

 (4) −σt∈∂[˙u]uϕ([˙u]u,α)on ΓF,uB.

Here, we used the decomposition of the stress field on the bottom side into its normal and tangential components

 σn =σn(u)=(σ(u)n)B⋅nu, σt =σt(u)=(σ(u)n)B−σn(u)nu,

respectively, and denotes the subdifferential of a state-dependent convex functional to be described below. Note that the stress vector is computed with respect to the reference normal , while its decomposition in tangential and normal components is computed with respect to the deformed configuration with the corresponding normal . This reflects the fact that we assume small deformations within the subdomains while the relative deformations of subdomains can be large.

For given relative slip rate , the evolution of the state is given by

 (5) −˙α=∂αψ(α,∣∣[˙u]u∣∣)on ΓF,uB,−˙α=0on ΓF∖ΓF,uB

with a second convex functional . Note that the state remains constant on where no contact occurs.

Assuming a visco–elastic Kelvin–Voigt material law, and fixing some final time , we are now ready to state the following formal description of the deformation of a body with a layered fault system and rate- and state-dependent friction.

###### Problem 2.1 (Layered fault system with rate- and state dependent friction).

Find

 u:Ω×[0,T0]→Rdandα:ΓF×[0,T0]→R

such that

 (6) σ(u) =Aε(˙u)+Bε(u) in Ω∖ΓF (Kelvin–Voigt material) (7) divσ(u)+f =ρ¨u in Ω∖ΓF (balance of momentum) with boundary conditions, u=˙u =0 on ΓD (Dirichlet condition) σ(u)n =fN on ΓN (Neumann condition) frictional contact conditions, (10) [u]u⋅nu =0 on ΓF,uB (closed-fault condition) (11) (σ(u)n)B =−ωu(σ(u)n)T∘πu on ΓF,uB (balance of normal forces) (12) −σt ∈∂[˙u]uϕ([˙u]u,α) on ΓF,uB (state-dependent friction law) contact state condition, (14) −˙α ∈∂αψ(α,∣∣[˙u]u∣∣), on ΓF,uB (rate-dependent state law) and non-contact interface conditions (16) −˙α =0 on ΓF∖ΓF,uB (non-contact state condition) (17) (σ(u)n)B =0 on ΓF∖ΓF,uB (bottom Neumann condition) (18) (σ(u)n)T =0 on ΓF∖ΓF,uT (top Neumann condition)

holds for all . Here, is a constant material density, and stand for the viscosity and elasticity tensor, respectively, and is the linearized strain or strain rate tensor. In addition, we impose initial conditions on the displacement , velocity , and state .

Throughout the following, we assume that the tensor fields und have the symmetry properties

 Aijkl =Aklij, Aijkl =Ajikl, Bijkl =Bklij, Bijkl =Bjikl

such that the stress tensor and the bilinear forms induced by and are symmetric.

Note that Problem 2.1 provides an extension of the model presented in [30] that describes unilateral frictional contact of a deformable body with a rigid foundation. The tangential velocity relative to the fixed rigid foundation appearing in [30] is now replaced by the relative tangential velocity of adjacent deformable bodies.

A further extension to fault opening can be performed by replacing (10) by the non-penetration condition together with dynamical freezing and thawing of rate- and state-dependent friction (12), (14) in case of opening or closing faults.

For ease of notation we will skip the superscript and mostly write in the sequel.

### 2.2. The Dieterich–Ruina model

The current form of the Dieterich–Ruina model of rate- and state-dependent friction goes back to [37] (see also [1, 6, 11, 25, 33, 40] and the papers cited therein). It is based on the following ansatz for the friction coefficient

 (19) μ∗(V,θ)=μ0+alog(VV0)+blog(V0θL)

that depends on the rate and involves positive parameters , , , and .

It is complemented by a suitable evolution of the state . Here, most popular choices are

 (20) ˙θ =1−VLθ (Dieterich’s law) and (22) ˙θ =−VLθlog(VLθ) (Ruina’s law).

Following [30], we briefly sketch how this completely phenomenological friction model translates into a corresponding state-dependent friction law (12) and a rate-dependent state evolution (14) as postulated above. Starting from collinearity of relative tangential velocity and stress

 −σt|[˙u]|=[˙u]|σt|,

we postulate the state equation with normal stress to obtain

 (23) −σt=|σn|μ∗(V,θ)[˙u]|[˙u]|.

In analogy to Tresca friction, we now replace the solution dependent normal stress by a given parameter  [30].

As becomes negative and thus meaningless for

 0≤V

we replace by its regularization

 μ(V,θ)={μ∗(V,θ) if V≥Vm(θ)0 otherwise .

Then elementary calculations show that (23) takes the form (12) with the convex functional defined by

 (24) ϕ∗(v,θ)={a|¯σn|(|v|log(|v|/Vm(θ))−|v|+Vm(θ)) if |v|≥Vm(θ)0 otherwise .

It remains to show that the rate evolutions (20) and (22) can be rewritten according to (14). Introducing the transformed state , Dieterich’s law (20) takes the form (14) with the scalar convex function

 (25) ψDieterich(α,V)=VLα+e−α.

Ruina’s law (22) is recovered in terms of (14) by the same transformation and the scalar convex function

 (26) ψRuina(α,V)=VL(12α2+log(V/L)α).

Inserting the transformed state into (24), we obtain the corresponding rate functional

 (27) ϕ(⋅,α)=ϕ∗(⋅,eα).

### 2.3. Weak formulation

We consider the Hilbert space with the canonical inner product , , , and introduce the closed linear subspace

 H0={v∈H|v=0 on ΓD}

of admissible displacements respecting the Dirichlet boundary conditions. The normal jump condition is incorporated in the closed affine subspace

 (28) Hu0={v∈H0|[v]u⋅nu=0}.

With the tensors , taken from (6), we introduce the bilinear forms

 (29) a(v,w)=∫Ω∖ΓFAε(v):ε(w)dx,b(v,w)=∫Ω∖ΓFBε(v):ε(w)dx,v,w∈H0,

involving the linear strain tensor , together with the linear functional

 (30) ℓ(v)=∫Ωfvdx+∫ΓNfNvds,v∈H0.

To ensure that the bilinear forms are well defined, we assume that the tensor fields und are uniformly elliptic in the sense that the bilinear forms induced by and on the space of symmetric matrices are elliptic with constants independent of .

By inserting the stress–strain relation (6) into the balance of momentum (7), testing with , integrating by parts, and exploiting the symmetry of together with the boundary conditions on and we then formally obtain

 (31) ⟨ρ¨u,v−˙u⟩+a(˙u,v−˙u)+b(u,v−˙u)−ℓ(v−˙u)=∫ΓF(σ(u)n)B⋅(v−˙u)Bds+∫ΓF(σ(u)n)T⋅(v−˙u)Tds

for all and . Here stands for the pairing of with its dual . Using the boundary conditions (17) and (18), integral transformation from to , and the normal force balance (11) we can rewrite the right hand side in (31) as

 ∫ΓF,uB(σ(u)n)B⋅(v−˙u)Bds+∫ΓF,uT(σ(u)n)T⋅(v−˙u)Tds = ∫ΓF,uB(σ(u)n)B⋅([v−˙u]u+(v−˙u)T∘πu)ds+∫ΓF,uBωu((σ(u)n)T⋅(v−˙u)T)∘πuds = ∫ΓF,uB(σ(u)n)B⋅[v−˙u]uds.

For given state , we introduce the convex functional on according to

 (32) Φu(⋅,α)=∫ΓF,uBϕ([⋅]u,α)ds

with the convex functional taken from the friction law (12). Now let satisfy the closed-fault condition (10), i.e. . Then, utilizing the decomposition together with the friction law (12), the closed-fault condition on , and the definition of subdifferentials, we find that

 (33)

Now we insert (33) into (31), in order to obtain the desired weak form of the rate equation

 (34) ⟨ρ¨u,v−˙u⟩+a(˙u,v−˙u)+b(u,v−˙u)+Φu(v,α)−Φu(˙u,α)≥ℓ(v−˙u)∀v∈Hu0.

Similarly, for given velocity and thus given rate , we define the convex functional on by

 (35) Ψu(⋅,˙u)=∫ΓF,uBψ(⋅,∣∣[˙u]u∣∣)ds

with the convex functional taken from the state law (14), and test the state evolution (14) with to obtain the weak formulation

 (36) (˙α,β−α)L2(ΓF)+Ψu(β,˙u)−Ψu(α,˙u)≥0∀β∈L2(ΓF).

This formulation automatically satisfies the non-evolution condition (16) for the state on the non-contact boundary since is defined on but only depends on values of on the contact boundary .

We are now ready to state the weak formulation of Problem 2.1

###### Problem 2.2 (Weak formulation).

Find

 u∈H1((0,T0),H0)∩H2((0,T0),H∗0) and α∈H1((0,T0),L2(ΓF))

such that and

 (37) ⟨ρ¨u,v−˙u⟩+a(˙u,v−˙u)+b(u,v−˙u)+Φu(v,α)−Φu(˙u,α) ≥ℓ(v−˙u) ∀ v∈Hu0, (38) (˙α,β−α)L2(ΓF)+Ψu(β,˙u)−Ψu(α,˙u) ≥0 ∀ β∈L2(ΓF)

holds for almost all together with initial conditions

 (39) u(0)=u0,˙u(0)=˙u0,α(0)=α0

with given and and .

It is natural to start the evolution out of an equilibrium configuration, i.e., with an initial displacement that solves the stationary problem

 (40) u0∈Hu00:b(u0,v)=ℓ(v)∀v∈Hu00.

In our numerical experiments to be reported below, (40) is solved iteratively by a fixed point iteration over the geometric nonlinearity, i.e., starting with , a new iterate is computed as the (up to tangential rigid body motions) unique solution of the corresponding linear problem on .

To our knowledge, existence and uniqueness of solutions of Problem 2.2 is widely open. In case of unilateral frictional contact with a rigid foundation and Dieterich’s law (20), long-time existence of solutions was established by Pipping [28].

## 3. Semi-discretization in time

Utilizing Rothe’s method [19, 36], we first perform a time discretization of Problem 2.2 leading to a sequence of continuous spatial problems to be (approximately) solved in each time step. To this end, the time interval is partitioned into time steps with given step size , and we write for notational convenience.

### 3.1. Rate problem with given state

We first consider the rate problem (34) for given state . Following [30], we apply the classical Newmark scheme

 (41) ˙un=˙un−1+τ2(¨un−1+¨un)un=un−1+τ˙un−1+(τ2)2(¨un−1+¨un),n=1,…,N,

which is well-known to be energy-conserving, consistent with second order, and unconditionally stable [17]. Utilizing (41), we eliminate

 (42) ¨un=2τ(˙un−˙un−1)−¨un−1,un=un−1+τ2(˙un+˙un−1),,n=1,…,N,

from (34) at fixed time and freeze the solution dependence in the closed-fault condition and in the friction law at to obtain the spatial variational inequality

 (43) ˙un∈Hun−10:an(˙un,v−˙un)+Φun−1(v,α)−Φun−1(˙un,α)≥ℓn(v−˙un),∀v∈Hun−10,

for . Here, we have set

 (44) an(v,w)=2τ(ρv,w)+a(v,w)+τ2b(v,w)

with denoting the canonical scalar product in and

 ℓn(v)=ℓ(v)+(ρ¨un−1,v)+2τ(ρ˙un−1,v)−τ2b(˙un−1,v)−b(un−1,v).

Note that is not given as an initial condition in the continuous Problem 2.2. Assuming initial acceleration towards equilibrium, is therefore computed from the auxiliar problem

 (45) ¨u0∈H0:(ρ¨u0,v)+b(u0,v)=ℓ(v)∀v∈H0.

Note that the jump terms , the contact boundary , and the contact mapping are all taken with respect to the last deformed state of the contact boundary. This eliminates the geometric nonlinearity associated with large (relative) deformations of the contact boundary, and we are left with the variational inequality (43) on the affine subspace of to be solved in each time step..

As is symmetric and positive definite and is convex, the variational inequality (43) can be equivalently written as the minimization problem

 (46) ˙un∈Hun−10:J(˙u,α)≤J(v,α)∀v∈Hun−10

for the corresponding energy functional

 J(v,α)=12an(v,v)+Φun−1(v,α)−ℓn(v).

The following lemma [12, Theorem 6.49] will be useful to show existence and uniqueness of a solution.

###### Lemma 3.1.

Assume that , , is a non-negative function, such that is lower semicontinuous for almost all . Then the induced functional

 ∫ΓFg(x,⋅)dx:L2(ΓF)→R∪{+∞}

is lower semicontinuous.

The convex functional defined in (32) for the Dieterich–Ruina model (27) is proper and lower semicontinuous by Lemma 3.1. Furthermore, by the assumptions on and , the bilinear form on the symmetric matrices is symmetric and uniformly elliptic with respect to . The following existence result therefore follows from Korn’s second inequality and [13, Lemma 4.1].

###### Proposition 3.2.

Let and . Assume that , , avoids self-penetration so that the contact mapping and thus are well-defined. Then the spatial rate problem (43) has a unique solution and any given state .

As a consequence of Proposition 3.2, the solution operator ,

 (47) L2(ΓF)∋α↦R(α)=˙un∈Hun−10,

of the spatial rate problem (43) is well-defined, if no self-penetration occurs in preceding time steps. This is a strong assumption, as the contact conditions are taken explicitly.

### 3.2. State problem with given rate

Discretizing the state problem (38) with given velocity by the implicit Euler method and freezing the state law at yields the variational inequality

 (48) αn∈L2(ΓF):(αn,β−αn)L2(ΓF)+τΨun−1(β,˙u)−τΨun−1(αn,˙u)≥(αn−1,β−αn)L2(ΓF)∀β∈