Finite volume discretization for poroelastic media with fractures modeled by contact mechanics

A fractured poroelastic body is considered where the opening of the fractures is governed by a nonpenetration law while slip is described by a Coulomb-type friction law. This physical model results in a nonlinear variational inequality problem. The variational inequality is rewritten as a complimentary function, and a semismooth Newton method is used to solve the system of equations. For the discretization, we use a hybrid scheme where the displacements are given in terms of degrees of freedom per element, and an additional Lagrange multiplier representing the traction is added on the fracture faces. The novelty of our method comes from combining the Lagrange multiplier from the hybrid scheme with a finite volume discretization of the poroelastic Biot equation, which allows us to directly impose the inequality constraints on each subface. The convergence of the method is studied for several challenging geometries in 2d and 3d, showing that the convergence rates of the finite volume scheme do not deteriorate when it is coupled to the Lagrange multipliers. Our method is especially attractive for the poroelastic problem because it allows for a straightforward coupling between the matrix deformation, contact conditions, and fluid pressure.


Introduction
Slip and opening of fractures due to fluid injection is of relevance to a number of subsurface engineering processes. In hydraulic reservoir stimulation, the effect is deliberately induced, while in storage operations and wastewater disposal, avoiding reactivation and opening of fractures is important for preserving caprock integrity. In any circumstance, triggering of larger slip events in the form of elevated levels of seismicity must be avoided. The mathematical model of fracture resistance, slip and opening results in a strongly coupled nonlinear variational inequality, which requires advanced numerical schemes to solve. The purpose of this work is to describe and implement a numerical method to solve this problem considering a poroelastic matrix. The fractures are a set of predefined surfaces in the domain, and the nucleation or growth of fractures is not considered.
The flow and mechanics of poroelastic media and the contact mechanics of elastic bodies are well-developed research fields. For a porous or poroelastic medium, we refer to the classical textbooks [1,2]. There exists an extensive number of discretizations for the elliptic equations describing fluid flow in a porous medium, and they all have different merits. The most popular discretizations are the so-called locally conservative discretizations [3], which include mixed finite elements [4], control-volume finite elements [5], and finite volume methods [6]. For the coupled poroelastic problem, without considering fractures, it is known that a naive discretization of the coupling terms of the fluid pressure and the solid displacement leads to stability issues for finite element schemes [7]. Several different methods have been proposed to remove these oscillations [8,9,10]. Recently, a finite volume method called the multipoint stress approximation was introduced for elastic deformations [11,12]. This method has been extended to the poroelastic Biot equations and shown to be stable without adding any artificial stabilization terms in the limit of incompressible fluids and small time steps [13].
The contact mechanics problem, i.e., contact between two elastic bodies, is also the topic of several textbooks [14,15]. A widely used solution strategy for the nonlinear variational inequalities resulting from the mathematical formulation is the penalty method [16]. The basic idea is to penalize a violation of the inequality by adding extra energy to the system. The solution depends then, in a very sensitive way, on the choice of the penalty parameter. If the value of the parameter is too small, the condition number of the algebraic system is extremely poor, and the nonlinear solver converges slowly. If the value is too large, the accuracy of the solution is very poor, and unphysical approximations can be obtained. Therefore, variationally consistent hybrid formulations have gained interest recently. The hybrid formulations are based on the discretization of additional unknown Lagrange multipliers added to the contact region. This method has been applied to, among others, the Signorini problem [17], frictional contact [18], and large deformations [19]; see the survey contribution [20] and the references therein.
For a poroelastic domain including fractures, different models for the contact problem are developed [21,22,23,24]. Most of these models, however, do not take into account the contact problem either by assuming the fractures stick together [21] or that the fluid pressure inside the fractures is so large that the fracture surfaces are never in contact [22,23]. The full contact problem for a fractured poroelastic domain is considered by Garipov et al [24], where they applied the penalty method to solve the nonlinear variational inequalities resulting from the contact problem.
In the current work, we present a different numerical solution approach for poroelastic media with contact mechanics. The discretization is based on a finite volume method for poroelasticity [13] combined with a variationally consistent hybrid discretization [25,20]. The hybrid formulation considered in this work can be regarded as a mortar formulation [26] using matching meshes with the displacement as the primal variable and the surface traction as the dual variable. The finite volume scheme has previously been extended to fracture deformation by adding additional displacement unknowns on the fracture faces [27]. This formulation was successfully used to implement a fixed-point type iteration to approximate the friction bound [28]; however, this formulation suffers from the fact that a step length parameter needs to be tuned and that it might require many iterations to converge [29]. An advantage of the scheme used in this work, where the Lagrange multiplier of the hybrid formulation is coupled with the surface traction obtained from the finite volume scheme, is that it gives a natural formulation of the contact condition per subface. This formulation allows us to rapidly solve the resulting nonlinear inequality problem by applying a semismooth Newton method; see the work by Hüeber et al [25], among others [20,30].
The remainder of this paper is structured as follows. First, we state the problem and give the governing equations. Then, the discretization is presented, which is divided into two parts: (i) the finite volume discretization for the Biot equations and (ii) the discrete hybrid formulation for the contact problem. We present three numerical examples. The first two consider the dry case where the coupling between fluid pressure and deformation of the rock is disregarded. The last example solves the poroelastic deformation of a 3d domain where the deformation of fractures is governed by a Coulomb friction law. Finally, we give concluding remarks.

