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 ) by numerical simulations.
In the underlying mathematical description, the Dieterich–Ruina model of rate- and state-dependent friction (RSF)  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. 
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  developed a numerical method for the dynamic simulation of slip events. This method was later generalized to three space dimensions  and cast into the software package SeisSol that was successfully utilized for the simulation, e.g., of the 2016 Kaikōura earthquake cascade . 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 . 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  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 . 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 , ,
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
For with , , we define the restrictions and of to with
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
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
In the following, we will parameterize the top reference domain over the bottom one by the bijective contact mapping
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
with or . Then, we define the jump of across the deformed contact boundary on the contact reference domain according to
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
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.
As a consequence, the jump of the relative tangential velocity satisfies
The closed-fault condition is complemented by the balance of normal forces
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
Here, we used the decomposition of the stress field on the bottom side into its normal and tangential components
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
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).
|(7)||in||(balance of momentum)|
|with boundary conditions,|
|frictional contact conditions,|
|(11)||(balance of normal forces)|
|(12)||(state-dependent friction law)|
|contact state condition,|
|(14)||(rate-dependent state law)|
|and non-contact interface conditions|
|(16)||(non-contact state condition)|
|(17)||(bottom Neumann condition)|
|(18)||(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
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  that describes unilateral frictional contact of a deformable body with a rigid foundation. The tangential velocity relative to the fixed rigid foundation appearing in  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  (see also [1, 6, 11, 25, 33, 40] and the papers cited therein). It is based on the following ansatz for the friction coefficient
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
Following , 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
we postulate the state equation with normal stress to obtain
In analogy to Tresca friction, we now replace the solution dependent normal stress by a given parameter .
As becomes negative and thus meaningless for
we replace by its regularization
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
Inserting the transformed state into (24), we obtain the corresponding rate functional
2.3. Weak formulation
We consider the Hilbert space with the canonical inner product , , , and introduce the closed linear subspace
of admissible displacements respecting the Dirichlet boundary conditions. The normal jump condition is incorporated in the closed affine subspace
With the tensors , taken from (6), we introduce the bilinear forms
involving the linear strain tensor , together with the linear functional
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
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
For given state , we introduce the convex functional on according to
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
Similarly, for given velocity and thus given rate , we define the convex functional on by
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).
such that and
holds for almost all together with initial conditions
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
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 .
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
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
for . Here, we have set
with denoting the canonical scalar product in and
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
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
for the corresponding energy functional
The following lemma [12, Theorem 6.49] will be useful to show existence and uniqueness of a solution.
Assume that , , is a non-negative function, such that is lower semicontinuous for almost all . Then the induced functional
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].
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 .
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