Algebraically stabilized Lagrange multiplier method for frictional contact mechanics with hydraulically active fractures

Accurate numerical simulation of coupled fracture/fault deformation and fluid flow is crucial to the performance and safety assessment of many subsurface systems. In this work, we consider the discretization and enforcement of contact conditions at such surfaces. The bulk rock deformation is simulated using low-order continuous finite elements, while frictional contact conditions are imposed by means of a Lagrange multiplier method. We employ a cell-centered finite-volume scheme to solve the fracture fluid mass balance equation. From a modeling perspective, a convenient choice is to use a single grid for both mechanical and flow processes, with piecewise-constant interpolation of Lagrange multipliers, i.e., contact tractions and fluid pressure. Unfortunately, this combination of displacement and multiplier variables is not uniformly inf–sup stable, and therefore requires a stabilization technique. Starting from a macroelement analysis, we develop two algebraic stabilization approaches and compare them in terms of robustness and convergence rate. The proposed approaches are validated against challenging analytical two-and three-dimensional benchmarks to demonstrate accuracy and robustness. These benchmarks include both pure contact mechanics problems and well as problems with tightly-coupled fracture flow.


Introduction
To accurately simulate the geomechanical response of a subsurface system, such as an aquifer or reservoir, it is often important to model faults and fractures [1].Phenomena such as micro-seismicity [2], fluid leakage [3], fault reactivation [4], and fracture propagation [5] depend strongly on coupled fluid-structure interaction.As a result, it is often necessary to explicitly model complex hydromechanical behavior at these surfaces [6].Experimental data confirms a strong dependence of fracture properties, like conductivity, on contact conditions [7][8][9].The core of the modeling challenge is dealing with a lubricated frictional contact problem [10].Specifically, we have fluid pressure acting as a normal forcing term on the surfaces of the fracture, while the conductivity of the fracture is a strong function of the effective aperture.This establishes a two-way and highly nonlinear coupling.A different but related challenge is to accurately model fracture propagation based on the resulting stress and strain fields near the fracture tip [11].
Many approaches have been developed over the years to solve the contact mechanics problem, and an extensive review is beyond the scope of this work.The largest difference between methods lies in how the discontinuity is fundamentally represented.A first class of models uses a Discrete Fracture Model (DFM), whereby a conformal computational grid is introduced, and the conservation equations are discretized using this grid.These methods have the advantage that standard discretization techniques may be directly applied.A potential drawback is that a conformal mesh for complex fracture networks may require distorted or excessively-refined meshes.With DFMbased models, the discontinuity can be explicitly discretized using zero-thickness interface elements [6,[12][13][14].In some models, thin layers of continuum finite elements with plastic behavior are used instead to mimic the fracture rheology [15][16][17][18].A second class employs an embedding strategy, in which a continuum discretization scheme is enriched to capture discontinuities cutting through continuum elements.For example, Embedded Discrete Fracture Models (EDFM) introduce additional, local degrees of freedom to capture opening and sliding modes [19][20][21][22] while Extended Finite Elements Methods (XFEM) introduce global degrees of freedom for this purpose [20,[23][24][25][26][27][28].A third class discards explicit discontinuity surfaces altogether, and instead represents a fracture via a regularized (smooth) field using a continuum discretization.This latter class includes Phase Field and damage-mechanics-based approaches [29][30][31][32].
In this work, we represent fractures explicitly using a conforming discretization.For these elements, the Karush-Kuhn-Tucker (KKT) conditions [33] for normal impenetrability and frictional compatibility must be enforced.The two most common classes of methods used to fulfill these conditions are penalty [13,14,34,35] and Lagrange multiplier methods [36][37][38][39][40].In penalty methods, constraints are satisfied in an approximate way using stiff "springs" connecting the two surfaces of a fracture, and no additional degree of freedom are introduced.However, penalty methods suffer from several drawbacks, including ill-conditioning of the resulting linear systems and a strong dependence of solution quality on the penalty parameters [41,42].In Lagrange multiplier-based methods, the KKT conditions are instead imposed directly.This advantage, however, comes at the cost of adding new global primary variables.Moreover, the resulting matrices have a generalized saddle-point structure that requires specialized solvers for efficiency and scalability [43][44][45].Other widely used techniques include Nitsche's method [46][47][48][49], regularized and augmented Lagrange multipliers [50][51][52][53], and mortar methods [54][55][56][57].In particular, mortar methods were originally developed to allow matching of different regions in a domain decomposition framework.They provide useful flexibility in dealing with large deformations and non-conforming discretization of the contact area [58,59].
In this paper, we focus entirely on a Lagrange multiplier approach to impose the normal and frictional compatibility conditions on the fracture surfaces.The displacement field in the bulk is approximated by lowest-order continuous finite elements.Using lubrication theory [7], the fracture flow discretization relies on a cell-centered finite volume approach using a two-point flux approximation (TPFA) scheme [60] for the numerical flux.We assume that the rock matrix is impermeable, though an extension of the proposed strategy to handle a poroelastic bulk is clear.Motivated by the coupling between fracture contact and flow processes, we use piecewise-constant equal-order interpolation of Lagrange multipliers, i.e. contact tractions, and fluid pressure.This is the natural and convenient choice given the mixed finite-element / finite-volume scheme adopted.Unfortunately, this combination of displacement and contact traction/pressure variables is not uniformly inf-sup stable [61, Section 3.1] and requires the implementation of a stabilization methodology.
To rectify this deficiency, possible remedies include enriching the displacement space with locally supported bubble functions [62][63][64][65] or using a coarser mesh for the traction/pressure space [66].Here, we instead introduce a stabilizing modification to the mass balance and constraint equations.We start by analyzing the pure contact mechanics problem, and explore several stabilization alternatives.First, we derive a local traction jump stabilization that depends on a stabilization parameter α defined locally at the macroelement level, following the macroelement analysis [67][68][69][70] originally proposed for the Stokes equation.Second, to circumvent difficulties associated with the definition of α in the presence of distorted grids and severe material heterogeneity, we develop an algebraic variant of the local traction jump stabilization technique and then generalize it to an algebraic global traction jump stabilization.Finally, we extend the stabilization procedure to problems that include fluid flow.
The paper is organized as follows.In Section 2, we present the general problem statement, providing the strong form of the governing equations.In Section 3, the discretized form and its linearization are provided.We also address the nonlinear optimization problem associated with the imposition of the KKT conditions.In Section 4, we describe three related stabilization strategies.We discuss the relative advantages of each, before settling on a global algebraic stabilization approach as our recommended strategy.We test the performance of this method on mechanical contact problems in Section 5 and hydromechanical contact problems in Section 6.We then conclude the paper with a few remarks regarding future work.

Problem statement
A fracture in a continuous medium is modeled as an internal discontinuity where additional constraints must be enforced.These constraints depend on the local state of the discontinuity-that is, whether the fracture is opening, sliding, or sticking.Inside the fracture, a single-phase fluid is present and a mass balance equation -expressed at the Darcy scale -is solved for pressure.The momentum and mass balance equations are then coupled through the fluid pressure and fracture aperture.For simplicity of exposition, we assume that the matrix is impermeable.The extension to a poroelastic medium, however, is clearly of interest in practical applications and could be included in a more general formulation [71].Similarly, extensions to multiphase flow and non-isothermal conditions may be required for some reservoir applications.
We consider an elastic, closed, polyhedral domain with Ω an open set, ∂Ω its boundary, and n Ω its outer normal vector.As usual, the boundary is subdivided into two non-overlapping subsets where ∂Ω = ∂Ω u ∪ ∂Ω σ such that ∂Ω u ∩ ∂Ω σ = ∅, where Dirichlet and Neumann boundary conditions for displacement and traction are applied (see Fig. 1(a)).From a mathematical standpoint, a fracture is described as an internal boundary Γ f embedded in Ω , consisting of two overlapping surfaces Γ − f and Γ + f as shown in Fig. 1(b).The fracture orientation is characterized by a unit vector orthogonal to the fracture plane.By convention, we choose On this lower-dimensional domain, the pressure field is defined.Let Γ f = Γ f ∪ ∂Γ f denote the closed domain occupied by the fracture, with Γ f a two-dimensional (2D) open surface and ∂Γ f a one-dimensional (1D) curve defining its boundary.The fracture boundary is subdivided into two non-overlapping subsets, ∂Γ f = ∂Γ f, p ∪ ∂Γ f,q such that ∂Γ f, p ∩∂Γ f,q = ∅, corresponding to the position of Dirichlet and Neumann boundary conditions for pressure and flux fields, respectively, as shown in Fig. 1(c).We define m f as the outer normal vector for Γ f .Finally, let T = (0, t max ) denote the time domain of interest.
We assume quasi-static conditions and infinitesimal strains.The fluid is taken to be incompressible, and we neglect body force and buoyancy effects.The strong form of initial boundary value problem (IBVP) can then be stated as [10,72,73] subject to the constraints Here, known boundary and initial conditions are given as ū : The following symbols, variables, and constitutive relationships are also introduced: • σ (u) = C : ∇ s u is the Cauchy stress tensor, which is related to the displacement vector u by the fourth-order elasticity tensor C, with ∇ s the symmetric gradient operator; µ ∇ p is the fluid volumetric flux in the fracture domain -assuming laminar flow and validity of Darcy's law [7] -with ∇ p the fluid pressure gradient, µ the fluid viscosity (constant), and C f the isotropic fracture hydraulic conductivity modeled as in [13]: Note that C f,0 captures the conductivity associated with two irregular surfaces that are in contact [74].This allows fluid to flow and pressure to propagate even if the fracture is nominally closed.From a physical viewpoint, the volume between two rough surfaces in contact is non-zero and fluid may infiltrate between asperities.For simplicity, here we assume a constant closed conductivity, though for some applications a normal-stress dependent model may be preferred.
is the traction vector over Γ f , with t N and t T = (t m 1 m 1 + t m 2 m 2 ) its normal and tangential component, respectively, with respect to the local reference system shown in Fig. 1

(b);
• τ max = c − t N tan(θ) is the limit value for the modulus of t T according to the Coulomb criterion, with c and θ the cohesion and friction angle, respectively; • * denotes the jump of a quantity across Γ f , with u = ( u| the relative displacement across Γ f , where g N and g T are normal and tangential components, respectively, and u| Γ + f and u| Γ − f are the restrictions of u on Γ + f and Γ − f .For additional details regarding the governing formulation, we refer the reader to [10,72,73]. Remark 1.In the literature, other constitutive laws are available that relate opening and fracture conductivity.For the laminar regime, they are based on improvements of the original Poiseuille law [75], describing flow between parallel smooth surfaces.Corrections are introduced to address surface irregularities.In (2), C f,0 is used to quantify the surface roughness, while other authors used a multiplicative coefficient 1/ f [7,76] or a hydraulic aperture, different from the nominal aperture g N [8].For the turbulent regime, the constitutive behavior can be based on more complex hydrodynamic relationships, like the Darcy-Weisbach equation and Moody's diagram [75].See also [27] for additional discussion on this topic.
Throughout this work, we will use subscripts N and T to identify the normal and tangential components of a vector quantity with respect to the discontinuity Γ f .Specifically, given the vector v ∈ R 3 , we have with ⊗ the dyadic product.Note that we consider the static Coulomb law.Therefore, in a discrete setting the tangential velocity ġT in (1n) can be replaced with the tangential displacement increment [61].This has to be done at every timestep if an implicit time-marching scheme is used.From now on, ġT will be replaced by ∆ n g T , where ∆ n denotes the backward difference operator such that ∆ n ( * ) = ( * ) n − ( * ) n−1 with the subscript n denoting the current discrete time level t n .The Coulomb frictional contact conditions equations (1m)-(1n) can then be rewritten as: In our framework, the fracture is explicitly modeled according to a Discrete Fracture Model [13,77], with Γ f encompassing the whole region where opening or contact may take place at any t ∈ T [73].That is, we assume Γ f is fixed and does not propagate.The only unknown is then its partitioning into stick, slip and open patches at a given point in time, . These three modes are associated with behavior regimes, namely: • Stick mode on Γ stick f : the discontinuity is fully closed and compressed with the Coulomb criterion satisfied, i.e., t N < 0 (Eq.(1j)) and ∥t T ∥ 2 < τ max (t N ) (Eq. (1m)).The three components of the traction are unknown and such that no relative movement is allowed between Γ + f and Γ − f .The conductivity is constant and equal to C f,0 .This implies linear, steady-state flow behavior within the stick portion.
• Slip mode on Γ slip f : the discontinuity is compressed, but a slip displacement g T between Γ + f and Γ − f is allowed for.Only the normal traction component t N is unknown, the tangential component having magnitude ∥t * T ∥ 2 = τ max (t N ) and direction collinear with ġT .The conductivity behavior is the same as for the stick mode, as no slip-induced dilation is modeled.The flow follows, as before, a linear, steady state model.
f and Γ − f are not in contact and a free relative displacement u is allowed for.The traction is known and equal to the zero vector in R 3 .In this case, the conductivity is related to the opening (see Eq. ( 2)) and the fluid behavior can be described as a nonlinear transient flow, with linear storage and nonlinear conductivity.
The numerical strategy to resolve this partitioning is described in the next section.

Discrete formulation
We solve numerically the IBVP (1) using a mixed finite-element/finite-volume approach for the spatial discretization combined with a fully-implicit time-marching scheme.In particular, the contact mechanics problem is addressed using the saddle point formulation [10,72,73] where traction vectors acting on Γ f are introduced as additional primary variables serving as Lagrange multipliers to enforce normal and frictional contact constraints.The simulation of the fluid flow through the discontinuity relies on a finite volume method [60].
We introduce a triangulation T of the domain consisting of nonoverlapping hexahedral cells that conforms to the discontinuity surface, i.e.Ω = ⋃ τ ∈T τ .Let us define F f as the set of quadrilateral faces ϕ in the triangulation such that Γ f = ⋃ ϕ∈F f ϕ.The edge normal unit vector to face ϕ K in the discontinuity plane is denoted by m K , with K the face global index.Let E f = (E int ∪ E p ∪ E q ) be the set of edges belonging to faces defining Γ f , with E int (respectively E p , E q ) the set of edges included in Γ f (respectively ∂Γ f, p , ∂Γ f,q ).An edge in E int shared by faces ϕ K and ϕ L is denoted as ε K ,L , with face global indices such that K < L. Similarly, an edge in E p ∪ E q belongs to a unique face ϕ K and is simply denoted as ε K .A unique orientation on Γ f is associated with each edge in E f through a unit vector m ε .We set m ε = m K for any edge ε K ,L ∈ E int and ε K ∈ E p ∪ E q .
The following discrete approximation spaces are used for the displacement and traction field [61,72] with U and U 0 the space of vector functions in [L 2 (Ω )] 3 with first derivatives in L 2 (Ω ) having trace on ∂Ω u equal to ū and 0, respectively, and M the dual space of the trace space W of U 0 restricted to Γ f .Here, C 0 (Ω ) and L 2 (Ω ) denote the space of continuous and square Lebesgue-integrable functions on Ω and Ω , Q x, y, z, x y, x z, yz, x yz} -namely, the mapping to τ of the space of trilinear polynomials on the unit cube [0, 1] 3 in R 3 -and P 0 (τ ) = [P 0 (τ )] 3 , with P 0 (τ ) the space of piecewise constant functions.Hereafter, for a given domain D, the compact notation ( * , ⋆) D is used to denote the L 2 -inner product of scalar, vector, or rank-2 tensor functions in L 2 (D), [L 2 (D)] 3 , or [L 2 (D)] 3×3 , as appropriate-for example, (∇ s η, σ ) Ω = ∫ Ω ∇ s η : σ dV .The pressure field is approximated in the space of piecewise constant functions on faces in F f , namely In this work the discretization of the mass balance equation (1b) is based on the classical two point flux approximation (TPFA) scheme.To allow for a unified presentation of the coupled mixed finite-element/finite-volume model, the TPFA scheme is written in weak form [60,78,79]. Therefore, for the space P h we introduce the following weighted inner product where is the harmonic average of positive one-sided transmissibility Υ K and Υ L associated to face ϕ K and ϕ L , respectively.In particular, the one-sided transmissibility is expressed as the product of a nonlinear (C f, β ) and a constant (Υ β ) term, respectively with |ε| the edge length, x ε a collocation point associated with each edge in E f , and x β and |ϕ β | the barycenter and area of face ϕ β .Notice that C f, β is the mean value of the fracture conductivity over ϕ β .The dependence of the hydraulic conductivity on normal displacement jump makes the transition from a closed to open state even more challenging.
Remark 2. In (7), the term qK,L = −Υ K L ( p h | ϕ L − p h | ϕ K ) represents the numerical flux approximation based on the TPFA scheme [60,80] to the volume flow rate per unit length through ε = ε K ,L ∈ E int , i.e. qK,L ≈ ∫ ε q • m ε dl.Similarly, for a boundary edge ε = ε K ∈ E p in the presence of homogeneous pressure conditions, such an approximation is given by qK = −Υ K p h | ϕ K .
Remark 3. In the derivation of transmissibility coefficients, the collocation point x ε in ( 8) is used to impose point-wise pressure continuity between faces sharing an edge .Different strategies can be used to choose x ε .Following [81], in our implementation we select x ε as the intersection of the edge ε K ,L and the line connecting the barycenters of faces ϕ K and ϕ L .For a boundary edge ε K ∈ E p ∪ E q , x ε is chosen as the orthogonal projection of x K .
Remark 4. Because of its robustness, simplicity and monotonicity-preserving properties, the classical TPFA scheme is the method-of-choice in industrial reservoir simulation.However, TPFA may lead to an inconsistent numerical flux on arbitrary grids [60,80], for which more sophisticated discretization methods like multipoint and/or nonlinear schemes should be considered, see [82][83][84] and references therein.Note that the stabilization scheme described below will work without modification for these alternative flux approximations.
Discretizing the time interval T into n T subintervals of size ∆ n t = (t n − t n−1 ), the mesh-dependent fully discrete weak form of (1) is: find with n ∈ {1, 2, . . ., n T }, and F F f (χ) and G F f (χ) two functionals that incorporate boundary conditions (1f)-(1g) Let {η i } i∈N u ∪N u be the standard vector nodal basis functions for the global finite element space of continuous piecewise-[Q 1 ] 3 functions associated with T , with N u and N u the set of indices of basis function vanishing on ∂Ω u and having support on ∂Ω u , respectively.Note that card(N u ) + card(N u ) is equal to three times the number of vertices in T .Let {χ k } k∈N p be the basis for P h , with χ k the characteristic function of the kth face in ∈ ϕ k , and N p = {1, . . ., card(F f )}.Let {µ j } j∈N t be the piecewiseconstant vector basis for M h (t N ), with N t = {1, . . ., 3 • card(F f )}-namely, the three basis functions χ j (x)n f (x j ), χ j (x)m 1 (x j ), and χ j (x)m 2 (x j ) associated to each face ϕ j ∈ F f .Discrete approximations for the displacement, traction, and pressure can then be expressed as The unknown nodal displacement components {u i,n }, face-centered traction components {t j,n }, and face-centered pressures { p k,n } at time level t n are gathered in algebraic vectors u n , t n , and p n .We emphasize that, at the righthand side of the expression for u h , the first sum represents an approximate displacement solution of the IBVP (1) satisfying homogeneous prescribed displacement, whereas the second sum is the discrete extension by zero (to the degrees of freedom) of the boundary datum ū over ∂Ω u .Hence, {η i } i∈N u is a basis for U h 0 .

Solution strategy
The discrete form of the IBVP (1) based on the weak form ( 9) consists of a nonlinear system of equations and inequalities in the unknowns u n , t n , and p n .However, if we postulate that active contact regions Γ stick f,n and Γ slip f,n are known at t n , the inequality (9b) can be replaced with a variational equality: where the four integrals correspond to the weak enforcement of the: (i) impenetrability condition on Γ stick f,n ∪Γ  11) into (9a), (12), and (9c) then allows for writing the discrete problem as a nonlinear system of equations that can be solved for the latest solution vectors u n , t n , and p n , with u n−1 the known discrete displacement solution from the previous timestep.To advance one timestep, since the partition of is of course not known a priori, an iterative procedure that includes a Newton method-based solver applied to residual equations ( 13) is needed.
In this work, we solve the contact problem using an active-set strategy, a numerical optimization technique employed in quadratic programming [85,86].Algorithm 3.1 summarizes the sequence of steps of an active-set algorithm applied to the contact problem.From a practical viewpoint, first we assign an initial status (lines 2-6) to all the Lagrange multipliers and solve the discrete nonlinear problem (13).We highlight that there are just two states, active and inactive, for both normal and frictional contact.The initial state is the previous time step solution, whenever available, or the stick state, i.e., Γ stick f,0 = Γ f , at the beginning of the simulation.Once the nonlinear problem is solved, the status of all Lagrange multipliers is checked and a new subdivision into stick, slip, and open portions is identified.If at least one multiplier changes status, we solve a new nonlinear system and check the consistency of the outcome again.The procedure stops when the starting subdivision of Γ f is consistent with the solution of the nonlinear problem.Fig. 2 provides an illustration of the resulting nonlinear convergence profile, with three consistency checks required before the final accepted solution.Remark 5. We note that there are no theoretical guarantees for the active-set convergence of the discrete version of the IBVP (1).It can occur that the exact solution may lie between two discrete solutions, i.e., two similar but different subdivisions of Γ f .In these cases, one solution must be chosen.In this work, we select the one with the smaller internal energy.
In Algorithm 3.1, Newton's method is used to drive the norm of the combined residual vector below a specified relative tolerance (line 9).To better identify each contribution associated with the contact process in the linearized problem at every iteration ℓ of the active-set algorithm, we further partition the vector containing the unknown traction degrees of freedom as , respectively.Subscripts N and T identify normal and tangential traction degrees of freedom on Γ ℓ,slip f,n .The set of indices of traction basis functions N t is also consistently expressed as union of disjoint sets such that . The solution of the nonlinear problem ( 13) is then computed as follows.Given an initial guess for displacement (u ℓ,(0) n ), traction (t ℓ,(0) β,n , β ∈ {S, N , T, O}), and pressure (p ℓ,(0) Algorithm 3.1 Active-set strategy ▷ set Lagrange multipliers to previous timestep status 4: else 5: Detailed expressions for the residual vectors and sub-matrices appearing in the linearized system (14) are given in Appendix.Note that matrices A T T and A O O are diagonal with entries equal to the area of the face the Lagrange multiplier degrees of freedom are associated with.Therefore, δt T and δt O can be eliminated through static condensation leading to the reduced system: with In the formulation proposed in [87], the Jacobian system is assembled in its reduced form with matrices B uu and B u N directly assembled as Expressions for the partial derivatives in (17) are provided in Appendix.
Remark 7. The focus of this work is on the nonlinear algorithm and stabilization only.Thus, the linear solution step in ( 14) is carried out using a direct solver.The design of a scalable solver for this linear system is clearly nontrivial, however, and is the subject of future research.

Stabilization
In this section we explore a family of stabilization techniques to enable the successful use of the spatial discretization presented in Section 3.1, We first elaborate the proposed stabilization procedures considering the simpler stick-contact problem, in the absence of frictional slip or fluid flow.The extension to more sophisticated physics will then follow seamlessly at the end of the section.
To clearly highlight the source of instability, consider a fracture entirely in stick mode, Γ f = Γ stick f .In this case, system (16) reduces to with K = A uu and C = A u S .For this saddle point system, stability and well-posedness of the discrete problem requires fulfillment of an inf-sup condition [88].In essence, "the trace space of the discrete displacement has to be well balanced with the finite-dimensional space for the surface traction" [61].Unfortunately, the Q 1 -displacement/ P 0 -Lagrange multiplier interpolation is not inf-sup stable, leading to unstable approximations [61, Section 3.1].3, without any stabilization (16 × 40 elements).In these plots, as well as in the following ones, the thick line represents the continuous solution obtained without the fault, while the points are the solution using the Q 1 -P 0 approach.
As we describe strategies to fix this deficiency, it is useful to apply them to an illustrative example for comparison purposes.We will use the 2D, plain strain problem shown in Fig. 3.The size of the domain is 8 × 20 m and the fracture has a dip of 10 • .Three material regions are considered, characterized by Young's modulus values E 1 = 3 GPa, E 2 = 15 GPa, and E 3 = (E 1 + E 2 )/2 GPa, and a homogeneous Poisson's ratio ν = 0.25.The horizontal and vertical loads are 2 MPa • m and 4 MPa • m, respectively.We compute solutions on a base grid with 16 × 40 elements, as well as on an anisotropically refined one with 16 × 200 elements.The resulting traction components orthogonal and parallel to the fracture are identified as t N and t T , respectively.We also compute a reference solution using a highly refined elastic model without a discontinuity (since pure stick conditions are assumed).In the remainder of the section, we will refer to this numerical solution as the continuous one.For each stabilized solution, we compute two integral relative differences with respect to the continuous one, E N and E T , for the two components of the traction.In Fig. 4, results using the Q 1 -P 0 interpolation without any stabilization are shown.We observe that the traction solution on average coincides with the continuous one, but it exhibits substantial checkerboard oscillations.We now consider three possible strategies to fix this issue: 1. Analytic macroelement stabilization 2. Algebraic macroelement stabilization 3. Algebraic global stabilization As we will see, the first method makes significant assumptions regarding the grid topology and material heterogeneity, while the latter methods provide greater generality.

Analytic macroelement stabilization
We first explore stabilizing the discretization using a macroelement approach [89].The method is most easily described using a 2D reference macroelement as shown in Fig. 5(a).This macroelement is formed by two interface elements and four quadrilaterals, two for each side of the fracture.To create a patch test, the ten displacement degrees of freedom (DOFs) located on the boundary are fixed.The objective is to derive a stabilization that produces a well-posed problem on an individual patch.Such a stabilization then implies well-posedness on a grid consisting of stabilized macroelements.
We can achieve this by investigating the kernel modes of the Schur complement for the Lagrange multipliers.For example, for the macroelement of Fig. 5(a), there are 4 displacement DOFs -the xand y-components for innermost nodes, above and below the fracture -and 4 traction DOFs-one normal and one tangential component for each element on Γ f .Assuming the ordering u = [ u ] for the unknowns, the explicit expressions of K and C are: with E the Young's modulus, ν the Poisson's ratio, h x and h y the mesh size in the xand y-direction, respectively, and R 2 a 2 × 2 rotation matrix from the local to global reference system.In the example of Fig. 5(a), R 2 reads with {e x , e y } the standard Euclidean basis in R 2 .By definition, the Schur complement S ∈ R 4×4 is Performing an eigen-decomposition, its complete set of eigenpairs We observe that S has rank 2. Eigenvectors v 1 and v 2 span the kernel of S and represent the checkerboard mode for the traction normal-and tangential-component, respectively, confirming the numerical results shown in Fig. 4. Conversely, the component-wise constant eigenvectors v 3 and v 4 span the column space of S.
The null eigenvectors v 1 and v 2 are the source of the macroelement instability and need to be removed from the null space of S. To do so, we introduce the symmetric and positive semi-definite matrix H * = V V T , with V the matrix having columns v 1 and v 2 .By construction, v 1 and v 2 are eigenvectors of H * , both corresponding to the eigenvalue 2, which span the range of H * .Also, v 3 and v 4 are now a basis for the kernel of H * .Let H = α H * be a scaled matrix, with α a scalar stabilization constant having units of squared length per pressure, i.e. the same units as the eigenvalues of S. The system to be solved is modified as where a stabilizing contribution now replaces the zero block of the original matrix.The resulting modified Schur complement is then S = −C T K −1 C −α H * .The eigenpairs of S are the same as those of S with the only difference that v 1 and v 2 are now associated with the same non zero eigenvalue 2α.The scaling constant α is chosen in such a way that, in a homogeneous case, the eigenvalues of S are bounded between min(λ 3 , λ 4 ) and max(λ 3 , λ 4 ), so that the matrix conditioning does not depend on the stabilization constant.For a regular Cartesian grid with h x = h y = h, from Eq. (22b) we can write thus, the "optimal" scaling factor simply reads: In case of geometric anisotropy, h can be computed as the average length of the interface elements composing the macroelement.
We emphasize that the stabilization matrix H represents an example of a minimal stabilization operator [90] that does not pollute the physical eigenpairs.It requires, however, explicit knowledge of the eigenvectors associated with non-zero eigenvalues.For the macroelement of Fig. 5(a), H can also be interpreted as a macroelement stabilization matrix for the α-weighted inter-element traction (component-wise) jump since the following relationship holds true Applying this stabilization to the 2D test case (Fig. 3) yields a substantial reduction of the oscillations observed in the unstabilized case as shown in Fig. 6, with the errors E N = 1.32•10 −1 , E T = 1.36•10 −1 and E N = 7.36•10 −2 , E T = 1.55 • 10 −1 for the base and the refined grids, respectively.Remark 8.The 2D reference macroelement (Fig. 5(a)) consists of four equal rectangular elements -hence, the stiffness matrix K , see Eq. ( 19), is diagonal -with the normal n f aligned with the y-axis.Consequently, normaland tangential-component of the traction are decoupled as the eigenvectors defined in Eqs.(22a)-(22b) reveal.This is not the case in a more general geometry where the two traction components are typically related.
Extension to 3D is based on the reference macroelement shown in Fig. 5(b), which consists of four quadrilateral interface elements and eight hexahedral elements.The derivation follows the same steps discussed above for the 2D case.We omit the computations and simply report the main results.The Schur complement S ∈ R 12×12 has rank three with column space spanned by three component-wise constant eigenvectors.There are nine spurious traction modes that need to be stabilized.A convenient basis for the kernel of S is shown in Fig. 7 and corresponds to checkerboard-like modes for each component of the traction vector with respect to three internal edge of the interface element patch.Note that the choice of the three edges is arbitrary.The optimal α value to obtain a stabilization contribution that falls in the already present lower/upper spectral Schur complement bounds is: As for the 2D case, this value is exact in case of h x = h x = h z = h and constant physical properties.In general, h can be computed as the cubic root of the average volume of the elements surrounding the fracture and sharing a face.
Fig. 7. Three different modes that need to the stabilized for the 3D case.Each of them applies to the three components (one normal and two tangential) of the traction vector, thus the kernel size is 9.Note that the remaining mode is a linear combination of these three, thus it is part of the same space, the Schur complement kernel.

Algebraic macroelement stabilization
The scaling factor α in the macroelement approach above is a scalar value collecting both mechanical and geometric information.Its definition can be quite difficult, especially for 3D problems with distorted grids and material heterogeneity.In this section, we propose an algebraic alternative for computing the stabilization matrix at the macroelement level that circumvents the need for introducing the factor α. The development is based on the following observation: α is fundamentally needed to scale matrix H * introduced in Section 4.1 so that the spectrum of the stabilized Schur complement S is bounded between the smallest and largest nonzero eigenvalues of S. Thus, whenever we provide a different but admissible scaling, α can be avoided.
Let us consider a general patch of interface elements in a non regular macroelement (Fig. 8).The only topological assumption we make is that the fracture elements within an individual macroelement are co-linear in 2D or co-planar in 3D.For the two dimensional case of Fig. 8(a), matrix C given in Eq. ( 19) reads where l1 and l2 are the interface element lengths associated to the vertex in common between elements ϕ 1 and ϕ 2 , i.e. the integral of the standard hat function associated to that vertex over ϕ 1 and ϕ 2 , respectively.The rows of the following matrix C, formed from C swapping block columns and changing sign to the first block column, are orthogonal by construction to C, i.e., C C T = 0. Thus, columns of C T belong to the kernel of the Schur complement S = −C T K −1 C. Also, C has rank 2, the number of modes that are known to require the stabilization.Hence, C T C is a stabilizing contribution to the Schur complement.It has to be scaled, but from the observation that C has the same entries of C, except for the order (they are swapped on a Lagrange multiplier base) and the sign, it Fig. 9. Comparison between the analytic (continuous line) and the algebraic (dashed line) macroelement stabilizations.On the y axis is the condition number κ( S) = λ max ( S)/λ min ( S) of the stabilized Schur complement.Panels from left to the right show (i) regular grid and homogeneous material properties, (ii) stretched grid with homogeneous material properties, and (iii) regular grid with heterogeneity in the material properties.
is natural to use the inverse of the stiffness matrix in (19), i.e.K −1 , to scale it.Indeed, a scaling that incorporates element size and material properties information is needed and K collects both of them.Numerical tests show that the inverse of its diagonal, denoted as D K ∈ R 4×4 , is enough.Therefore, a macroelement stabilization matrix can be defined as H = C T D −1 K C. Using the same notation of Section 4.1, such a stabilization can also be expressed in the following compact form: with D ∈ R 2×2 a diagonal matrix given by the sum of the inverse of the diagonal portion of the stiffness matrix associated with the two nodes (on both the discontinuity surfaces) shared by the two interface elements forming the macroelement (i.e., N + 1 and N − 1 in Fig. 5(a)).Note that whenever applied to a regular grid, with h x = h y , and homogeneous material properties, this stabilization reduces to the α-based method of the previous section.
To compare the α-based and algebraic macroelement stabilizations, in Fig. 9 we report the behavior of the condition number κ( S) = λ max ( S)/λ min ( S) of the stabilized Schur complement for different macroelement configurations.The x-axis is normalized with respect to the "optimal" value according to Eq. ( 25), the dot is the conditioning arising from the α-based stabilization and the dashed line is the conditioning of the algebraic stabilization.From the left to the right, we have three cases: • Regular grid, with homogeneous material properties.It can be observed that the two techniques offer the same optimal result.
• Stretched grid (Fig. 3) with h y,τ 1 = 5h y,τ 3 .The material is still homogeneous.In this case, the stabilization based on α is not able to get the minimal conditioning, while the algebraic one is.
• Regular grid (Fig. 5(a)), with heterogeneity in the Young's modulus: Again, the stabilization that uses a scalar value is not able to provide the minimal conditioning, while the algebraic one is.
Observing the plots in Fig. 9, it can be noticed that the conditioning has a unique minimum when the grid is regular, but there is a set of minima in the middle case, with a stretched grid.The former occurrence implies that the two eigenvalues λ 3 and λ 4 are the same, while in the latter that they are different.This finding is consistent with the analytical definition of the Schur complement eigenvalues (see Eqs. (22b)), where the aspect ratio h x / h y is a parameter.Using a more refined definition of α * , considering also the geometric anisotropy, would improve the conditioning number of the α-based S for this specific case.In general, however, as assumptions regarding mesh geometry and material properties are relaxed, an algebraic approach is increasingly appealing in its simplicity.
To further test this approach, we focus again on the model problem depicted in Fig. 3. Numerical results are provided in Fig. 10, with the errors E N = 1.33 • 10 −1 , E T = 1.35 • 10 −1 and E N = 7.49 • 10 −2 , E T = 1.63 • 10 −1 for the base and the refined grids, respectively.Comparing this technique with the α-based approach, we observe quite similar performance.We conclude this section by discussing the extension to 3D problems.Matrix C for a 3D macroelement reads where the areas Ãi , i = {1, . . ., 4}, are the counterpart of lengths l1 and l2 (Fig. 8(b)), and R 3 is a 3D rotation matrix.Using the same arguments as in the 2D case, matrix H can be obtained by stabilizing the traction jumps across internal edges of the macroelement as Remark 9.No explicit assembly process is needed in the construction of V and D both in 2D and 3D.A local gathering from the global stiffness and coupling matrix is enough to extract the required blocks.

Algebraic global stabilization
In this section, we extend the algebraic approach of Section 4.2 to overcome the key limitation of any macroelement-based approach, the topological restriction that the mesh be partitioned into macroelements in the first place.As previously stated, the key point is to stabilize the jump between two interface elements.Rather than doing this for macroelement-internal edge alone, we now simply stabilize all possible jumps between each pair of interface elements.
form C by swapping columns of C loc and changing sign 5: compute stabilization contribution To do so, it is sufficient to compute the matrix Cloc for all degrees-of-freedom shared between two elements discretizing the fracture, gather the local matrix K loc , and then assemble individual H loc contributions into a global stabilization matrix H . Algorithm 4.1 describes the approach for the 2D case -see Fig. 11 -using a Matlab-style notation.Note that the proposed algorithm assumes, implicitly, that fractures are planar.Also recall that multiple DOFs are associated to each geometric object (nodes and faces) and thus the gather/scatter indexing refers to vector sets of component indices.
We note that the resulting stabilization matrix has the same stencil as a two-point flux approximation, having been assembled using points (edges) shared between pairs of elements in 2D (3D).For a quadrilateral (hexahedral) mesh, it is a 3-point (5-point) block Laplacian in 2D (3D).By block Laplacian, we mean that each matrix super-entry is composed of a 2 × 2 or 3 × 3 dense block, according to the space dimension, with one row/column for each component of the traction vector.
This Laplacian stencil can be easily built using the pre-assembled global C T matrix, because its pattern provides the DOF-connectivity between faces and nodes for the fracture discretization.In 3D, we must take care to avoid the assembly of elements sharing just one node, instead of an edge.We can do so by computing the product C T 1 C 1 , where C 1 is a matrix with the sparsity pattern of C but all non-zero entries set to one.Each row/column entry of the product matrix is then the number of shared nodes between the respective elements.We may then readily filter out any element pairs that share just one node.
As with other approaches, Fig. 12 shows the performance on this method on our comparison problem.The resulting errors are E N = 7.49 • 10 −2 , E T = 1.46 • 10 −1 and E N = 5.27 • 10 −2 , E T = 1.16 • 10 −1 for the base and the refined grids, respectively.Comparing this technique with the other two approaches, this method produces smaller errors.To better highlight the different behavior of the three methods, we solve again the model problem of Fig. 3 using a more refined grid with 64 × 160 elements.We report results in terms of normal component t N only, being representative of overall performance.In Fig. 3, the black circle identifies the portion of the fracture represented in Fig. 13.The global relative error E N , shown therein, is computed as before.
Given its straightforward computation and good performance, we adopt the algebraic global stabilization as our default approach.
Remark 10.We observe that, for all three approaches, the stabilization matrix H is symmetric and positive semi-definite (SPSD).We recall that, being the Schur complement of a negative definite matrix, it is added with the minus sign, i.e., S = S − H .  Remark 11.The idea of using a global stabilization based on jump penalization is not completely new, having been suggested for other discretization and interpolation schemes [92][93][94].

Inclusion of friction and fluid flow
The methods described for the stick problem form the building blocks for stabilizing the full physical model, allowing for frictional sliding and fluid flow.In particular, all that is required is to modify the original system ( 14) to its stabilized version, where five stabilizing H i have added to particular components of the multiphysics problem.In particular, the traction components t S and t N and the pressure field p must all receive stabilizing contributions.Stability is primarily provided by the matrices H SS , H N N , and H P P , which are assembled on macroelements using the same procedures as described in the previous section, but now applied to different traction and pressure components as warranted.The off-diagonal terms H S N and H N S are necessary to capture contributions from stabilized edges between elements in a "mixed" state -that is, where an element in a stick state is adjacent to one in frictional sliding state.We observe that t T -i.e., the friction component of the traction on Γ slip f -is a function of the relative displacement rate ∆g T and therefore does not need to be stabilized.Similarly, for Γ open t , there is no need for stabilization, as the traction t O is a priori known to be zero.As before, because the matrices A T T and A O O are diagonal, a Schur-complement reduction can be performed to eliminate t T and t O .The reduced system will closely resemble Eq. ( 16), but with the stabilizing blocks above included in the appropriate places.

Analytical benchmarks for contact mechanics without fluid flow
In this section, we validate the formulation, the discretization, and the stabilization strategies for pure contact mechanics problems, without fluid flow.

Analytical benchmark for shear behavior: Constant solution
This analytical benchmark was originally proposed in [91].The representation of the 2D model domain is reported in Fig. 14(a), where the lower boundary is y-constrained, the circled corner is constrained also in the x direction and on the upper boundary a uniform displacement is imposed ( ūy = 0.10 m).The material is homogeneous with elastic parameters E = 5000 MPa and ν = 0.25.Coulomb's frictional parameters are c = 0 and θ = 5.71 • , such that the friction coefficient is 0.1.The solution is a constant sliding on the fracture of value  We highlight that the simulation is carried out on a 3D mesh (Fig. 14(b)), because in a 2D setting the direction of the shear vector is known; thus, the Coulomb frictional contact condition is not needed and the shear behavior is linear.To test the nonlinearity introduced by the dependency of the direction on the relative displacement rate, when the failure condition is matched and the fracture slides, we need to work in a 3D setting.To match a 2D analytical solution, we respect the plane strain assumption and build the domain extruding the surface shown in Fig. 14(a) by 0.5 m and z-constraining the two surfaces parallel to the x-y plane.
From Figs. 14(c)-15(a), it is clear that the model is able to match the analytic linear solution even with a small number of elements (4680 nodes, 3610 hexahedra and 95 quadrilaterals for the fracture).We report also the convergence profile of the nonlinear solution algorithm 3.1 in Fig. 15(b).In this specific case, after the elastic step, i.e., the first solution with Γ f = Γ stick f , the final configuration of Γ f is obtained and no other outer loop iterations are required.

Single crack under compression
The second example is a single crack in a 2D infinite domain under a constant uniaxial compression.The benchmark geometry is described in detail by [95] and reproduced in Fig. 16(a).Being E and ν the linear elastic parameters, σ the compressive stress, ψ the fracture inclination, 2b its length and θ the friction angle for Coulomb's criterion (with zero cohesion), the analytical solution provides the normal traction on the fracture and the sliding on it: respectively, where 0 ≤ ξ ≤ 2b is a local coordinate on the fracture.A plane strain status is assumed.For the simulation, we used the following values: E = 25000 MPa and ν = 0.25, the friction angle is θ = 30 • , the fracture is tilted by ψ = 20 • and extends for 2b = 2 m, and |σ : (e x ⊗ e x )| = 100 MPa.
The model is discretized with different resolutions in the x-y plane, from 5K to 68K quadrilaterals and from 108 to 432 interface elements, for the coarsest and the finest grid, respectively.The final 3D domain is obtained by extrusion.The less refined computational domain is shown in Fig. 16(b), where a zoom on the mesh around the fracture is also provided.The boundary conditions for u x and u y are set in order to respect the symmetry of the expected solution (see Fig. 16(b)).The two faces parallel to the x-y plane are constrained in the z direction, because of the plane strain assumption.
Comparisons between numerical and analytical solution are provided for the frictional traction (Fig. 17(a)) and the relative displacement (Fig. 17(b)).It can be noticed that the computed displacement is in good agreement with the expected one everywhere, while the traction is quite different very close to the fracture tip, where some oscillations are observed.Refining the mesh, these oscillations are still present, but always closer to the fracture tip.From a physical viewpoint, these can be explained as a locking phenomenon [91].Given the imposed external load, the fracture tries to slip, but the elements close to the tip are not allow for because of the tip itself.To accommodate the expected behavior, the fracture has to open near the tip, even if there is a compressive normal traction on the crack.With a less refined mesh, the interface element with a fixed edge on the tip opens, yielding a positive normal traction.In agreement with these considerations, the rightmost and leftmost traction values are much lower than the average, tending to zero, i.e., the near-tip elements tend to open.
Using the previous example, where an analytical solution for displacement and traction field on the fracture is available, we study the convergence rate, i.e., the error dependence on the mesh size.First of all, we observe that the oscillations on t N (see Fig. 17(a)) will prevent any traction error norm from convergence, thus, as done in [96], to compute a meaningful norm we neglect the extreme portions of the domain, i.e., the traction related norms are computed on the central 90% of the fracture trace.In Fig. 18, we show the convergence of the two error norms, on the normal Lagrange multiplier (Fig. 18(a)) and on the sliding component of the displacement (Fig. 18(b)).While the first one is a properly defined norm on the domain Γ f , where t is defined, the latter one is just the norm of a continuous field projected on a surface, indeed u is defined on Ω but the norm is computed on Γ f only.The order of convergence m is slightly larger than 1 for the Lagrange multiplier norm, while it is close to 1.4 for the projection of the displacement.There is no well-defined value for the convergence rate of the mixed finite elements space used in this work.The rate 1.11 obtained here is in agreement with the theoretical results given in [61, Section 4].

Zipper crack problem
To complete the validation of the model, we use the line crack problem as described in [97].The zipper crack case is also known as Griffith problem [98] and consists of an infinite plane (x-y plane) with a single linear crack of length 2l in the x-direction.Inside the fracture, there is a fluid with a given pressure p(x).The assumption is of plane strain.Fig. 19(a) represents a sketch of the setup.The analytical solution provides the opening of the fracture  for every location x and the stress in the continuous medium on the x direction (y = 0), for x > l.The far-field stress orthogonal to the fracture is σ 0 .We impose a constant pressure p 0 only on one part of the fracture, in such a way that the tip closes smoothly.The pressurized length ends at x 0 , with: Introducing ⏐ ⏐ , the analytical solutions for opening and stress in y direction are [97]: In Fig.In Figs.20(a) and 20(b), we report the comparison between the outcomes of the numerical model and the analytical solutions.It can be notice, for x > l, the fracture is closed and the transition is smooth, as predicted by the theoretical solution.The stress behavior is very similar to the analytical one, except very close to the tip, where the numerical model shows a smaller jump between the fracture and the continuous material.Nevertheless, we have a smooth transition between the solutions computed with two different discretizations-on a portion of the closed length interface elements are used, while on the other classical finite elements are employed.Overall, we can see a good agreement, with an integral relative error of 1.0% for the displacements and 2.5% for the stress.For the error evaluation related to the stress, we neglected the near-tip portion of the domain, i.e., we considered only x > 5/4l = 12.5 m.
As a concluding remark for this analytical benchmark, Fig. 20(c) shows the convergence rate for the L 2 -norm of the error on the fracture aperture g N (x) only, as the stress field is unbounded close to the tip, so unsuited for such a test.As expected, the order of convergence is around 1 because of the discontinuity in the stress field, indeed, the asymptotic behavior is almost always lost whenever there are singularities in the solution [99].The coarser mesh has 1K finite elements and 60 interface finite elements, while the finest has 194K and 900 finite and interface finite elements, respectively.

KGD problem
We consider a 2D hydraulic fracture propagation assuming plane strain conditions.The medium is isotropic, homogeneous, impermeable, and is fully described using a linear elastic model.An incompressible Newtonian fluid is injected from a fixed point, at a constant rate Q 0 .The fracture propagates in the direction that is orthogonal to the maximum principal direction of the stress tensor in the surrounding medium.In Fig. 21(a), we represent the set-up of the problem and introduce the quantities of interest: • g N (x, t): the fracture opening for any time and position x ≤ l; • p(x, t): the net fluid pressure inside the fracture for any time and position x ≤ l; • l(t): the fracture half-length.The analytical solution is provided in terms of g N (x, t), p(x, t) and l(t), but the complete expressions require the definition of some dimensionless quantities [103]: (i) an opening Ω , (ii) a net pressure Π and (iii) a fracture length γ .Given these dimensionless variables, the quantities of interest become: where ξ = x/l is the similarity variable, i.e., a dimensionless fracture coordinate.For the zero toughness case, we have [102]: In (39), we introduced Ω = Ω γ , µ ′ = 12µ l , with µ l the fluid viscosity, and E ′ = E 1−ν 2 , i.e., the plane strain modulus.The two functions Ω (ξ ) and Π (ξ ), called self-similar fracture opening and self-similar net fluid pressure, respectively, are approximated through polynomial series, based on the Gegenbauer polynomials [108] for the former one, while the latter uses Euler's beta functions and Gauss's hyper-geometric functions [108].The full expression of these two dimensionless functions, with the numerical coefficients of the series expansion, is provided in [102].This analytical solution is based on zero-toughness and zero-lag assumptions.The first hypothesis implies that the energy dissipated by the fracture propagation is negligible compared to the energy dissipated in the fluid by viscous flow.In our current implementation, the fracture surface is pre-defined and no energy is dissipated in fracturing.Regarding the second assumption, for real fractures a fluid lag or non-wetted zone may appear.From physical considerations, the fluid pressure at the tip has to be finite, as well as the stress in the rock surrounding the fracture tip.In the regime of interest, however, the influence of this region on the global pressure and displacement solution is confined in a small region near the tip (see numerical results [109,110]).Given these considerations, the analytical solution used herein is a good approximation of the main physical behavior.We emphasize the fact that lim ξ →1 Π = −∞, thus, the analytical solution predicts an infinite pressure at the fracture tip.This nonphysical result is due to the assumption that the fluid reaches the tip and fill every empty space at the same density, without allowing for cavitation.Note that both the fluid density and the initial stress regime do not affect any quantity of interest.
In Fig. 21(b), we show the mesh used for the simulation.It is composed of 10209 nodes, 6600 hexahedra and 76 interface elements.The domain is 300 × 600 × 8 m, a fracture of 150 m and an average element size of h x ≈ 4 m.We highlight that currently our model does not handle fracture propagation, but we know a priori the fracture trajectory and can pre-discretize a surface of sufficient length.In some sense this is then a fracture "reactivation" problem.To minimize the boundary effects, we imposed the symmetric condition on the face parallel to the y axis and intersecting the fracture and z-constrained the two surfaces parallel to the x-y plane, to reflect the plane strain assumption.The 2D node highlighted in Fig. 21(b), that corresponds to a column of nodes in 3D, is the only one that is y-constrained.The material parameters are E = 30 GPa and ν = 0.25.The fracture frictional behavior is governed by Coulomb's criterion, characterized by θ = 30 • and zero cohesion.The fluid viscosity is µ l = 10 −9 MPa • s.The injection rate and the confining stress are Q 0 = 2 • 10 −3 (m 3 /s)/m and σ 0 = 10 MPa, respectively.Finally, according to [74], the conductivity initial value (see Eq. • m.The simulated time is 100 s, with ∆t = 0.01 s from the beginning to t = 0.2 s, then ∆t = 0.1 s up to t = 2 s and finally ∆t = 1 s up to the end of the simulation.We highlight that for the current flow rate, the Reynolds number is about 2000, thus the regime is laminar and the cubic law assumption is reasonable. Figs. 22(a) and 22(b) show the outcomes of the model in terms of opening and pressure for two different time steps, i.e., at half (t = 50 s) and at the end of the simulation (t = 100 s).The continuous line is the analytical  solution.Overall, there is a good agreement, for both the aperture (with an integral relative error of 4.4% for t = 100 s) and the pressure (with an integral relative error of 3.8% at the same time step).At time t = 100 s, the pressure is slightly different close to the tip, where the theoretical behavior tends to −∞.We emphasize that the analytical solution for the pressure is constrained in an integral sense by the propagation criterion, being In our model, we use a Dirichlet boundary condition on the pressure value at the fracture end, that is not the fracture "tip", where p = 0 is imposed.To compare the two behaviors, our outcome in terms of pressure is shifted by a constant value.In Fig. 22(c), we report the fracture length as function of time.The model is able to predict quite accurately the fracture length, with an average relative error of 2.1%.Finally, in Fig. 22(d), the pressure at the injection location is shown for all time steps.Because of the different strategies used to impose the pressure boundary condition, we cannot directly compare this profile with the analytical solution.We limit ourselves to a visual comparison with similar works, e.g.[14,27,111], observing profiles that are in excellent agreement.

Penny-shaped crack
The problem consists in an axisymmetric hydraulic fracture in an infinite medium, that is isotropic, homogeneous, impermeable and with a linear elastic behavior.From the center of the fracture, an incompressible Newtonian fluid is injected, at a constant rate Q 0 .The fracture propagation does not depend on the far-field stress status and, as in the KGD example, it is enough to solve for the net fluid pressure.In Fig. 23(a), we represent the setup for the problem.The quantities of interest are the same as before, except the fracture half-length, that is now substitute by R(t), i.e., the fracture radius.
The behavior in space and time, represented by g N (x, t), p(x, t) and L(t), is provided by the analytical solution through some dimensionless quantities [103]: (i) an opening Ω , (ii) a net pressure Π and (iii) a fracture radius γ .
The main quantities can be expressed from these dimensionless functions as: where the similarity variable ρ = r/R is the dimensionless fracture coordinate.For the zero toughness case, we have [106]: In (42), we introduced Ω = Ω γ , µ ′ = 12µ l and E ′ = E 1−ν 2 as in the KDG study.The self-similar fracture opening Ω (ρ) and the self-similar net fluid pressure Π (ρ) are expressed as sum of a general and a particular solution.The first respects the governing equations, while the second represents the inlet asymptotic behavior.According to [106], the chosen basis function is a combination of Jacobi polynomials of arbitrary order in the interval [0, 1].In this reference, there is the complete expression for the two self-similarity functions.In the current work, we use just the first order expansion for both Ω (ρ) and Π (ρ).The pressure solution is unbounded, being lim ρ→0 Π = +∞ and lim ρ→1 Π = −∞, and the model can struggle in the approximation of this nonphysical values happening at the extremes of the domain.As for the KGD solution, we highlight that the analytical solutions are not affected by neither the fluid density nor the initial stress regime.
The computational domain simulates one fourth of the problem, composed of 16061 nodes, 13896 hexahedra and 640 interface elements, as shown in Figs.•m.As for the previous case, the simulated time is 100 s, with ∆t = 0.01 s from the beginning to t = 0.2 s, then ∆t = 0.1 s up to t = 2 s and finally ∆t = 1 s up to the end of the simulation.According to [112], the effective Reynolds number is 1.1 • 10 −1 after 1 s and 3.2 • 10 −3 at the end of the simulation.Thus, the laminar flow assumption required by the cubic law holds.We emphasize that the proposed mesh is not really suited for a TPFA finite volumes solution scheme, nevertheless, the results prove to be accurate enough to verify the analytical solution.
For two different simulated times, i.e., at half (t = 50 s) and at the end of the simulation (t = 100 s), we show the opening and pressure profiles for the y = 0 fracture trace in Figs.24(a) and 24(b), where the continuous line is the analytical solution.Overall, there is a quite good agreement, for both the aperture (with an integral relative error of 5.2% for t = 100 s) and the pressure (with an integral relative error of 4.2% at the same time step).The pressure is slightly different close to the extremes of the domain, i.e., ρ = 0 and ρ = 1, because the theoretical behavior diverges, tending to ±∞.As in the KGD case, the propagation criterion provides an integral constraint for the pressure analytical solution, being Using a simple Dirichlet boundary condition on the pressure value at the fracture maximum radius, where p = 0 is imposed, we need to shift our outcome for the sake of comparison.Fig. 24(c) represents the fracture radius as function of time.The model is able to predict quite accurately the fracture length, with an average relative error of 2.2%.Finally, in Fig. 24(d), the pressure behavior for the entire simulation at the injection point is shown.As for the KGD example, we rely on a visual comparison with other solutions available in the literature, see for instance [14].

Conclusions
In this work, we have presented a stabilized displacement-Lagrange multiplier-pressure formulation for quasistatic contact mechanics coupled with fracture fluid flow.Our discretization is based on a finite element method for the contact mechanics subproblem combined with a finite volume scheme for the flow subproblem.The global nonlinear problem is solved using an active set strategy.We employ lowest-order continuous finite elements for the displacement field, piecewise constant functions for both Lagrange multiplier (traction) and pressure field, and a linear two-point flux approximation for intercell numerical fluxes on the discontinuity surface.
The mixed space we adopt does not automatically satisfy the discrete inf-sup stability condition.To stabilize the formulation, we began by revisiting the well-known macroelement approach, originally proposed for the Stokes equation.We then developed two new algebraic strategies that work under less restrictive assumptions while providing good performance.The resulting approach, both without and with fluid flow in the fracture, has been benchmarked against analytical solutions to validate its behavior in case of normal and frictional activation, as well as in two-and three-dimensional fracture reactivation.
Future developments will deal with the simulation of fluid flow in the matrix coupled with the structural mechanics problem through porosity variation and fracture flow through leak-off.Networks of fractures, for which the computation of the hydraulic conductivity requires particular care, will also be investigated.Finally, a key step to enable high-performance computing applications is the design of an efficient preconditioning strategy for the particular Jacobian linear systems encountered here.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.with the partial derivatives expanded as: η j,N dA (A.3e)

Fig. 1 .
Fig. 1.Conceptual scheme for fracture modeling (a); local reference frame on the discontinuity surface (b); and the scheme for the flow domain (c).
slip f,n , (ii) no slip condition on Γ stick f,n , (iii) known tangential traction on Γ slip f,n given by Eq. (1n), and (iv) zero traction condition on Γ open f,n .Introducing expressions (

Fig. 2 .
Fig. 2. Example of convergence profile using the active set strategy.It is clear the jump in the residual after each consistency check.

Fig. 3 .
Fig.3.Sketch of the 2D test case used to compare stabilizations.The black circle locates the area where a more detailed zoom is taken.

Fig. 4 .
Fig. 4. Numerical results for example in Fig.3, without any stabilization (16 × 40 elements).In these plots, as well as in the following ones, the thick line represents the continuous solution obtained without the fault, while the points are the solution using the Q 1 -P 0 approach.

Fig. 5 .
Fig. 5. Location of displacement and traction degrees of freedom (DOFs) for a 2D (a) and a 3D (b) macroelement, respectively.Diamonds ( ) show the location of displacement DOFs, while squares ( ) represent the position of traction DOFs.Black ( ) and gray ( ) mean fixed and active DOFs, respectively.

Fig. 6 .
Fig. 6.Numerical results for example in Fig. 3 using the analytic macroelement stabilization.First row: base grid, last row: vertically refined grid.

Fig. 10 .
Fig. 10.Numerical results for example in Fig. 3 using the algebraic macroelement stabilization.First row: base grid, last row: vertically refined grid.

Fig. 11 .
Fig. 11.Sketch of two elements sharing an edge.Common nodes on both surfaces of the fracture are highlighted.

Fig. 12 .
Fig. 12. Numerical results for example in Fig. 3 using the algebraic global stabilization.First row: base grid, last row: vertically refined grid.

Fig. 13 .
Fig. 13.Comparison of different stabilization techniques for a refined regular grid.From the left to the right: analytic macroelement, algebraic macroelement, and algebraic global stabilization approach.

Fig. 16 .
Fig. 16.Single crack under compression in infinite domain: (a) sketch of setup; (b) computational mesh; (c) zoom close to the crack.Diamonds ( ) and squares ( ) indicate u x = 0 and u y = 0 Dirichlet boundary conditions, respectively.

Fig. 17 .
Fig.17.Single crack under compression in an infinite domain: comparison between analytical and numerical solutions.

Fig. 18 .
Fig. 18.Single crack under compression in an infinite domain.Convergence rates for: (a) normal traction component, and (b) tangential component of the relative displacement across the crack.The order of convergence m is computed through a least-square linear interpolation.

Fig. 19 .
Fig. 19.Zipper crack problem: (a) sketch of the setup; (b) computational domain with the fracture highlighted as a gray segment.

2 √ 3 ≈
19(b), we show the computational domain used to reproduce this solution.The 3D domain, whose size is 150×300×0.3m, is discretized with 17772 nodes, 13050 hexahedra and 228 interface finite elements.The objective is to simulate a fracture length of l = 10 m, but we discretize a longer fracture, l 1 = 15 m, to verify if the remaining length l 1 − l = 5 m remains closed.We need a large domain to attenuate the boundary effect as the analytical solution is obtained for an infinite domain.It is 10 × 20 times the actual fracture length.Referring to Fig.19(b), the only constrained boundary is the one intersecting the fracture, where symmetric conditions are imposed.The only fixed point is the one highlighted in the figure.The 3D domain is obtained extruding a 2D domain, and the faces parallel to the x-y plane are z-constrained, to fulfill the conditions of the assumed plane strain state.Regarding material properties, we have E = 25 GPa and ν = 0.25.Fluid pressure and far-field stress are p 0 = 15 MPa and σ 0 = 10 MPa.With the chosen set of parameters, we have lim x→l σ y (x, 0) = −15 + 30 π arctan −6.816 MPa, thus there is a discontinuity in the stress field at the fracture tip.

Fig. 20 .
Fig. 20.Zipper crack problem: (a) comparison between numerical and analytical solution for fracture opening; (b) comparison between numerical and analytical solution for σ y -stress in the continuous domain; (c) convergence rate for fracture opening with the order of convergence m computed through a least-square linear interpolation.In (b), the gray diamond ( ) indicates where the discretization with interface elements ends and classical finite elements are used.

Fig. 21 .
Fig. 21.KGD problem: (a) sketch of the setup; (b) computational domain with the fracture highlighted as a gray segment.The square ( ) represents the projection on the x-y plane of the y-constrained nodes.

Fig. 22 .
Fig. 22. Results for the KGD simulation.From the left to the right: fracture opening, fluid pressure, fracture length and pressure behavior at the injection location.

Fig. 23 .
Fig. 23.Penny-shaped crack problem: (a) sketch of the setup; (b) x-y view of computational domain; (c) and full 3D mesh.In (b), the maximum area of the fracture is highlighted, while in (c) the black line indicating the fracture trace.
23(b) and 23(c).The global size is 100 × 100 × 100 m, with a fracture radius of about 60 m and an average element area of 2.4 m 2 .The discretized fracture surface is large enough to allow the propagation of the fracture without geometric constraints.The symmetric boundary conditions are imposed on the two symmetry plane, i.e., the boundaries parallel to x and y axes containing the well.The other two faces parallel to x and y axes are z-constrained.The material parameters are E = 30 GPa and ν = 0.25.The fracture frictional behavior is governed by Coulomb's criterion, characterized by θ = 30 • and zero cohesion.The fluid viscosity is µ l = 10 −9 MPa • s.The injection rate and the confining stress are Q 0 = 1 • 10 −2 m 3 /s and σ 0 = 10 MPa, respectively.Finally, we use the same value of initial conductivity as in the KGD example, i.e., C f,0 = 10 mD•m = 9.87•10 −15 m 2

Fig. 24 .
Fig.24.Results for the penny-shaped crack simulation.From the left to the right: fracture opening, fluid pressure, fracture radius and pressure behavior at the injection location.