Problem statement
Let Ω be a fractured deformable porous body. The boundaries of the domain ∂Ω are divided into three disjoint open sets, ΓD, ΓN , and ΓC , as illustrated in Figure 1: for the first set, a Dirichlet boundary condition is assigned; for the second, a Neumann boundary condition is assigned; and the last is the internal fracture boundary. We consider the Biot model for a poroelastic medium [31]: in Ω, where we have assumed backward Euler time stepping for a given time step η. All parameters are, in general, functions of space, e.g., C = C(x), x ∈ Ω; however, the explicit dependence is suppressed to keep the notation simple. Parameters associated with the pressure p and displacement u are given a subscript with the same symbol. The vector fu is a given body force, while fp is a given source term including the contribution from the previous time step. The stiffness tensor is denoted C, the Biot coupling coefficient α, the storage coefficient c0, and the permeability K. Indicated by subscripts, g represents Dirichlet and Neumann boundary conditions for displacement and pressure, while n is the outwards pointing normal vector. In this work, we use C : (∇u + (∇u) )/2 = G(∇u + (∇u) ) + Λtr(∇u)I, where G and Λ are the Lamé parameters. Traction can also be derived for other material laws. The fracture boundary, ΓC , is divided into a positive side Γ + and a negative side Γ + . The choice of which side is positive and which is negative is arbitrary and will only make a difference in the implementation. For Figure 1: A domain Ω where the external boundary is divided into two parts: Γ D and Γ N . Included in the domain are two internal boundaries, or fractures, Γ C . The two sides of the internal boundaries are labeled Γ + and Γ + , as shown in the magnified circular region of the domain. The function g(x), x ∈ Γ + gives the initial gap between the two fracture sides. The left fracture has an initial gap g > 0, while the top right fracture has an initial gap g = 0.
the fracture segments, a nonpenetration condition is enforced in the normal direction, meaning that the positive and negative sides cannot penetrate each other. In the tangential direction, a Coulomb friction law divides the contact region into a sliding part and a sticking part. To formulate these contact conditions, the normal vector for the contact region is defined as the normal vector of the positive side n(x) = n + (x). Further, let be a mapping that projects a point from the positive boundary onto the negative boundary as given by the normal vector. The gap function, which will appear in the nonpenetration condition, is then defined as where · is the Euclidean norm. Due to Newton's third law, the surface traction, T = σ · n, on the contact boundaries must be equal up to the sign and we use the notation TC = T + . The surface traction is divided into a normal and tangential part by and the displacement jump is defined as [u(x)] = u(x) − u(R(x)) for x ∈ Γ + . The normal and tangential displacement jump is defined analogously to Equation (4): The nonpenetration condition can now be formulated as where the first condition ensures that the two sides of the fracture cannot penetrate, the second ensures that either the normal traction is zero or the fracture sides are in contact, and the last enforces a negative normal component of the surface traction.
x K Figure 2: Notation used to describe the mesh. For a cell K, face π and vertex v of the mesh, we associate a subcell (K, v) and subface (π, v), as well as a cell center x K and continuity point x v π .
The tangential part of the surface traction is governed by a Coulomb friction law: where F is the coefficient of friction, andu is the displacement velocity. The first equation gives the friction bound, the second ensures that if the friction bound is not reached, then the surface is sticking, and the last equation ensures that if the friction bound is reached, then the tangential sliding velocity is parallel to the tangential traction. In the static case, e.g., for the purely mechanical problem when α = 0, the notion of a velocity does not exist. For these cases, it is common to replace the sliding velocity, [u]τ , by the displacement jump, [u]τ , in Equation (6) [20]. In the formulation of Biot's equations given in Equation (1), a backward Euler time-stepping is used, and the velocity at time step n + 1 is approximated byu n+1 = (u n+1 − u n )/η. The sliding in the dynamic case can therefore be considered as the static formulation with an additional (known) contribution from the previous time step. When the discretization is presented in the following section, the static case is used to avoid an additional index for the time stepping. However, this is easily extended to the dynamic case by inserting the discrete velocity (u n+1 − u n )/η. For the fluid, the fractures are modeled as impermeable. This means that the fluid cannot flow in or across the fractures, i.e., q(x) · n(x) = 0, x ∈ ΓC .

Discretization
We define the triplet (T , F, V) as the cells, faces and vertices of our mesh. Before the discretization is described, we need to define some notation, and we start by giving the relation between cells, faces and vertices using the standard notation for finite volume methods [32,13]: • For a cell K ∈ T , we denote its faces by FK and its vertices by VK .
• For a face π ∈ F , we denote the neighboring cells as Tπ and its vertices as Vπ.
• For a vertex v ∈ V, we denote the adjacent cells by Tv and the faces by Fv.
In addition to the mesh triplet (T , F, V), we define the so-called subcells and subfaces illustrated in Figure 2: • For a vertex v ∈ Vπ, we associate a subface identified by (π, v) with an area m v π such that v∈Vπ m v π = mπ = π dx.
Note that in an abuse of notation, we use K for both indexing and the geometric object so that both VK and K dx make sense. All subfaces (π, v), π ∈ F , v ∈ Vπ are divided into three disjoint sets P, N , and R, where P contains all subfaces located on the positive boundary Γ + , N contains all subfaces located on the negative boundary Γ + , and R contains the remaining subfaces.
Finally, for each element K ∈ T , a cell center xK ∈ K is defined, and for each subface (π, v), we associate a continuity point x v π located somewhere between the vertex v and the face center of π. The unit normal for each face is denoted by nπ, which is equal to the subface normal of the face n v π . When it is necessary to distinguish the direction of the normal, it is defined as the outward pointing normal n π K of a cell K ∈ Tπ. Note that for a face π, we have Tπ = {K, L}, n π K = −n π L .

Finite volume discretization
We use a finite volume discretization [13] to discretize the Biot Equations (1). This is based on two discrete variables, uK and pK , which are the cell-centered displacement and pressure, respectively. Within each subcell (K, v), K ∈ T , v ∈ VK , it is assumed that the displacements and pressures are piecewise linear, and the discrete gradients are denoted by (∇u) v K and (∇p) v K , where the bar over the gradient operator is added to distinguish it from the continuous gradients. For the mechanical stress, we adapt the notion of weak symmetry [12]; given the volume weighted average associated with a vertex v, the discrete weakly symmetric mechanical stress is given by This is referred to as weak symmetry because To simplify notation, the tensor C v K is referred to as the stress tensor, which acts to weakly symmetrize the stress: θ should not be interpreted as a single tensor vector product but as a weighted sum of products given by Equation (7).
Using the weak symmetry, the flux and traction over each subface given by the discrete variables can be stated as For a spatially varying permeability and stress tensor, we use the cell-center value to evaluate the parameters KK = K(xK ) for each cell. The finite volume scheme will be constructed such that the gradient unknowns can be eliminated by performing a local static condensation. In the following, we first present the local systems that are used to do so. The finite volume structure of our discretization is then obtained by enforcing mass/momentum conservation for each cell. The final scheme will be locally conservative and given by the cellcentered displacement and pressure. A detail that will be important when we introduce the hybrid discretization is the possibility of exactly reconstructing the discrete gradients, and thus also the flux and traction, from the cell-centered variables u and p.
The discrete fluid flux given in Equation (8) does not contain any dependence on the displacement, and it is identical to the fluid flux for the uncoupled fluid pressure, i.e., α = 0. To discretize the flux, we can use any of the standard methods for flow in porous media, and we have chosen the well-known MPFA scheme [6]. Each subcell gradient (∇p) v K is associated with a fluid flux as given in (8). Conservation of mass is enforced for each subface (π, v), π ∈ F, v ∈ Vπ. This requires the fluid flux for cells (K, L) ∈ Tπ sharing a face π to be equal and opposite over each of their shared subfaces; that is, The pressure is not required to be continuous across the whole subface. Instead, pressure continuity is enforced at the continuity points, x v π , that is, Here, we have made use of the assumption that the pressure is linear in each subcell to write the pressure at the continuity point x v π as a function of the cell center pressure pK and gradient (∇p) v K . For subfaces on the boundary, the right-hand side of Equation (10) and Equation (11) is replaced by the Neumann and Dirichlet conditions, respectively. Around each vertex v we can now form a local linear system of equations from which the gradients (∇p) v K , K ∈ Tv can be eliminated: The first block Qp(∇p)v = gp,N in this linear system is the collection of all flux balance Equations (10) for the vertex v. The next block Dp,G(∇p)v = gp,D − Dpp collects all the pressure continuity Equations (11). Thus, (∇p)v is the vector of the subcell gradients (∇p) v K , the matrix Qp represents products of the form m v π n π K KK , the matrix Dp,G represents the distances x v π − xK , the vectors gp,N and gp,D are possible boundary conditions, and Dp has entries 1 for p v K and −1 for p v L . The elimination of the displacement gradients (∇u) v K is similar to the elimination of the pressure gradients ∇p v K . First, the continuity of traction gives us for each subface It is worth pointing out that, for internal faces, the averaging part of the operator C v K : (∇u) v K is the same on the right-and left-hand sides. Thus, the balance of traction can be written as However, for boundary faces, the complete Equation (13) must be used. Unlike the fluid fluxes in (8), the traction is different from the uncoupled system due to the term αpK I. It is important to include the Biot stress in the local systems to obtain the correct force balance in our method [13]. We will see later that this approach also gives a higher-order term in the mass balance for the fluid, which acts analogously to the stabilization terms in other colocated schemes. For the fluid pressure, displacement continuity is enforced at the continuity points Again, Neumann and Dirichlet subfaces are incorporated by changing the right-hand side of Equation (13) and (14), respectively. A local elimination of the displacement gradients (∇u) v K can now be done around each vertex to express them in terms of the cell-center displacement and pressure: For the pressure gradients, (∇u)v is the vector of the displacement gradients, (∇u) v K around vertex v, the matrix Qu represents products of the form m v π n C v K , the matrix Du,G represents the same distance vectors as in (12), the vectors gu,N and gu,D are possible boundary conditions, and Du is a matrix with entries ±1. The term P is the only difference between the coupled and uncoupled system and contains products of the form m v π αIn v π . The finite volume discretization will conserve fluid mass over each cell K, The pressure gradient (∇p) v K and displacement divergence (∇ · u) v K = tr(∇u) v K are obtained as linear functions of the cell-centered pressures and displacements from the local systems given in (12) and (15). The appearance of the pressure in the discrete displacement divergence is essential for the consistency of the method and is similar to the artificially introduced stability terms in other methods; see, e.g., Gaspar et al [33].
For the mechanics, momentum is conserved for all cells K, Note that the term αpK I from the Biot stress in (9) sums to zero over a cell due to Gauss's theorem; however, the pressure dependence on the subscale gradients gives the correct fluid pressure contribution to the mechanics. To summarize, the finite volume scheme is constructed by defining a set of discrete pressure and displacement gradients for each subcell. Flux and pressure continuity is enforced over each subface for the fluid, and traction and displacement continuity is enforced for each subface for the solid. This defines a small local system around each node from which the pressure and displacement gradients can be expressed as a linear combination of the cell-centered pressure and displacements and then eliminated. A stable coupling between the fluid and solid is achieved by considering the Biot stress, i.e., C : (∇u + (∇u) )/2 − αpI, for traction balance of the local systems.

Hybrid formulation
To solve the contact conditions (5) and (6), we apply the active-set strategy, which is equivalent to a semismooth Newton method described by Hüeber et al [25]. See also the paper by Wohlmuth [20]. We recapitulate the solution strategy in this section for the completeness of this paper. The main difference in our approach is how the Lagrange multipliers, which represent the surface traction, are coupled to the displacement unknowns in the surrounding domain. In our finite volume scheme, the Lagrange multipliers will enter into the local equations for the displacement gradients.
A set of Lagrange multipliers is defined on the positive subface boundaries The normal λ v πn and tangential λ λ λ v πτ components of the Lagrange multiplier are defined analogously to (4). The displacement u v π is obtained as in Equation (14) for local systems. The discrete formulation of the nonpenetration condition (5) can for each subface be written as and the static Coulomb friction (6) as Recall that for the dynamic case, the displacement jump is replaced by the sliding velocity, [u v π ]τ . The nonpenetration condition can be rewritten as the nonlinear complementary function and similarly for the Coulomb friction Here c > 0 is a given numerical parameter and max{·, ·} is the maximum function. The solution pair (u v π , λ λ λ v π ) satisfies the nonpenetrating condition (17) and Coulomb law (18) if and only if Cn([u v πn ], λ v πn ) = 0 and Cτ ([u v π ], λ λ λ f ) = 0. We apply a semismooth Newton method to this problem, which results in an active set method. Given the solution (u k , λ λ λ k ) from the previous Newton iteration, we define b v,k ) and divide the contact subfaces into three disjoint sets: The first set contains the subfaces not in contact. The second set contains the subfaces in contact whose friction bound is not reached, i.e., they are sticking. The third set contains the subfaces in contact where the friction bound is reached, i.e., they are sliding. The new iterates (u v,k+1 πτ , λ λ λ v,k+1 πτ ) in the semismooth Newton scheme are then calculated depending on which set the subface belongs to. The update is found by calculating the derivative of the complementary functions Cn and Cτ for each of the three sets. For the subfaces not in contact, zero traction is enforced λ λ λ v,k+1 For the subfaces in contact and sticking, we enforce This enforces the condition that the negative and positive subfaces must be in contact in the next iteration k + 1. In the tangential direction the normal component of the Lagrange multiplier only adds a contribution if the subface was sliding in the previous iteration, k. For subfaces in contact and sliding, we enforce Again, this enforces the condition that the negative and positive subfaces be in contact at the next iteration k + 1. In the tangential direction, the condition approximates the sliding direction and distance. The matrices and vectors are: .

Regularization
For the subfaces in the inactive set I k+1 n , i.e., the subfaces not in contact, the Newton update gives a homogeneous Neumann boundary condition. For the subfaces in the contact sets I k+1 τ and A k+1 , the Newton update gives a Dirichlet condition in the normal direction and a Robin boundary condition in the tangential direction. This Robin condition guarantees positive definiteness of the system only if L v,k π , defined by Equation (25), is positive definite. In the converged limit, the matrix L v,k π tends to a positive definite matrix. However, during the iterations, there is no guarantee that this will hold. We therefore add a regularization to the Robin conditions by replacing , which is only different from Q v,k π when the inequalities in Equation (18) are violated. Further, we define otherwise.
Using the notation that tilde (·) denotes the regularization, we haveM v,k π = e v,k π (I d−1 −Q v,k π ) and replace the matrix L v,k π from Equation (25) which guarantees its positive definiteness.

Contribution to local system
The Robin condition appears in the local system for the displacement gradient by replacing the displacement continuity Equation (14) by where the matrices A v,k π andL v,k π and vectorr v,k π are defined according to Equation (22)- (24). The Lagrange multiplier also appears on the right-hand side for the traction balance in Equation (13): The contribution to the negative side −λ λ λ(R −1 (x v π )) is just the mapping onto the Lagrange multiplier on the corresponding positive subface as given by Equation (2). The local systems for the displacement gradients are then written in the form where the matrices A † and L † and the vector r † are collections of the matrices A v,k π andL v,k π and the vectorr from Equation (26). The matrix M± contains the positive areas m v π for the positive subfaces and the negative areas −m v π for the negative subfaces and represents the Lagrange multiplier contribution to the traction balance in Equation (27). The matrix D λ,G represents the distances from cell centers to the continuity point, and D λ contains 1 for the positive subfaces and −1 for the negative subfaces. These two matrices are equivalent to Du,G and DU and can be interpreted as the displacement at the continuity points. From a computational point of view, it is worth noting that during the Newton iteration, only the matrices L † , A † and the vector r † will change. This means that updating the discretization is inexpensive as it is only a local update for the subcells bordering the fractures.

Numerical examples
Three numerical examples are given. For the first two, we neglect the fluid contribution to the mechanical stress to investigate the performance of the numerical approach for the purely mechanical contact problem, i.e., we set α = 0. In all of the examples, Young's modulus is E0 = 4 GPa, the Poisson ratio is ν = 0.2, and the initial gap of the fractures is g = 0. In our experience, the algorithm is quite robust with respect to the numerical parameter c, and in the examples, it is fixed to c = 100 GPa/m.
We assign a space varying coefficient of friction so that the slip of the fractures will arrest before it reaches the fracture tips. This choice of the friction coefficient is done to obtain a solution with high enough regularity to study the convergence in stress. If the slip of the fractures reaches the fracture tips, the solution will contain singularities in the stress, which reduces the regularity of our solution. Note that our method is not restricted to the regularized solution, as discussed more thoroughly in Appendix A.
The discrete solution is denoted u h , which is interpreted as the piecewise constant function over each cell K x such that u h (x) = uK . The discrete solution λ λ λ h for the Lagrange multiplier is defined as piecewise constant on each face π on Γ + and is equal to the area weighted sum of the subface values, λ λ λ h (x)mπ = v∈Vπ m v π λ v π , x ∈ π. Likewise, the displacement jump is defined as the piecewise constant on each face, π, on Γ + corresponding to the subface average, where |Vπ| is the number of subfaces of the face, π, which is equal to three if π is a triangle. The continuous solution is denoted by the pair (u, λ λ λ).
We define the relative error of a discrete variable ξ h in a domain γ as where ξ is a reference solution. The Newton iteration is terminated when the change in the solution is below a given stopping criterion: where k is the Newton iteration index. The implementation has been implemented in the open source Python toolbox Porepy [34], which has an interface for meshing in Gmsh [35]. The run scripts for the examples are open source [36]. ParaView [37] was used to make

Example 1
The first example is a domain 2 m × 1 m with six fractures, as depicted in Figure 3. This example includes difficult cases such as a fracture with a kink and a fracture reaching the boundary. An advantage of our finite volume method is that no special treatment is needed to handle these cases because the degrees of freedom are located in the cell and subface centers and not on the nodes. In this example, we do not consider any fluid and solve only for the linear elasticity in Equation (1) coupled to the contact conditions given in Equations (5) and (6). The bottom boundary is fixed, the two vertical boundaries are free, and at the top boundary a Dirichlet condition gu,D = [0.005, −0.002] m is assigned. The initial guess in the Newton iteration is u = 0 m, λn = −100 Pa and λ λ λτ = 0 Pa, i.e., zero displacement and all fractures in contact and sticking. The coefficient of friction is for each fracture i = 1 . . . 6 set to Fi(x) = 0.5(1 + exp(−Di(x) 2 /0.005m 2 ), x ∈ Γ + i , where Di(x) is the distance from x to the tips of fracture i. Note that the bend in Fracture 1 and the right end of Fracture 5 are not considered tips for the distance function D, and thus the coefficient of friction at these points is F ≈ 0.5.
A contour plot of the solution is shown in Figure 3 where the discontinuous displacement over the sliding or opening fractures can clearly be seen. To better visualize the different behaviors of the fractures, the fracture regions that are slipping, sticking, and opening are plotted in different colors in Figure 3. For Fracture 1, the top boundary is sliding to the right, while the bottom boundary is sliding to the left. This situation causes the fracture to open in a small segment after the bend. Figure 4 shows the shear component of the Lagrange multiplier as well as the friction bound and displacement jump. At the bend of fracture 1, there is a singularity in the stress that causes the sharp increase in the Lagrange multiplier. For Fracture 2, we observe a change in the shear and normal component of the Lagrange multiplier at approximately the midpoint that is caused by the opening of Fracture 6. In the vicinity of the fracture tips, there is a sharp increase in the shear component of the Lagrange multiplier as the fractures change behavior from sliding to sticking.
As a reference solution (u, λ λ λ), we use the solution calculated for a fine mesh using 1.7 million degrees of freedom. The second finest mesh has 270 thousand degrees of freedom and is the mesh used for the results in Figure 4 and 3 . In Figure 5, the relative errors [u]) and ε Γ + i (λ λ λ h , λ λ λ), given by Equation (29), are plotted for each fracture i = 1, . . . , 6. For the displacement jump, the convergence is of first order for all fractures except Fracture 4, which is correctly predicted to be sticking (and thus, the error is zero). For the Lagrange multiplier λ λ λ h , we observe first-order convergence for Fractures 4 and 5, while the error for Fracture 6 is zero.   Table 1: Fracture geometry in Example 2 and 3. The strike angle is the rotation from x-axis in the x-y-plane defining the strike line. The dip angle is rotation around the strike line.
[0, 0, −4.5] MPa. The coefficient of friction is for the two fractures, i = 1, 2: where Ri is the radius of fracture i and Di(x) the distance from the center of the fracture to x. Figure 7 shows the displacement jump [u h ]τ and the shear component of the Lagrange multiplier λ λ λ hτ . The fractures are in contact, i.e., the normal displacement jump [u h ]n = 0 is zero. Going from two dimensions to three adds an additional challenge to the contact problem as we have to find not only the magnitude of the slip but also the direction. The advantage of the hybrid formulation in combination with a semismooth Newton scheme is that the same computer code can be used for any dimension, and as observed in the figure, the correct sliding direction (parallel to the Lagrange multiplier) is found by the algorithm.
The errors are calculated by comparison to a reference solution that has 500 thousand degrees of freedom. The relative errors ε Γ + i ([u h ], [u]) and ε Γ + i (λ λ λ h , λ λ λ) for the two fractures, i = 1, 2, are shown in Figure 8. We observe first-order convergence for the displacement jump, while the Lagrange multiplier shows a somewhat reduced order of convergence.

Example 3
In this example, the same setup as in Example 2 is used, but a fluid is included. The domain is sealed for the fluid, i.e., homogeneous Neumann conditions, for all sides except the top boundary, which is given a Dirichlet condition gp,D = 0 Pa. The permeability is K = 10 −8 m 2 Pa −1 s −1 , the storage coefficient c0 = 1 · 10 −10 Pa −1 , and the Biot coefficient α = 1. The initial displacement and pressure is set to zero.
Without the fractures, this setup is equivalent to a consolidation problem, which can be found in standard textbooks [38]. When the load is applied to the top surface at time t = 0, there is an instantaneous increase in the pore pressure in the domain. The fluid will then drain slowly out from the top surface and finally relax back to the initial condition. As this process occurs, the domain will continue to deform vertically increasing the mechanical load on the fractures, which causes them to slip. Twenty time steps are taken, and the simulation is stopped after 625 minutes, at which time, for practical purposes, equilibrium is reached.   The slip over time is plotted in Figure 9. Initially, the pore pressure carries most of the applied load, and the fractures are not sliding. As the fluid drains and the domain deforms, the tangential part of the Lagrange multiplier on the fractures increases, and after approximately 150 minutes, the fractures start to slide. The sliding then gradually slows down and qualitatively reaches the solution of the drained medium, i.e., the solution from Example 2. There are small differences between the solution from this example at the final time and the solution of Example 2, which are caused by the use of a dynamic friction model in this example and a static friction model in Example 2.
The number of iterations needed for convergence of the Newton solver at each time step is shown in Figure 9. For the first time step, 6 Newton iterations are needed, which is twice as many as for any of the other time steps. It is well known that the Newton strategy is very sensitive to the initial guess. A naive choice generally results in an increase in the required number of Newton iterations for smaller mesh sizes. However, either in a dynamic or a multilevel context, there are good options to set the initial guess [25,20]. In this case, the initial condition at t = 0 is (λ λ λ = 0 and u = 0), which assigns all subfaces to the noncontact set, In, while those at the first time step belong to the sticking set Iτ (see Equation (21)). During the dynamic sliding, the initial guess (the solution from the previous time step) gives a good approximation of the solution in the current time step, and thus, fewer iterations are needed. As the fractures start to slide at time step six, a few Newton iterations are needed for convergence. However, when approaching steady state, the algorithm predicts the correct slip in just one iteration.

Conclusion
In this paper, we present an approach for solving the poroelastic Biot equations in a fractured domain. A classical hybrid formulation for contact mechanics is combined with a finite volume discretization for poroelasticity. The fractures are modeled as internal contact boundaries and are governed by a nonpenetration condition in the normal direction and a Coulomb friction law in the tangential direction. The inequalities in the contact conditions are handled by a semismooth Newton method. The finite volume discretization has several advantages for these types of problems. The cell-center colocation of the discrete displacement and pressure variables gives a sparse linear system, efficient data structures, and no need for staggered grids. Moreover, the contact conditions are obtained naturally in the discretization as a condition per subface in the local systems. Thus, these conditions can be treated in an equivalent manner to boundary conditions on the external boundary. Finally, there is no need for special treatment of the contact conditions in the poroelastic case versus the purely elastic case, as the correct pressure contribution to the effective stress is obtained in the local system.
We showed that the hybrid formulation coupled with the finite volume discretization handles a given spatially varying coefficient of friction. The formulation is also suitable for other friction models such as rate and state friction or temperature-dependent coefficient of friction.
Three numerical examples illustrate the method's robustness and applicability to difficult cases. By comparison to a reference solution, the discrete solution shows first-order convergence in displacements and slightly less than first-order convergence for the Lagrange multipliers. We also show that the method handles singularity in the solution resulting from a piecewise linear fracture with a kink. Finally, a 3d example is presented where we study the effect of the fluid pressure on the solution.
The model presented in this work is limited to fluid flow in the matrix. A natural extension is to include fluid flow also in the fractures. The fluid pressure in the fracture will then act as a force on the fracture sides, effectively reducing the normal traction. Experiments have also shown that asperities along fracture surfaces can have a very important effect on both the opening and sliding of fractures. These effects can be included by adding a nonlinear deformation model to the fractures. The advantage of our framework is that any nonlinear extensions to the model can be included in the same Newton iteration, which might be crucial for the convergence of the resulting scheme. where D(x) is the distance from x to the tips of the fracture. As seen in Figure 10, this arrests the fracture before the tip, and the added regularity gives first-order convergence in both the Lagrange multiplier and displacements, as shown in Figure 11. The worst oscillations that we have encountered in 3d using our finite volume scheme coupled with the hybrid formulation are shown in Figure 12. The setup in this example is the same as the setup in Section 4.2 but with only Fracture 1 and a constant coefficient of friction, F = 0.5. Thus, we have sliding reaching the tip of the fractures. The oscillations have an amplitude of approximately 5 percent from the mean traction and grow larger as we approach the fracture tips. As in the 2d case, the displacement jump [u h ] is not effected significantly by these oscillations.
Note that the singularity at the fracture tips is a challenge for any numerical method. Similar oscillations for first-and second-order Galerkin finite elements are reported, for example, by Garipov et al [24], Fig. 8, for a setup where they study a single sliding fracture.