<figure><img src="image_1.png" /></figure>

HAL
open science

# Analysis of a VEM-fully discrete polytopal scheme with bubble stabilisation for contact mechanics with Tresca friction

Jérôme Droniou, Ali Haidar, Roland Masson

► To cite this version:

Jérôme Droniou, Ali Haidar, Roland Masson. Analysis of a VEM-fully discrete polytopal scheme with bubble stabilisation for contact mechanics with Tresca friction. 2024. hal-04542460

HAL Id: hal-04542460

https://hal.science/hal-04542460

Preprint submitted on 12 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Analysis of a VEM-fully discrete polytopal scheme with bubble stabilisation for contact mechanics with Tresca friction

Jérôme Droniou, Ali Haidar, Roland Masson

► To cite this version:

Jérôme Droniou, Ali Haidar, Roland Masson. Analysis of a VEM-fully discrete polytopal scheme with bubble stabilisation for contact mechanics with Tresca friction. ESAIM: Mathematical Modelling and Numerical Analysis, 2025, 59 (2), pp.1043-1074. 10.1051/m2an/2025013. hal-04542460v2

HAL Id: hal-04542460

https://hal.science/hal-04542460v2

Submitted on 2 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Analysis of a VEM-fully discrete polytopal scheme with
bubble stabilisation for contact mechanics with Tresca friction

Jérôme Droniou*², Ali Haidar†¹,³, and Roland Masson‡¹

¹Université Côte d'Azur, Inria, Laboratoire J.A. Dieudonné, Nice, France

²IMAG, Univ. Montpellier, CNRS, Montpellier, France. School of Mathematics, Monash University, Australia

³IFP Energies nouvelles, department of mathematics, 92500 Rueil-Malmaison, France

Abstract

This work performs the convergence analysis of the polytopal nodal discretisation of contact-
mechanics (with Tresca friction) recently introduced in [18] in the framework of poro-elastic
models in fractured porous media. The scheme is based on a mixed formulation, using face-
wise constant approximations of the Lagrange multipliers along the fracture network and a fully
discrete first order nodal approximation of the displacement field. The displacement field is
enriched with additional bubble degrees of freedom along the fractures to ensure the inf-sup
stability with the Lagrange multiplier space. It is presented in a fully discrete formulation, which
makes its study more straightforward, but also has a Virtual Element interpretation.

The analysis establishes an abstract error estimate accounting for the fully discrete frame-
work and the non-conformity of the discretisation. A first order error estimate is deduced for
sufficiently smooth solutions both for the gradient of the displacement field and the Lagrange
multiplier. A key difficulty of the numerical analysis is the proof of a discrete inf-sup condition,
which is based on a non-standard $H^{-1/2}$-norm (to deal with fracture networks) and involves the
jump of the displacements, not their traces. The analysis also requires the proof of a discrete
Korn inequality for the discrete displacement field which takes into account fracture networks.
Numerical experiments based on analytical solutions confirm our theoretical findings.

**Keywords:** Contact-mechanics, fracture networks, polytopal method, fully discrete approach, virtual element method, bubble stabilisation, error estimates, discrete inf-sup condition, discrete Korn inequality.

# 1 Introduction

The simulation of poromechanical models in fractured (or faulted) porous rocks plays an important role in many subsurface applications such as the assessment of fault reactivation risks in CO2 storage or the hydraulic fracture stimulation in deep geothermal systems. These models couple the flow along the fractures and the surrounding matrix to the rock mechanical deformation and the mechanical behavior of the fractures. Fractures are classically represented as a network of planar surfaces connected to the surrounding matrix domain, leading to the so-called mixed-dimensional models which have been the object of many recent works in poromechanics [2, 5, 19, 22, 23, 18, 21, 24, 27].

Polytopal discretisations are motivated in subsurface applications to cope with the complexity of the
geometries representing geological structures including faults/fractures, layering, erosions and hetero-
geneities. Different classes of polytopal methods have been developed in the field of mechanics such

<footnote>* jerome.droniou@umontpellier.fr</footnote>

<footnote>† Corresponding author, ali.haidar@univ-cotedazur.fr</footnote>

<footnote>‡ Corresponding author, roland.masson@univ-cotedazur.fr</footnote>
Analysis of a VEM-fully discrete polytopal scheme with
bubble stabilisation for contact mechanics with Tresca friction

Jérôme Droniou*¹,², Ali Haidar†³,⁴, and Roland Masson‡³

¹IMAG, Univ. Montpellier, CNRS, Montpellier, France.

²School of Mathematics, Monash University, Australia

³Université Côte d'Azur, Inria, CNRS, Laboratoire J.A. Dieudonné, team Galets, Nice, France

⁴IFP Energies nouvelles, department of mathematics, 92500 Rueil-Malmaison, France

Abstract

This work performs the convergence analysis of the polytopal nodal discretisation of contact-
mechanics (with Tresca friction) recently introduced in [18] in the framework of poro-elastic
models in fractured porous media. The scheme is based on a mixed formulation, using face-
wise constant approximations of the Lagrange multipliers along the fracture network and a fully
discrete first order nodal approximation of the displacement field. The displacement field is
enriched with additional bubble degrees of freedom along the fractures to ensure the inf-sup
stability with the Lagrange multiplier space. It is presented in a fully discrete formulation, which
makes its study more straightforward, but also has a Virtual Element interpretation.

The analysis establishes an abstract error estimate accounting for the fully discrete frame-
work and the non-conformity of the discretisation. A first order error estimate is deduced for
sufficiently smooth solutions both for the gradient of the displacement field and the Lagrange
multiplier. A key difficulty of the numerical analysis is the proof of a discrete inf-sup condition,
which is based on a non-standard $H^{-1/2}$-norm (to deal with fracture networks) and involves the
jump of the displacements, not their traces. The analysis also requires the proof of a discrete
Korn inequality for the discrete displacement field which takes into account fracture networks.
Numerical experiments based on analytical solutions confirm our theoretical findings.

**Keywords:** Contact-mechanics, fracture networks, polytopal method, fully discrete approach, virtual element method, bubble stabilisation, error estimates, discrete inf-sup condition, discrete Korn inequality.

# 1 Introduction

The simulation of poromechanical models in fractured (or faulted) porous rocks plays an important role in many subsurface applications such as the assessment of fault reactivation risks in CO₂ storage or the hydraulic fracture stimulation in deep geothermal systems. These models couple the flow along the fractures and the surrounding matrix to the rock mechanical deformation and the mechanical behavior of the fractures. Fractures are classically represented as a network of planar surfaces connected to the surrounding matrix domain, leading to the so-called mixed-dimensional models which have been the object of many recent works in poromechanics [5–9, 21–23, 27, 29].

Polytopal discretisations are motivated in subsurface applications to cope with the complexity of
the geometries representing geological structures including faults/fractures, layering, erosions and

<footnote>* jerome.droniou@umontpellier.fr</footnote>

<footnote>† Corresponding author, ali.haidar@univ-cotedazur.fr, ali.haidar@aum.edu.kw</footnote>

<footnote>‡ Corresponding author, roland.masson@univ-cotedazur.fr</footnote>
as Discontinuous Galerkin [24], Hybrid High Order (HHO) [16], MultiPoint Stress Approximation (MPSA) [26], Hybrid Mimetic Methods [15] and Virtual Element Methods (VEM) [3, 14]. Some of them have been extended to account for contact-mechanics as in [5] for the MSPA based on facewise constant approximations of the surface tractions and displacement jump along the fracture network, in [12] for HHO combined with a Nitsche's contact formulation, and in [30] for VEM based on node to node contact conditions. Among these polytopal methods, VEM, as a natural extension of the Finite Element Method (FEM) to polyhedral meshes, has received a lot of attention in the mechanics community since its introduction in [3] and has been applied to various problems including in the context of geomechanics [1], poromechanics [13, 10, 20] and fracture mechanics [31].

In [18], an extension of the first order VEM to contact-mechanics is introduced based on a fully discrete framework with vector space of discrete unknowns and reconstruction operators in the spirit of Hybrid High Order discretisations [16]. Following [4, 19, 21, 7] in the FEM case, the contact problem is expressed in mixed form with face-wise constant Lagrange multipliers imposing the contact conditions in average on each face of the fracture network. This approach enables the handling of fracture networks (including corners, tips and intersections), the use of efficient semi-smooth Newton nonlinear solvers, and the preservation at the discrete level of the dissipative properties of the contact terms. On the other hand, the combination of a first order nodal discretisation of the displacement field with a face-wise constant approximation of the Lagrange multiplier requires a stabilisation to ensure the inf-sup compatibility condition. This is achieved in [18] by extending to the polytopal framework the $P^1$-bubble FEM discretisation [4] based on the enrichment of the displacement space by an additional bubble unknown on one side of each fracture face. Numerical evidence is obtained that this enrichment achieves the sought stabilisation, but no proof or analysis thereof is provided in this reference.

The goal of the present work is to precisely perform the numerical analysis of the polytopal discretisation proposed in [18]. To simplify the presentation, we focus on a static isotropic linear elastic mechanical model with Tresca frictional contact at matrix-fracture interfaces. The key new difficulty is related to the proof of the inf-sup condition between the discrete displacement and Lagrange multiplier spaces, which must account for the polytopal nature of the scheme and the geometrical complexity of the fracture network including tips, corners and intersections. Previous results [4, 19] established for FEM cannot be used since they are adapted neither to the fully discrete framework, nor to fracture networks which necessitates the introduction of a specific $H_{0,j}^{-1/2}(\Gamma)$-norm. The stability and convergence analysis also requires the proof of a discrete Korn inequality for the discrete displacement field accounting for the bubble stabilisation and the fracture network. The error estimate relies on an abstract estimate taking into account our fully discrete framework and the non-conformity of the displacement discretisation.

The remaining of this paper is organised as follows. Section 2 introduces the static contact-mechanical model with Tresca friction and its mixed formulation. Section 3 recalls the main ingredients of the discretisation from [18] with the mesh described in Section 3.1, the discrete spaces of displacement and Lagrange multiplier unknowns in Section 3.2, the function, jump and gradient reconstruction operators in Section 3.3, the definition of the interpolation operators in Section 3.4, and the discrete mixed formulation in Section 3.5. Section 4 states the main results regarding the well-posedness and convergence analysis of the scheme, while the proofs are reported in Section 5, with additional details in the Appendix. Section 6 investigates the numerical behavior of the scheme on analytical solutions in order to assess our theoretical results.
heterogeneities. Different classes of polytopal methods have been developed in the field of mechanics such as Discontinuous Galerkin [24], Hybrid High Order (HHO) [16], MultiPoint Stress Approximation (MPSA) [26], Hybrid Mimetic Mixed Methods [15] and Virtual Element Methods (VEM) [3, 14]. Some of them have been extended to account for contact-mechanics as in [5] for the MPSA based on face-wise constant approximations of the surface tractions and displacement jump along the fracture network, in [12] for HHO combined with a Nitsche contact formulation, and in [30] for VEM based on node to node contact conditions. Among these polytopal methods, VEM, as a natural extension of the Finite Element Method (FEM) to polyhedral meshes, has received a lot of attention in the mechanics community since its introduction in [3] and has been applied to various problems including in the context of geomechanics [1], poromechanics [10, 13, 20] and fracture mechanics [31].

In [18], an extension of the first order VEM to contact-mechanics is introduced based on a fully discrete framework with vector space of discrete unknowns and reconstruction operators in the spirit of Hybrid High Order discretisations [16]. Following [4, 7, 19, 21] in the FEM case, the contact problem is expressed in mixed form with face-wise constant Lagrange multipliers imposing the contact conditions in average on each face of the fracture network. This approach enables the handling of fracture networks (including corners, tips and intersections), the use of efficient semi-smooth Newton nonlinear solvers, and the preservation at the discrete level of the dissipative properties of the contact terms. On the other hand, the combination of a first order nodal discretisation of the displacement field with a face-wise constant approximation of the Lagrange multiplier requires a stabilisation to ensure the inf-sup compatibility condition. This is achieved in [18] by extending to the polytopal framework the $P^1$-bubble FEM discretisation [4] based on the enrichment of the displacement space with an additional bubble unknown on one side of each fracture face. Numerical evidence is obtained that this enrichment achieves the sought stabilisation, but no proof or analysis thereof is provided in this reference.

The goal of the present work is to precisely perform the numerical analysis of the polytopal discretisation proposed in [18]. To simplify the presentation, we focus on a static isotropic linear elastic mechanical model with Tresca frictional contact at matrix-fracture interfaces. The key new difficulty is related to the proof of the inf-sup condition between the discrete displacement and Lagrange multiplier spaces, which must account for the polytopal nature of the scheme and the geometrical complexity of the fracture network including tips, corners and intersections. Previous results [4, 19] established for FEM cannot be used since they are adapted neither to the fully discrete framework, nor to fracture networks which necessitates the introduction of a specific dual norm on the network. The stability and convergence analysis also require the proof of a discrete Korn inequality for the discrete displacement field accounting for the bubble stabilisation and the fracture network. The error estimate relies on an abstract estimate taking into account our fully discrete framework and the non-conformity of the displacement discretisation.

## 2 Model

We consider a Discrete Fracture Matrix (DFM) model on the polyhedral domain $\Omega \subset \mathbb{R}^d$ including a fracture network $\Gamma \subset \Omega$ with co-dimension 1. We assume that the network does not disconnect any part of $\Omega$ from the boundary of the domain, that is,

The intersection of $\partial\Omega$ with the boundary of each connected component of $\Omega \setminus \bar{\Gamma}$
is of non-zero $(d-1)$-dimensional measure. \hfill (1)
## 2 Model

We consider a Discrete Fracture Matrix (DFM) model on the polyhedral domain $\Omega \subset \mathbb{R}^d$ including a fracture network $\Gamma$ with co-dimension 1 defined by

$$ \bar{\Gamma} = \bigcup_{i \in I} \bar{\Gamma}_i . $$

We assume that $\Omega \setminus \Gamma$ is connected. Each fracture $\Gamma_i \subset \Omega$ ($i \in I$) is a polygonal simply connected open subdomain of a plane of $\mathbb{R}^d$. Without restriction of generality, it is assumed that fractures may only intersect at their boundaries. The two sides of a given fracture of $\Gamma$ are denoted by $\pm$ in the matrix domain $\Omega \setminus \bar{\Gamma}$. The two unit normal vectors $\mathbf{n}^{\pm}$, oriented outward from the sides $\pm$, satisfy $\mathbf{n}^{+} + \mathbf{n}^{-} = \mathbf{0}$. Given, for simplicity, homogeneous Dirichlet boundary conditions, the space for the displacement

$$ \mathbf{U}_0 = H_0^1(\Omega \setminus \bar{\Gamma})^d $$

is endowed with the norm $\|\mathbf{v}\|_{U_0} = \|\nabla \mathbf{v}\|_{L^2(\Omega \setminus \bar{\Gamma})}$ (which is, indeed, a norm since $\Omega \setminus \bar{\Gamma}$ is connected). The oriented jump operator on $\Gamma$ for functions $\mathbf{u} \in U_0$ is defined by

$$ [[\mathbf{u}] = \gamma^{+} \mathbf{u} - \gamma^{-} \mathbf{u}, $$

where $\gamma^{\pm}$ are the trace operators on both sides of $\Gamma$. Its normal and tangential components are denoted respectively by $[[\mathbf{u}]_n = [\mathbf{u}] \cdot \mathbf{n}^+$ and $[[\mathbf{u}]_{\tau} = [\mathbf{u}] - [[\mathbf{u}]_n \mathbf{n}^+]$. Note that $[[\mathbf{u}]_{\tau}$ depends on the orientation, while $[[\mathbf{u}]_n$ does not. The sided normal trace operator on $\Gamma$ oriented outward to the side $\pm$, applied to $H_{\text{div}}(\Omega \setminus \bar{\Gamma})$ functions, is denoted by $\gamma_n^{\pm}$. The symmetric gradient operator $\mathfrak{C}$ is defined on $U_0$ by $\mathfrak{C}(\mathbf{v}) = \frac{1}{2}(\nabla \mathbf{v} + t \nabla \mathbf{v})$.

The model we consider accounts for the mechanical equilibrium equation with a linear elastic constitutive law and a Tresca frictional contact model at matrix–fracture interfaces. In its strong form, it is written

$$ \left\{ \begin{array}{ll} -\operatorname{div}\sigma(\mathbf{u}) = \mathbf{f} & \text{on } \Omega \setminus \bar{\Gamma}, \\ \sigma(\mathbf{u}) = 2\mu\mathfrak{C}(\mathbf{u}) + \lambda\operatorname{div}\mathbf{u} & \text{on } \Omega \setminus \bar{\Gamma}, \\ \gamma_n^+\sigma(\mathbf{u}) + \gamma_n^-\sigma(\mathbf{u}) = \mathbf{0} & \text{on } \Gamma, \\ T_n(\mathbf{u}) \le 0, \quad [[\mathbf{u}]_n \le 0, \quad [[\mathbf{u}]_n T_n(\mathbf{u}) = 0 & \text{on } \Gamma, \\ |T_{\tau}(\mathbf{u})| \le g & \text{on } \Gamma, \\ T_{\tau}(\mathbf{u}) \cdot [[\mathbf{u}]_{\tau} + g|[[\mathbf{u}]_{\tau}| = 0 & \text{on } \Gamma, \end{array} \right. \quad (1) $$

with Tresca threshold $g \ge 0$, Lamé coefficients $\mu$ and $\lambda$, and normal and tangential surface tractions $T_n(\mathbf{u}) = \gamma_n^+\sigma(\mathbf{u}) \cdot \mathbf{n}^+$ and $T_{\tau}(\mathbf{u}) = \gamma_n^+\sigma(\mathbf{u}) - T_n(\mathbf{u})\mathbf{n}^+$. It is assumed in the following that the external force term $\mathbf{f}$ belongs to $L^2(\Omega)^d$.

The weak formulation of the mechanical model with Tresca frictional-contact (1) is written in mixed form using a vector-valued Lagrange multiplier $\lambda : \Gamma \to \mathbb{R}^d$ at matrix–fracture interfaces. Define the displacement jump space by

$$ H_{0,j}^{1/2}(\Gamma) = \{ [[\mathbf{v}] : \mathbf{v} \in U_0 \} $$

and denote by $H_{0,j}^{-1/2}(\Gamma)$ its dual space; the duality pairing between these two spaces is written $\langle \cdot, \cdot \rangle_{\Gamma}$. We will also use this notation for the duality pairing between $H_0^{1/2}(\Gamma)$, trace space of $H^1(\Omega \setminus \Gamma)$, and
We also assume that the network can be written as

$$ \bar{\Gamma} = \bigcup_{i \in I} \bar{\Gamma}_i, $$

where each fracture $\Gamma_i \subset \Omega$ ($i \in I$) is a polygonal simply connected open subset of a plane of $\mathbb{R}^d$. Without restriction of generality, the fractures are chosen such that any two of them may only intersect at their boundaries. The two sides of a given fracture of $\Gamma$ are denoted by $\pm$ in the matrix domain $\Omega \setminus \bar{\Gamma}$. The two unit normal vectors $\mathbf{n}^{\pm}$, oriented outward from the sides $\pm$, satisfy $\mathbf{n}^{+} + \mathbf{n}^{-} = \mathbf{0}$. Here and in the following, we use boldface for vectors, vector-valued functions, their spaces, and the corresponding discrete spaces; so, for example, $H^1(X) = H^1(X)^d$. Given, for simplicity, homogeneous Dirichlet boundary conditions, the space for the displacement is

$$ \mathbf{U}_0 = \{ \mathbf{v} \in H^1(\Omega \setminus \bar{\Gamma}) : \mathbf{v} = \mathbf{0} \text{ on } \partial\Omega \}. $$

It is endowed with the norm $\|\mathbf{v}\|_{\mathbf{U}_0} = \|\nabla \mathbf{v}\|_{L^2(\Omega \setminus \bar{\Gamma})}$ (the assumption (1) ensures that this is indeed a norm). The oriented jump operator on $\Gamma$ for functions $\mathbf{u} \in \mathbf{U}_0$ is defined by

$$ \llbracket \mathbf{u} \rrbracket = \gamma^{+} \mathbf{u} - \gamma^{-} \mathbf{u}, $$

where $\gamma^{\pm}$ are the trace operators on both sides of $\Gamma$. Its normal and tangential components are denoted respectively by $\llbracket \mathbf{u} \rrbracket_{\mathbf{n}} = \llbracket \mathbf{u} \rrbracket \cdot \mathbf{n}^{\pm}$ and $\llbracket \mathbf{u} \rrbracket_{\tau} = \llbracket \mathbf{u} \rrbracket - \llbracket \mathbf{u} \rrbracket_{\mathbf{n}} \mathbf{n}^{\pm}$. Note that $\llbracket \mathbf{u} \rrbracket_{\tau}$ depends on the orientation, while $\llbracket \mathbf{u} \rrbracket_{\mathbf{n}}$ does not. The sided normal trace operators on $\Gamma$ oriented outward from the side $\pm$, applied to $H_{\text{div}}(\Omega \setminus \bar{\Gamma})$ functions, are denoted by $\gamma_{\mathbf{n}}^{\pm}$. The symmetric gradient operator $\mathfrak{c}$ is defined on $\mathbf{U}_0$ by $\mathfrak{c}(\mathbf{v}) = \frac{1}{2}(\nabla \mathbf{v} + {}^t \nabla \mathbf{v})$.

**Remark 2.1** (Regularity of $\Omega \setminus \bar{\Gamma}$). We note that $\Omega \setminus \bar{\Gamma}$ is not a Lipschitz domain, and thus some properties commonly used on Sobolev spaces (density of smooth functions, etc.) do not hold for $\mathbf{U}_0$. However, this does not impact our analysis as we do not rely on such properties. The only practical impact, from our point of view, of the geometry of $\Omega \setminus \bar{\Gamma}$ is the existence of two different traces on each side of $\Gamma$, as mentioned above.

The model we consider accounts for the mechanical equilibrium equation with a linear elastic constitutive law and a Tresca frictional contact model at matrix–fracture interfaces. In its strong form, it is written

$$ \left\{ \begin{array}{ll} -\mathrm{div}\sigma(\mathbf{u}) = \mathbf{f} & \text{on } \Omega \setminus \bar{\Gamma}, \\ \sigma(\mathbf{u}) = \frac{E}{1+\nu}(\mathfrak{c}(\mathbf{u}) + \frac{\nu}{1-2\nu}(\mathrm{div}\,\mathbf{u})\mathbf{I}) & \text{on } \Omega \setminus \bar{\Gamma}, \\ \gamma_{\mathbf{n}}^{+}\sigma(\mathbf{u}) + \gamma_{\mathbf{n}}^{-}\sigma(\mathbf{u}) = \mathbf{0} & \text{on } \Gamma, \\ T_{\mathbf{n}}(\mathbf{u}) \le 0, \llbracket \mathbf{u} \rrbracket_{\mathbf{n}} \le 0, \llbracket \mathbf{u} \rrbracket_{\mathbf{n}} T_{\mathbf{n}}(\mathbf{u}) = 0 & \text{on } \Gamma, \\ |\mathbf{T}_{\tau}(\mathbf{u})| \le g & \text{on } \Gamma, \\ \mathbf{T}_{\tau}(\mathbf{u}) \cdot \llbracket \mathbf{u} \rrbracket_{\tau} + g |\llbracket \mathbf{u} \rrbracket_{\tau}| = 0 & \text{on } \Gamma, \\ \mathbf{u} = \mathbf{0} & \text{on } \partial\Omega, \end{array} \right. \quad (2) $$

with Tresca threshold $g \in L^2(\Gamma)$, $g \ge 0$, Young modulus $E$ and Poisson ratio $\nu$, and normal and tangential surface tractions $T_{\mathbf{n}}(\mathbf{u}) = \gamma_{\mathbf{n}}^{+}\sigma(\mathbf{u}) \cdot \mathbf{n}^{\pm}$ and $\mathbf{T}_{\tau}(\mathbf{u}) = \gamma_{\mathbf{n}}^{+}\sigma(\mathbf{u}) - T_{\mathbf{n}}(\mathbf{u})\mathbf{n}^{\pm}$. It is assumed in the following that the external force term $\mathbf{f}$ belongs to $L^2(\Omega)$.

The weak formulation of the mechanical model with Tresca frictional-contact (2) is written in mixed form using a vector-valued Lagrange multiplier $\lambda : \Gamma \to \mathbb{R}^d$ at matrix–fracture interfaces. Define the displacement jump space by

$$ \mathbf{W}_j(\Gamma) = \{ \llbracket \mathbf{v} \rrbracket : \mathbf{v} \in \mathbf{U}_0 \} $$
its dual space $H_0^{-1/2}(\Gamma)$. We note that $L^2(\Gamma)^d \subset H_{0,j}^{-1/2}(\Gamma)$ and that $\langle \mu, v \rangle_\Gamma = \int_\Gamma \mu \cdot v$ whenever $\mu \in L^2(\Omega)^d$. The dual cone is then defined by

$$ C_f = \{ \mu \in H_{0,j}^{-1/2}(\Gamma) : \langle \mu, v \rangle_\Gamma \le \langle g, |v_\tau| \rangle_\Gamma \text{ for all } v \in H_{0,j}^{1/2}(\Gamma) \text{ with } v \cdot n^+ \le 0 \}. $$

The weak mixed-variational formulation of (1) reads: find $\mathbf{u} \in \mathbf{U}_0$ and $\lambda \in C_f$ such that, for all $\mathbf{v} \in \mathbf{U}_0$ and $\mu \in C_f$,

$$ \int_{\Omega} \phi(\mathbf{u}) : \varepsilon(\mathbf{v}) + \langle \lambda, \llbracket \mathbf{v} \rrbracket \rangle_{\Gamma} = \int_{\Omega} \mathbf{f} \cdot \mathbf{v}, \quad (2a) $$

$$ \langle \mu - \lambda, \llbracket \mathbf{u} \rrbracket \rangle_{\Gamma} \le 0. \quad (2b) $$

Note that, based on the variational formulation, the Lagrange multiplier satisfies

$$ \lambda = -\gamma_{\mathbf{n}}^{+} \phi(\mathbf{u}) = \gamma_{\mathbf{n}}^{-} \phi(\mathbf{u}). $$

# 3 Scheme

## 3.1 Mesh

We take a polyhedral mesh of the domain $\Omega$ that is conforming with the fracture network $\Gamma$. We also assume that the Tresca threshold $g \ge 0$ is piecewise constant on the trace of the mesh on $\Gamma$. For each cell $K$ (resp. face $\sigma$), we denote by $h_K$ (resp. $h_\sigma$) and $|K|$ (resp. $|\sigma|$) its diameter and its measure, and we set

$$ h_{\mathcal{D}} = \max_{K \in \mathcal{M}} h_K. $$

The set of cells $K$, the set of faces $\sigma$, the set of nodes $s$ and the set of edges $e$ are denoted respectively by $\mathcal{M}$, $\mathcal{F}$, $\mathcal{V}$ and $\mathcal{E}$. By conformity of the mesh with respect to $\Gamma$, there exists a subset of faces $\mathcal{F}_\Gamma \subset \mathcal{F}$ such that

$$ \bar{\Gamma} = \bigcup_{\sigma \in \mathcal{F}_\Gamma} \bar{\sigma}. $$

We denote by $\mathcal{M}_\sigma$ the set of cells neighboring a face $\sigma \in \mathcal{F}$; thus, $\mathcal{M}_\sigma = \{K, L\}$ for interior face $\sigma \in \mathcal{F}^{\text{int}}$ (in which case we write $\sigma = K|L$) and $\mathcal{M}_\sigma = \{K\}$ for boundary face $\sigma \in \mathcal{F}^{\text{ext}}$. Since $\Gamma \subset \Omega$, we have $\mathcal{F}_\Gamma \subset \mathcal{F}^{\text{int}}$. For a face $\sigma \in \mathcal{F}_\Gamma$, $K$ and $L$ in the notation $\sigma = K|L$ are labelled such that $\mathbf{n}_{K\sigma} = \mathbf{n}^+$ and $\mathbf{n}_{L\sigma} = \mathbf{n}^-$, where $\mathbf{n}_{K\sigma}$ (resp. $\mathbf{n}_{L\sigma}$) is the unit normal vector to $\sigma$ oriented outward of $K$ (resp. $L$). We denote by $\mathcal{V}^{\text{ext}}$ the boundary nodes, and by $\mathcal{V}_\sigma$ the set of nodes of $\sigma$, $\mathcal{E}_\sigma$ the set of edges of $\sigma$, $\mathcal{F}_K$ the set of faces of $K$, $\mathcal{V}_K$ the set of nodes of $K$, and $\mathcal{M}_s$ the set of cells that contain the vertex $s$. For each $\sigma \in \mathcal{F}$, $\mathbf{n}_{\sigma e}$ is the unit normal vector to $e \in \mathcal{E}_\sigma$ in the plane $\sigma$ oriented outward to $\sigma$. For each $K \in \mathcal{M}$ and $\sigma \in \mathcal{F}_K$ we denote by $\gamma^{K\sigma}$ the trace operator on $\sigma$ for functions in $H^1(K)$ (or their vector-valued versions).

Throughout this paper we suppose that mesh regularity assumptions of [16, Definition 1.9] hold, and we write $a \lesssim b$ (resp. $a \ge b$) as a shorthand for $a \le Cb$ (resp. $Ca \le b$) with $C > 0$ depending only on $\Omega$, $\Gamma$, on the mesh regularity parameter, and possibly on the Lamé coefficients and $\mathbf{f}$.

If $X \in \mathcal{M} \cup \mathcal{F}$ and $\ell \in \mathbb{N}$, we denote by $\mathbb{P}^\ell(X)$ the space of polynomials of degree $\le \ell$ on $X$. For $X = \mathcal{M}$ or $X = \mathcal{F}_\Gamma$, we use the notation $\mathbb{P}^\ell(X)$ for the space of piecewise-polynomials of degree $\le \ell$ on $X$.
and denote by $W'_j(\Gamma)$ its dual space; the duality pairing between these two spaces is written $\langle \cdot, \cdot \rangle_\Gamma$.
We note that $L^2(\Gamma) \subset W'_j(\Gamma)$ and that $\langle \mu, w \rangle_\Gamma = \int_\Gamma \mu \cdot w$ whenever $\mu \in L^2(\Gamma)$ and $w \in W_j(\Gamma)$. The
space of admissible Lagrange multipliers is then defined by

$$ C_f = \{ \mu \in W'_j(\Gamma) : \langle \mu, w \rangle_\Gamma \le \int_\Gamma g |w_\tau| \text{ for all } w \in W_j(\Gamma) \text{ with } w \cdot n^+ \le 0 \}. $$

The weak mixed-variational formulation of (2) reads: find $\mathbf{u} \in \mathbf{U}_0$ and $\lambda \in C_f$ such that, for all $\mathbf{v} \in \mathbf{U}_0$ and $\mu \in C_f$,

$$ \int_{\Omega} \sigma(\mathbf{u}) : \mathcal{E}(\mathbf{v}) + \langle \lambda, \llbracket \mathbf{v} \rrbracket \rangle_{\Gamma} = \int_{\Omega} \mathbf{f} \cdot \mathbf{v}, \quad (3a) $$

$$ \langle \mu - \lambda, \llbracket \mathbf{u} \rrbracket \rangle_{\Gamma} \le 0. \quad (3b) $$

Note that, based on the variational formulation, the Lagrange multiplier satisfies

$$ \lambda = -\gamma_{\mathbf{n}}^{+} \sigma(\mathbf{u}) = \gamma_{\mathbf{n}}^{-} \sigma(\mathbf{u}). $$

# 3 Scheme

## 3.1 Mesh

We take a polytopal mesh of the domain $\Omega$ that is compliant with the fracture network $\Gamma$, in the sense that no mesh cell cuts $\Gamma$ (note that, in the framework of polytopal methods, this is actually not a very restrictive assumption, since each cell that would be across $\Gamma$ could simply be cut in two cells on each side of $\Gamma$ – possibly with some agglomeration strategy to avoid generating badly shaped cells). We also assume that the Tresca threshold $g \ge 0$ is piecewise constant on the trace of the mesh on $\Gamma$. For each cell $K$ (resp. face $\sigma$), we denote by $h_K$ (resp. $h_\sigma$) and $|K|$ (resp. $|\sigma|$) its diameter and its measure, and we set

$$ h_{\mathcal{D}} = \max_{K \in \mathcal{M}} h_K. $$

The set of cells $K$, the set of faces $\sigma$, the set of nodes $s$ and the set of edges $e$ are denoted respectively by $M$, $\mathcal{F}$, $\mathcal{V}$ and $\mathcal{E}$. By compliance of the mesh with $\Gamma$, there exists a subset of faces $\mathcal{F}_\Gamma \subset \mathcal{F}$ such that

$$ \bar{\Gamma} = \bigcup_{\sigma \in \mathcal{F}_\Gamma} \bar{\sigma}. $$

Let $\mathcal{F}^{\text{bd}}$ and $\mathcal{F}^{\text{in}}$ be respectively the sets of boundary and interior faces of $\mathcal{F}$ with respect to the domain $\Omega$. We denote by $M_\sigma$ the set of cells neighbouring a face $\sigma \in \mathcal{F}$; thus, $M_\sigma = \{K, L\}$ for an interior face $\sigma \in \mathcal{F}^{\text{in}}$ (in which case we write $\sigma = K|L$) and $M_\sigma = \{K\}$ for a boundary face $\sigma \in \mathcal{F}^{\text{bd}}$. Since $\Gamma \subset \Omega$, we have $\mathcal{F}_\Gamma \subset \mathcal{F}^{\text{in}}$. For a face $\sigma \in \mathcal{F}_\Gamma$, $K$ and $L$ in the notation $\sigma = K|L$ are labelled such that $\mathbf{n}_{K\sigma} = \mathbf{n}^+$ and $\mathbf{n}_{L\sigma} = \mathbf{n}^-$, where $\mathbf{n}_{K\sigma}$ (resp. $\mathbf{n}_{L\sigma}$) is the unit normal vector to $\sigma$ oriented outward from $K$ (resp. $L$). We denote by $\mathcal{V}^{\text{bd}}$ the boundary nodes, and by $\mathcal{V}_\Gamma$ the set of nodes on $\Gamma$, $\mathcal{V}_\sigma$ the set of nodes of $\sigma$, $\mathcal{E}_\sigma$ the set of edges of $\sigma$, $\mathcal{F}_K$ the set of faces of $K$, $\mathcal{V}_K$ the set of nodes of $K$, and $M_s$ the set of cells that contain the vertex $s$. For each $\sigma \in \mathcal{F}$, $\mathbf{n}_{\sigma e}$ is the unit normal vector to $e \in \mathcal{E}_\sigma$ in the plane $\sigma$ oriented outward from $\sigma$. For each $K \in \mathcal{M}$ and $\sigma \in \mathcal{F}_K$ we denote by $\gamma^{K\sigma}$ the trace operator on $\sigma$ for functions in $H^1(K)$ (or their vector-valued versions).

Throughout this paper we suppose that the mesh regularity assumptions of [16, Definition 1.9] hold, and we write $a \lesssim b$ (resp. $a \ge b$) as a shorthand for $a \le Cb$ (resp. $a \ge Cb$) with $C > 0$ depending
## 3.2 Spaces

The degrees of freedom (DOFs) for the displacement are nodal (attached to the vertices of the mesh). To account for the discontinuity of the discrete displacement field at matrix fracture interfaces, these nodal DOFs can be discontinuous across the fracture network – each vertex on the fracture network (that is not an interior tip) has two or more values attached to it, one for each local connected component of the matrix in a neighborhood of the vertex. To use specific notations, for each $s \in \mathcal{V}_\Gamma$ and each $K \in \mathcal{M}_s$, we denote by $v_{\mathcal{K}_s}$ the nodal unknown corresponding to the side of $K$ in the matrix. This unknown is identical for all cell lying on the same side as $K$: if $L \in \mathcal{M}_s$ is on the same side as $K$, then $v_{\mathcal{K}_s} = v_{\mathcal{L}_s}$. The notation $\mathcal{K}_s$ is used for the set of all cells containing $s$ and lying on the same side of $\Gamma$ as $K$. If $s \notin \mathcal{V}_\Gamma$, there is a unique displacement unknown at $s$, which is denoted by $v_{\mathcal{K}_s}$; in that case, $\mathcal{K}_s = \mathcal{M}_s$. Figure 1 visually illustrates this idea.

Additionally, on the positive side of the fracture, we attach to each face a vector DOF corresponding to a “correction” of the face values generated by the nodal DOFs; this correction plays the role, at the discrete level represented by the space of DOFs, of a bubble function and is introduced to ensure a suitable inf-sup condition between the space of Lagrange multipliers and the jump of the displacements.

<figure><img src="image_2.png" /><figcaption>Figure 1. Nodal degrees of freedom.</figcaption></figure>

For each $K \in \mathcal{M}$, we denote by

$$ \mathcal{F}_{\Gamma,K}^{+} = \{ \sigma \in \mathcal{F}_{\Gamma} \cap \mathcal{F}_{K} : \mathbf{n}_{K\sigma} \cdot \mathbf{n}^{+} > 0 \} $$

the set of faces of $K$ that lie on the positive side of $\Gamma$ (that set is empty if $K$ does not touch $\Gamma$ or only touches it from its negative side). According to the discussion above, the discrete space of displacement is

$$ \mathbf{U}_{0,\mathcal{D}} = \left\{ \mathbf{v}_{\mathcal{D}} = ((\mathbf{v}_{\mathcal{K}_s})_{K \in \mathcal{M}}, s \in \mathcal{V}_K, (\mathbf{v}_{K\sigma})_{K \in \mathcal{M}}, \sigma \in \mathcal{F}_{\Gamma,K}^{+}) : \right. \\ \mathbf{v}_{\mathcal{K}_s} \in \mathbb{R}^d, \mathbf{v}_{K\sigma} \in \mathbb{R}^d, \mathbf{v}_{\mathcal{K}_s} = 0 \text{ if } s \in \mathcal{V}^{\text{ext}} \\ \left. \mathbf{v}_{\mathcal{K}_s} = \mathbf{v}_{\mathcal{L}_s} \text{ if } K, L \in \mathcal{M}_s \text{ are on the same side of } \Gamma \right\}. $$

The Lagrange multiplier plays the role of approximations of $-\gamma_{\mathbf{n}}^{+}\sigma(\mathbf{u})$ on $\Gamma$. Its space is made of piecewise constant vectors:

$$ \mathbf{M}_{\mathcal{D}} = \{ \lambda_{\mathcal{D}} \in L^2(\Gamma)^d : \lambda_{\sigma} := (\lambda_{\mathcal{D}})|_{\sigma} \text{ is constant for all } \sigma \in \mathcal{F}_{\Gamma} \}. $$

For $\lambda_{\mathcal{D}} \in \mathbf{M}_{\mathcal{D}}$, we define its normal and tangential components by

$$ \lambda_{\mathcal{D},\mathbf{n}} = \lambda_{\mathcal{D}} \cdot \mathbf{n}^{+}, \quad \lambda_{\mathcal{D},\tau} = \lambda_{\mathcal{D}} - \lambda_{\mathcal{D},\mathbf{n}}\mathbf{n}^{+}, $$
only on $\Omega$, $\Gamma$, on the mesh regularity parameter, and possibly on the Young modulus and Poisson ratio and $f$.

If $X \in \mathcal{M} \cup \mathcal{F}$ and $\ell \in \mathbb{N}$, we denote by $\mathbb{P}^\ell(X)$ the space of polynomials of total degree $\le \ell$ on $X$. For $\mathcal{X} = \mathcal{M}$ or $\mathcal{X} = \mathcal{F}_\Gamma$, we use the notation $\mathbb{P}^\ell(\mathcal{X})$ for the space of piecewise-polynomials of degree $\le \ell$ on $\mathcal{X}$.

## 3.2 Spaces

The degrees of freedom (DOFs) for the displacement are nodal (attached to the vertices of the mesh). To account for the discontinuity of the discrete displacement field at matrix fracture interfaces, these nodal DOFs can be discontinuous across the fracture network – each vertex on the fracture network (that is not an interior tip) has two or more values attached to it, one for each local connected component of the matrix in a neighbourhood of the vertex. To use specific notations, for each $s \in \mathcal{V}_\Gamma$ and each $K \in \mathcal{M}_s$, we denote by $v_{\mathcal{K}_s}$ the nodal unknown corresponding to the side of $K$ in the matrix. This unknown is identical for all cell lying on the same side as $K$: if $L \in \mathcal{M}_s$ is on the same side as $K$, then $v_{\mathcal{K}_s} = v_{\mathcal{L}_s}$. The notation $\mathcal{K}_s$ is used for the set of all cells containing $s$ and lying on the same side of $\Gamma$ as $K$. If $s \notin \mathcal{V}_\Gamma$, there is a unique displacement unknown at $s$, which is denoted by $v_{\mathcal{K}_s}$; in that case, $\mathcal{K}_s = \mathcal{M}_s$. Figure 1 visually illustrates this idea.

<figure><img src="image_2.png" /><figcaption>Figure 1. Nodal degrees of freedom.</figcaption></figure>

Additionally, on the positive side of the fracture, we attach to each face a vector DOF corresponding to a “correction” of the face values generated by the nodal DOFs; this correction plays the role, at the discrete level represented by the space of DOFs, of a bubble function and is introduced to ensure a suitable inf-sup condition between the space of Lagrange multipliers and the jump of the displacements.

For each $K \in \mathcal{M}$, we denote by

$$ \mathcal{F}_{\Gamma,K}^{+} = \{ \sigma \in \mathcal{F}_{\Gamma} \cap \mathcal{F}_{K} : \mathbf{n}_{K\sigma} \cdot \mathbf{n}^{+} > 0 \} $$

the set of faces of $K$ that lie on the positive side of $\Gamma$ (that set is empty if $K$ does not touch $\Gamma$ or only touches it from its negative side). According to the discussion above, the discrete space of displacements is

$$ \mathbf{U}_{\mathcal{D}} = \left\{ \mathbf{v}_{\mathcal{D}} = ((\mathbf{v}_{\mathcal{K}_s})_{K \in \mathcal{M}}, s \in \mathcal{V}_K, (\mathbf{v}_{K\sigma})_{K \in \mathcal{M}}, \sigma \in \mathcal{F}_{\Gamma,K}^{+}) : \right. \\ \mathbf{v}_{\mathcal{K}_s} \in \mathbb{R}^d, \mathbf{v}_{K\sigma} \in \mathbb{R}^d, \\ \mathbf{v}_{\mathcal{K}_s} = \mathbf{v}_{\mathcal{L}_s} \text{ if } K, L \in \mathcal{M}_s \text{ are on the same side of } \Gamma \right\} $$
and the discrete dual cone by

$$C_D = \{ \lambda_D \in \mathbf{M}_D : \lambda_{D,n} \ge 0, |\lambda_{D,\tau}| \le g \} \subset C_f.$$

(the inclusion being easy to check).

### 3.3 Reconstruction operators in $U_{0,D}$

We first define, for each $K \in \mathcal{M}$ and $\sigma \in \mathcal{F}_K$, a tangential face gradient $\nabla^K \sigma : U_{0,D} \to \mathbb{P}^0(\sigma)^{d \times d}$ and tangential displacement reconstruction $\Pi^K \sigma : U_{0,D} \to \mathbb{P}^1(\sigma)^d$. First, we choose nonnegative weights $(\omega_s^\sigma)_{s \in \mathcal{V}_\sigma}$ to express the center of mass $\bar{\mathbf{x}}_\sigma$ of $\sigma$ in terms of its vertices:

$$\bar{\mathbf{x}}_{\sigma} = \sum_{s \in \mathcal{V}_{\sigma}} \omega_{s}^{\sigma} \mathbf{x}_{s}, \quad \sum_{s \in \mathcal{V}_{\sigma}} \omega_{s}^{\sigma} = 1. \qquad (4)$$

Then, for $\mathbf{v}_D \in U_{0,D}$, we set

$$\nabla^{K\sigma} \mathbf{v}_D = \frac{1}{|\sigma|} \sum_{e=s_1 s_2 \in \mathcal{E}_\sigma} |e| \frac{\mathbf{v}_{Ks_1} + \mathbf{v}_{Ks_2}}{2} \otimes \mathbf{n}_{\sigma e}, \qquad (5) \\ \Pi^{K\sigma} \mathbf{v}_D(\mathbf{x}) = \nabla^{K\sigma} \mathbf{v}_D(\mathbf{x} - \bar{\mathbf{x}}_\sigma) + \bar{\mathbf{v}}_{K\sigma} \quad \forall \mathbf{x} \in \sigma, \quad \text{where } \bar{\mathbf{v}}_{K\sigma} = \sum_{s \in \mathcal{V}_\sigma} \omega_s^\sigma \mathbf{v}_{Ks}.$$

Above, we have noted $e = s_1 s_2$ to indicate that the edge $e$ has vertices $s_1, s_2$.

If $\sigma \in \mathcal{F}_\Gamma$ is a fracture face, and $K$ (resp. $L$) is the cell on the positive (resp. negative) side of $\sigma$, we define the displacement jump operator on $\sigma$ as $[\cdot]_\sigma : U_{0,D} \to \mathbb{P}^0(\sigma)^d$ such that, for all $\mathbf{v}_D \in U_{0,D}$,

$$[\mathbf{v}_D]_\sigma = \frac{1}{|\sigma|} \int_\sigma (\Pi^{K\sigma} \mathbf{v}_D - \Pi^{L\sigma} \mathbf{v}_D) + \mathbf{v}_{K\sigma} = \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma} + \mathbf{v}_{K\sigma}. \qquad (6)$$

The normal component of that jump is denoted by $[\cdot]_{\sigma,n} = [\cdot]_\sigma \cdot \mathbf{n}_{K\sigma}$.

For each cell $K \in \mathcal{M}$, we select nonnegative weights $(\omega_s^K)_{s \in \mathcal{V}_K}$ of a linear decomposition of the center of mass $\bar{\mathbf{x}}_K$ of $K$ in terms of its vertices

$$\bar{\mathbf{x}}_K = \sum_{s \in \mathcal{V}_K} \omega_s^K \mathbf{x}_s, \quad \sum_{s \in \mathcal{V}_K} \omega_s^K = 1,$$

and we design a gradient reconstruction $\nabla^K : U_{0,D} \to \mathbb{P}^0(K)^{d \times d}$ and a displacement reconstruction $\Pi^K : U_{0,D} \to \mathbb{P}^1(K)^d$ by setting, for $\mathbf{v}_D \in U_{0,D}$,

$$\nabla^K \mathbf{v}_D = \frac{1}{|K|} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_{K\sigma} \otimes \mathbf{n}_{K\sigma} + \sum_{\sigma \in \mathcal{F}_{\Gamma,K}^+} \frac{|\sigma|}{|K|} \mathbf{v}_{K\sigma} \otimes \mathbf{n}_{K\sigma}, \qquad (7)$$

$$\Pi^K \mathbf{v}_D(\mathbf{x}) = \nabla^K \mathbf{v}_D(\mathbf{x} - \bar{\mathbf{x}}_K) + \bar{\mathbf{v}}_K \quad \forall \mathbf{x} \in K, \quad \text{where } \bar{\mathbf{v}}_K = \sum_{s \in \mathcal{V}_K} \omega_s^K \mathbf{v}_{Ks}. \qquad (8)$$

These local jump, gradient and displacement reconstructions are patched together to create their global piecewise polynomial counterparts $[\cdot]_D : U_{0,D} \to \mathbb{P}^0(\mathcal{F}_\Gamma)^d$, $\nabla^D : U_{0,D} \to \mathbb{P}^0(\mathcal{M})^{d \times d}$ and $\Pi^D : U_{0,D} \to \mathbb{P}^1(\mathcal{M})^d$: for all $\mathbf{v}_D \in U_{0,D}$,

$$([\mathbf{v}_D]_D)_\sigma = [\mathbf{v}_D]_\sigma \quad \forall \sigma \in \mathcal{F}_\Gamma,$$

$$(\nabla^D \mathbf{v}_D)_K = \nabla^K \mathbf{v}_D \quad \forall K \in \mathcal{M},$$
and its subspace incorporating homogeneous Dirichlet boundary conditions on $\partial\Omega$ is

$$ \mathbf{U}_{0,\mathcal{D}} = \{ \mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{\mathcal{D}} : \mathbf{v}_{\mathcal{K}_s} = \mathbf{0} \text{ if } s \in \mathcal{V}^{\text{bd}} \}. \quad (4) $$

The Lagrange multiplier plays the role of approximations of $-\gamma_{\mathbf{n}}^{+}\sigma(\mathbf{u})$ on $\Gamma$. Its space is made of piecewise constant vectors:

$$ \mathbf{M}_{\mathcal{D}} = \{ \lambda_{\mathcal{D}} \in L^2(\Gamma) : \lambda_{\sigma} := (\lambda_{\mathcal{D}})|_{\sigma} \text{ is constant for all } \sigma \in \mathcal{F}_{\Gamma} \}. $$

For $\lambda_{\mathcal{D}} \in \mathbf{M}_{\mathcal{D}}$, we define its normal and tangential components by

$$ \lambda_{\mathcal{D},\mathbf{n}} = \lambda_{\mathcal{D}} \cdot \mathbf{n}^{+}, \quad \lambda_{\mathcal{D},\tau} = \lambda_{\mathcal{D}} - \lambda_{\mathcal{D},\mathbf{n}}\mathbf{n}^{+}, $$

and the discrete dual cone by

$$ \mathbf{C}_{\mathcal{D}} = \{ \lambda_{\mathcal{D}} \in \mathbf{M}_{\mathcal{D}} : \lambda_{\mathcal{D},\mathbf{n}} \ge 0, |\lambda_{\mathcal{D},\tau}| \le g \} \subset \mathbf{C}_f. $$

The last inclusion is easy to check.

### 3.3 Reconstruction operators in $U_{0,\mathcal{D}}$

We first define, for each $K \in \mathcal{M}$ and $\sigma \in \mathcal{F}_K$, a tangential face gradient $\nabla^{K\sigma} : U_{0,\mathcal{D}} \to \mathbb{P}^0(\sigma)^{d \times d}$ and tangential displacement reconstruction $\Pi^{K\sigma} : U_{0,\mathcal{D}} \to \mathbb{P}^1(\sigma)^d$. First, we choose nonnegative weights $(\omega_s^\sigma)_{s \in \mathcal{V}_\sigma}$ to express the center of mass $\bar{\mathbf{x}}_\sigma$ of $\sigma$ in terms of its vertices:

$$ \bar{\mathbf{x}}_{\sigma} = \sum_{s \in \mathcal{V}_{\sigma}} \omega_s^{\sigma} \mathbf{x}_s, \quad \sum_{s \in \mathcal{V}_{\sigma}} \omega_s^{\sigma} = 1. \quad (5) $$

Then, for $\mathbf{v}_{\mathcal{D}} \in U_{0,\mathcal{D}}$, we set

$$ \nabla^{K\sigma} \mathbf{v}_{\mathcal{D}} = \frac{1}{|\sigma|} \sum_{e=s_1 s_2 \in \mathcal{E}_{\sigma}} |e| \frac{\mathbf{v}_{\mathcal{K} s_1} + \mathbf{v}_{\mathcal{K} s_2}}{2} \otimes \mathbf{n}_{\sigma e}, \quad (6) \\ \Pi^{K\sigma} \mathbf{v}_{\mathcal{D}}(\mathbf{x}) = \nabla^{K\sigma} \mathbf{v}_{\mathcal{D}}(\mathbf{x} - \bar{\mathbf{x}}_{\sigma}) + \bar{\mathbf{v}}_{K\sigma} \quad \forall \mathbf{x} \in \sigma, \quad \text{where } \bar{\mathbf{v}}_{K\sigma} = \sum_{s \in \mathcal{V}_{\sigma}} \omega_s^{\sigma} \mathbf{v}_{\mathcal{K} s}. $$

Above, we have noted $e = s_1 s_2$ to indicate that the edge $e$ has vertices $s_1, s_2$, and $\otimes$ is the tensor product of two vectors (that is, $\mathbf{a} \otimes \mathbf{b}$ is the matrix such that $(\mathbf{a} \otimes \mathbf{b})_{i,j} = \mathbf{a}_i \mathbf{b}_j$).

If $\sigma \in \mathcal{F}_{\Gamma}$ is a fracture face, and $K$ (resp. $L$) is the cell on the positive (resp. negative) side of $\sigma$, we define the displacement jump operator on $\sigma$ as $[\cdot]_{\sigma} : U_{0,\mathcal{D}} \to \mathbb{P}^0(\sigma)^d$ such that, for all $\mathbf{v}_{\mathcal{D}} \in U_{0,\mathcal{D}}$,

$$ [\mathbf{v}_{\mathcal{D}}]_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} (\Pi^{K\sigma} \mathbf{v}_{\mathcal{D}} - \Pi^{L\sigma} \mathbf{v}_{\mathcal{D}}) + \mathbf{v}_{K\sigma} = \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma} + \mathbf{v}_{K\sigma}. \quad (7) $$

The normal component of this jump is denoted by $[\cdot]_{\sigma,\mathbf{n}} = [\cdot]_{\sigma} \cdot \mathbf{n}_{K\sigma}$.

For each cell $K \in \mathcal{M}$, we select nonnegative weights $(\omega_s^K)_{s \in \mathcal{V}_K}$ of a linear decomposition of the center of mass $\bar{\mathbf{x}}_K$ of $K$ in terms of its vertices

$$ \bar{\mathbf{x}}_K = \sum_{s \in \mathcal{V}_K} \omega_s^K \mathbf{x}_s, \quad \sum_{s \in \mathcal{V}_K} \omega_s^K = 1, $$
$$ (\Pi^{\mathcal{D}}\mathbf{v}_{\mathcal{D}})|_{K} = \Pi^{K}\mathbf{v}_{\mathcal{D}} \quad \forall K \in \mathcal{M}. $$

We also define the cellwise constant reconstruction operator $\tilde{\Pi}^{\mathcal{D}}\mathbf{v}_{\mathcal{D}} : \mathbf{U}_{0,\mathcal{D}} \to \mathbb{P}^{0}(\mathcal{M})^{d}$ such that $(\tilde{\Pi}^{\mathcal{D}}\mathbf{v}_{\mathcal{D}})|_{K} = \bar{\mathbf{v}}_{K}$. Finally, the discrete symmetric gradient $\mathcal{E}_{\mathcal{D}}$, divergence $\mathrm{div}_{\mathcal{D}}$ and stress tensor $\sigma_{\mathcal{D}}$ are deduced from the previous operators:

$$ \mathcal{E}_{\mathcal{D}} = \frac{1}{2}(\nabla^{\mathcal{D}} + \mathbf{T}^{\mathcal{D}}), \quad \mathrm{div}_{\mathcal{D}} = \mathrm{Tr}(\mathcal{E}_{\mathcal{D}}) \quad \text{and} \quad \sigma_{\mathcal{D}}(\cdot) = 2\mu\mathcal{E}_{\mathcal{D}}(\cdot) + \lambda\mathrm{div}_{\mathcal{D}}(\cdot)\mathbb{I}. $$

**Remark 3.1** (Non planar faces). At this point, we have only considered meshes with planar faces. In numerical simulations in geosciences, meshes with non-planar faces are however naturally encountered, for example in Corner Pointed Geometries (CPG) situations. One approach to handle such meshes is to cut the non-planar faces into two or more planar subfaces; this leads to polytopal meshes that can be handled by our method. Another approach was introduced in [11]. It consists in adding a barycentric center $\mathbf{x}_{\sigma} = \frac{1}{|\mathcal{H}|} \sum_{s \in \mathcal{V}_{\sigma}} \mathbf{x}_{s}$ to the non-planar face $\sigma$, which is artificial in the sense that it is not counted as a geometric node of the mesh. The discrete displacement $\mathbf{u}_{\sigma}$ at $\mathbf{x}_{\sigma}$ is then defined by barycentric linear combination of the nodal displacements at the vertices of the face. Additionally, a triangulation of the face $\sigma$, centered at $\mathbf{x}_{\sigma}$, is defined. In this context, the gradient reconstruction operator (7) in each cell must be adapted as follows

$$ \nabla^K \mathbf{v}_D = \frac{1}{|K|} \sum_{\sigma \in \mathcal{T}_K} \sum_{e = s_1 s_2 \in \mathcal{E}_\sigma} \frac{|T_e|}{3} (\mathbf{v}_{Ks_1} + \mathbf{v}_{Ks_2} + \mathbf{v}_\sigma) \otimes \mathbf{n}_{KT_e} + \sum_{\sigma \in \mathcal{T}_K^+} \frac{|\sigma|}{|K|} \mathbf{v}_{K\sigma} \otimes \mathbf{n}_{K\sigma}, $$

where $T_e$ is the triangle defined by $\mathbf{x}_{\sigma}$ and the edge $e$, $\mathbf{n}_{KT_e}$ is the normal vector to triangle $T_e$ pointing out of $K$, and $\mathbf{v}_{\sigma} = \frac{1}{|\mathcal{H}|} \sum_{s \in \mathcal{V}_{\sigma}} \mathbf{v}_{Ks}$. Note that the fracture faces are still assumed to be planar, hence the bubble terms are unchanged. Moreover, for cells with planar faces, this new gradient is identical to (7). With this approach, the same discrete space of displacement $\mathbf{U}_{0,\mathcal{D}}$ is retained, as the unknowns remain unchanged.

## 3.4 Interpolators

The space $C_0^0(\bar{\Omega}\setminus\Gamma)$ is spanned by functions that are continuous on $\bar{\Omega}\setminus\Gamma$, have limits on each side of $\Gamma$, and vanish on $\partial\Omega$. The interpolator $\tilde{I}_{U_{0,\mathcal{D}}}: C_0^0(\bar{\Omega}\setminus\Gamma)^d \to \mathbf{U}_{0,\mathcal{D}}$ is defined through its components by setting, for $\mathbf{v} \in C_0^0(\bar{\Omega}\setminus\Gamma)^d$,

$$ \begin{aligned} (\tilde{I}_{U_{0,\mathcal{D}}} \mathbf{v})_{Ks} &= \mathbf{v}|_K(\mathbf{x}_s) & \forall K \in \mathcal{M}, \forall s \in \mathcal{V}_K, \\ (\tilde{I}_{U_{0,\mathcal{D}}} \mathbf{v})_{K\sigma} &= \frac{1}{|\sigma|} \int_{\sigma} (\gamma^{K\sigma} \mathbf{v} - \Pi^{K\sigma} (\tilde{I}_{U_{0,\mathcal{D}}} \mathbf{v})) & \forall K \in \mathcal{M}, \forall \sigma \in \mathcal{F}_{\Gamma,K}^{+}. \end{aligned} \quad (9) $$

This definition is seemingly recursive, since the DOF corresponding to $K\sigma$ are built using the interpolator $\tilde{I}_{U_{0,\mathcal{D}}}$ itself. However, the definition (5) of $\Pi^{K\sigma}$ shows that $\Pi^{K\sigma}(\tilde{I}_{U_{0,\mathcal{D}}} \mathbf{v})$ only depends on the nodes unknowns $(\tilde{I}_{U_{0,\mathcal{D}}} \mathbf{v})_{Ks}$, which are properly defined without any self-reference to $\tilde{I}_{U_{0,\mathcal{D}}}$. We also note that, since $\mathbf{v} = 0$ on $\partial\Omega$, this operator indeed defines an element in $\mathbf{U}_{0,\mathcal{D}}$.

The interpolator $\tilde{I}_{M_{\mathcal{D}}}: L^2(\Gamma) \to M_{\mathcal{D}}$ on the Lagrange space simply corresponds to averaging on each face: for $\lambda \in L^2(\Gamma)^d$,

$$ (\tilde{I}_{M_{\mathcal{D}}} \lambda)_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} \lambda \quad \forall \sigma \in \mathcal{F}_{\Gamma}. $$
and we design a gradient reconstruction $\nabla^K : \mathbf{U}_{0,\mathcal{D}} \to \mathbb{P}^0(\mathcal{K})^{d \times d}$ and a displacement reconstruction $\Pi^K : \mathbf{U}_{0,\mathcal{D}} \to \mathbb{P}^1(\mathcal{K})^d$ by setting, for $\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}$,

$$ \nabla^K \mathbf{v}_{\mathcal{D}} = \frac{1}{|K|} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_{K\sigma} \otimes \mathbf{n}_{K\sigma} + \sum_{\sigma \in \mathcal{F}_{\Gamma,K}^+} \frac{|\sigma|}{|K|} \mathbf{v}_{K\sigma} \otimes \mathbf{n}_{K\sigma}, \quad (8) $$

$$ \Pi^K \mathbf{v}_{\mathcal{D}}(\mathbf{x}) = \nabla^K \mathbf{v}_{\mathcal{D}}(\mathbf{x} - \bar{\mathbf{x}}_K) + \bar{\mathbf{v}}_K \quad \forall \mathbf{x} \in K, \quad \text{where } \bar{\mathbf{v}}_K = \sum_{s \in \mathcal{V}_K} \omega_s^K \mathbf{v}_{Ks}. \quad (9) $$

These local jump, gradient and displacement reconstructions are patched together to create their global piecewise polynomial counterparts $[\cdot]_{\mathcal{D}} : \mathbf{U}_{0,\mathcal{D}} \to \mathbb{P}^0(\mathcal{F}_{\Gamma})^d$, $\nabla^{\mathcal{D}} : \mathbf{U}_{0,\mathcal{D}} \to \mathbb{P}^0(\mathcal{M})^{d \times d}$ and $\Pi^{\mathcal{D}} : \mathbf{U}_{0,\mathcal{D}} \to \mathbb{P}^1(\mathcal{M})^d$. We additionally define the cell-wise constant reconstruction operator $\tilde{\Pi}^{\mathcal{D}} : \mathbf{U}_{0,\mathcal{D}} \to \mathbb{P}^0(\mathcal{M})^d$ as the projection of $\Pi^{\mathcal{D}}$ on piecewise constant functions. These operators therefore satisfy: for all $\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}$,

$$ (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket)_{|\sigma} = \llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\sigma} \quad \forall \sigma \in \mathcal{F}_{\Gamma}, $$

$$ (\nabla^{\mathcal{D}} \mathbf{v}_{\mathcal{D}})_{|K} = \nabla^K \mathbf{v}_{\mathcal{D}} \quad \forall K \in \mathcal{M}, $$

$$ (\Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}})_{|K} = \Pi^K \mathbf{v}_{\mathcal{D}} \quad \forall K \in \mathcal{M}, $$

$$ (\tilde{\Pi}^{\mathcal{D}} \mathbf{v}_{\mathcal{D}})_{|K} = \bar{\mathbf{v}}_K \quad \forall K \in \mathcal{M}. $$

Finally, the discrete symmetric gradient $\epsilon^{\mathcal{D}}$, divergence $\mathrm{div}^{\mathcal{D}}$ and stress tensor $\sigma^{\mathcal{D}}$ are deduced from the previous operators: with Tr the matrix trace operator,

$$ \epsilon^{\mathcal{D}} = \frac{1}{2}(\nabla^{\mathcal{D}} + {}^t\nabla^{\mathcal{D}}), \quad \mathrm{div}^{\mathcal{D}} = \mathrm{Tr}(\epsilon^{\mathcal{D}}) \quad \text{and} \quad \sigma^{\mathcal{D}}(\cdot) = \frac{E}{1+\nu} \left( \epsilon^{\mathcal{D}}(\cdot) + \frac{\nu}{1-2\nu} \mathrm{div}^{\mathcal{D}}(\cdot) \mathbb{I} \right). $$

**Remark 3.1** (Non planar faces). Up to this point, we have only considered meshes with planar faces. In numerical simulations in geosciences, meshes with non-planar faces are however naturally encountered, for example in Corner Pointed Geometries (CPG) situations. One approach to handle such meshes is to cut the non-planar faces into two or more planar subfaces; this leads to polytopal meshes that can be handled by our method. Another approach was introduced in [11]. It consists in adding a barycentric center $\mathbf{x}_{\sigma} = \frac{1}{|\mathcal{F}_{\sigma}|} \sum_{s \in \mathcal{V}_{\sigma}} \mathbf{x}_s$ to the non-planar face $\sigma$, which is artificial in the sense that it is not counted as a geometric node of the mesh. The discrete displacement $\mathbf{u}_{\sigma}$ at $\mathbf{x}_{\sigma}$ is then defined as the barycentric linear combination of the nodal displacements at the vertices of the face. Additionally, a triangulation of the face $\sigma$, centered at $\mathbf{x}_{\sigma}$, is defined. In this context, the gradient reconstruction operator (8) in each cell must be adapted as follows

$$ \nabla^K \mathbf{v}_{\mathcal{D}} = \frac{1}{|K|} \sum_{\sigma \in \mathcal{F}_K} \sum_{e = s_1 s_2 \in \mathcal{E}_{\sigma}} \frac{|T_e|}{3} (\mathbf{v}_{Ks_1} + \mathbf{v}_{Ks_2} + \mathbf{v}_{\sigma}) \otimes \mathbf{n}_{KT_e} + \sum_{\sigma \in \mathcal{F}_{\Gamma,K}^+} \frac{|\sigma|}{|K|} \mathbf{v}_{K\sigma} \otimes \mathbf{n}_{K\sigma}, $$

where $T_e$ is the triangle defined by $\mathbf{x}_{\sigma}$ and the edge $e$, $\mathbf{n}_{KT_e}$ is the normal vector to triangle $T_e$ pointing out of $K$, and $\mathbf{v}_{\sigma} = \frac{1}{|\mathcal{F}_{\sigma}|} \sum_{s \in \mathcal{V}_{\sigma}} \mathbf{v}_{Ks}$. Note that the fracture faces are still assumed to be planar, hence the bubble terms are unchanged. Moreover, for cells with planar faces, this new gradient is identical to (8). With this approach, the same discrete space of displacements $\mathbf{U}_{0,\mathcal{D}}$ is retained, as the unknowns remain unchanged.
## 3.5 Mixed-variational formulation

We now introduce the numerical scheme for the mixed-variational formulation of the mechanics contact problem (2): Find $(\mathbf{u}_{\mathcal{D}}, \lambda_{\mathcal{D}}) \in \mathbf{U}_{0,\mathcal{D}} \times \mathbf{C}_{\mathcal{D}}$ such that, for all $(\mathbf{v}_{\mathcal{D}}, \mu_{\mathcal{D}}) \in \mathbf{U}_{0,\mathcal{D}} \times \mathbf{C}_{\mathcal{D}}$,

$$ \int_{\Omega} \sigma_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}) : \mathcal{E}_{\mathcal{D}}(\mathbf{v}_{\mathcal{D}}) + S_{\mu,\lambda,\mathcal{D}}(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) + \int_{\Gamma} \lambda_{\mathcal{D}} \cdot [[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}} = \int_{\Omega} \mathbf{f} \cdot \tilde{\Pi}^{\mathcal{D}} \mathbf{v}_{\mathcal{D}}, \quad (10a) $$

$$ \int_{\Gamma} (\mu_{\mathcal{D}} - \lambda_{\mathcal{D}}) \cdot [[\mathbf{u}_{\mathcal{D}}]]_{\mathcal{D}} \le 0. \quad (10b) $$

Here, $S_{\mu,\lambda,\mathcal{D}}$ is the scaled stabilisation bilinear form defined by

$$ S_{\mu,\lambda,\mathcal{D}}(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) = \sum_{K \in \mathcal{M}} (2\mu_K + \lambda_K)S_K(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) $$

where, for each $K \in \mathcal{M}$, $\mu_K$ and $\lambda_K$ are respectively the average of $\lambda$ and $\mu$ on $K$ and the local stabilisation bilinear form $S_K: \mathbf{U}_{0,\mathcal{D}} \times \mathbf{U}_{0,\mathcal{D}} \to \mathbb{R}$ is given by

$$ S_K(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) = h_K^{d-2} \sum_{s \in \mathcal{V}_K} (\mathbf{u}_{\mathcal{K}s} - \Pi^K \mathbf{u}_{\mathcal{D}}(\mathbf{x}_s)) \cdot (\mathbf{v}_{\mathcal{K}s} - \Pi^K \mathbf{v}_{\mathcal{D}}(\mathbf{x}_s)) \\ + h_K^{d-2} \sum_{\sigma \in \mathcal{F}_{\Gamma,K}^{+}} \mathbf{u}_{K\sigma} \cdot \mathbf{v}_{K\sigma}. \quad (11) $$

Let us also introduce the unscaled stabilisation bilinear form

$$ S_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) = \sum_{K \in \mathcal{M}} S_K(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}). $$

*Remark 3.2.* Thanks to the fracture face-wise Lagrange multiplier, the variational inequality (10b) together with $\lambda_{\mathcal{D}} \in \mathbf{C}_{\mathcal{D}}$ can be equivalently replaced by the following non linear equations (see [2]), in which $[s]_{\mathbb{R}^+} = \max(0, s)$, and $[\xi]_g$ is the projection of $\xi \in \mathbb{R}^d$ on the ball of center 0 and radius g:

$$ \left\{ \begin{array}{l} \lambda_{\mathcal{D},\mathbf{n}} = [\lambda_{\mathcal{D},\mathbf{n}} + \beta_{\mathcal{D},\mathbf{n}} [[\mathbf{u}_{\mathcal{D}}]]_{\mathcal{D},\mathbf{n}}]_{\mathbb{R}^+} \\ \lambda_{\mathcal{D},\tau} = [\lambda_{\mathcal{D},\tau} + \beta_{\mathcal{D},\tau} [[\mathbf{u}_{\mathcal{D}}]]_{\mathcal{D},\tau}]_g \end{array} \right. \quad (12) $$

(these equations can easily be expressed locally to each fracture face). Here, $\beta_{\mathcal{D},\mathbf{n}} > 0$, $\beta_{\mathcal{D},\tau} > 0$ are given face-wise constant functions along $\Gamma$.

*Remark 3.3* (Virtual element interpretation). It is shown in [18] that the scheme (10) can be re-interpreted in a Virtual Element presentation.

## 4 Main results

To carry out the convergence analysis of the scheme, we need to define norms on the spaces of unknowns: an $H^1$-like norm on the space of displacement and an $H_{0,j}^{-1/2}(\Gamma)$-like norm on the space of Lagrange multipliers.
## 3.4 Interpolators

The space $C_0^0(\bar{\Omega}\setminus\Gamma)$ is spanned by functions that are continuous on $\bar{\Omega}\setminus\Gamma$, have limits on each side of $\Gamma$, and vanish on $\partial\Omega$. The interpolator $\mathcal{I}_{U_{0,\mathcal{D}}}: C_0^0(\bar{\Omega}\setminus\Gamma) \to U_{0,\mathcal{D}}$ is defined by setting, for $\mathbf{v} \in C_0^0(\bar{\Omega}\setminus\Gamma)$,

$$ (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{v})_{\mathcal{K}\sigma} = \mathbf{v}_{|\mathcal{K}}(\mathbf{x}_s) \qquad \forall \mathcal{K} \in \mathcal{M}, \forall s \in \mathcal{V}_K, \\ (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{v})_{K\sigma} = \frac{1}{|\sigma|} \int_{\sigma} (\gamma^{K\sigma} \mathbf{v} - \Pi^{K\sigma} (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{v})) \qquad \forall \mathcal{K} \in \mathcal{M}, \forall \sigma \in \mathcal{F}_{\Gamma,K}^+. \tag{10} $$

This definition is seemingly recursive, since the DOF corresponding to $K\sigma$ are built using the interpolator $\mathcal{I}_{U_{0,\mathcal{D}}}$ itself. However, the definition (6) of $\Pi^{K\sigma}$ shows that $\Pi^{K\sigma}(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{v})$ only depends on the nodal unknowns $(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{v})_{\mathcal{K}\sigma}$, which are properly defined without any self-reference to $\mathcal{I}_{U_{0,\mathcal{D}}}$. We also note that, since $\mathbf{v} = \mathbf{0}$ on $\partial\Omega$, this operator indeed defines an element in $U_{0,\mathcal{D}}$.

The interpolator $\mathcal{I}_{M_{\mathcal{D}}}: L^2(\Gamma) \to M_{\mathcal{D}}$ on the Lagrange multiplier space simply corresponds to averaging on each face: for $\lambda \in L^2(\Gamma)$,

$$ (\mathcal{I}_{M_{\mathcal{D}}} \lambda)_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} \lambda \qquad \forall \sigma \in \mathcal{F}_{\Gamma}. $$

## 3.5 Mixed-variational formulation

We now introduce the numerical scheme for the mixed-variational formulation of the contact mechanics problem (3): Find $(\mathbf{u}_{\mathcal{D}}, \lambda_{\mathcal{D}}) \in U_{0,\mathcal{D}} \times C_{\mathcal{D}}$ such that, for all $(\mathbf{v}_{\mathcal{D}}, \mu_{\mathcal{D}}) \in U_{0,\mathcal{D}} \times C_{\mathcal{D}}$,

$$ \int_{\Omega} \phi^{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}) : \mathcal{E}^{\mathcal{D}}(\mathbf{v}_{\mathcal{D}}) + S_{E,\nu,\mathcal{D}}(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) + \int_{\Gamma} \lambda_{\mathcal{D}} \cdot [[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}} = \int_{\Omega} \mathbf{f} \cdot \tilde{\Pi}^{\mathcal{D}} \mathbf{v}_{\mathcal{D}}, \tag{11a} $$

$$ \int_{\Gamma} (\mu_{\mathcal{D}} - \lambda_{\mathcal{D}}) \cdot [[\mathbf{u}_{\mathcal{D}}]]_{\mathcal{D}} \le 0. \tag{11b} $$

Here, $S_{E,\nu,\mathcal{D}}$ is the scaled stabilisation bilinear form defined by

$$ S_{E,\nu,\mathcal{D}}(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) = \sum_{K \in \mathcal{M}} \zeta(E_K, \nu_K) S_K(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) $$

with local stabilisation bilinear form $S_K: U_{0,\mathcal{D}} \times U_{0,\mathcal{D}} \to \mathbb{R}$ given by

$$ S_K(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) = h_K^{d-2} \sum_{s \in \mathcal{V}_K} \left( \mathbf{u}_{\mathcal{K}s} - \Pi^K \mathbf{u}_{\mathcal{D}}(\mathbf{x}_s) \right) \cdot \left( \mathbf{v}_{\mathcal{K}s} - \Pi^K \mathbf{v}_{\mathcal{D}}(\mathbf{x}_s) \right) \\ + h_K^{d-2} \sum_{\sigma \in \mathcal{F}_{\Gamma,K}^+} \mathbf{u}_{K\sigma} \cdot \mathbf{v}_{K\sigma}. \tag{12} $$

Above, for each $K \in \mathcal{M}$, $E_K$ and $\nu_K$ respectively denote the average of $E$ and $\nu$ on $K$, and $\zeta$ is such that $\zeta(E, \nu) > 0$ scales as $\frac{E}{1-2\nu}$ (or, equivalently, as the sum of the cell-averaged Lamé coefficients).

Let us also introduce the unscaled stabilisation bilinear form

$$ S_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) = \sum_{K \in \mathcal{M}} S_K(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}). $$
**Definition 4.1** (Discrete $H^1$-like semi-norm on $U_D$). The semi-norm $\|\cdot\|_{1,D}$ on $U_D$ is defined by: for all $u_D \in U_D$,

$$ \|u_D\|_{1,D} := \left( \sum_{K \in \mathcal{M}} \|u_D\|_{1,K}^2 \right)^{1/2} \quad \text{with} \quad \|u_D\|_{1,K} = \left( \|\nabla^K u_D\|_{L^2(K)}^2 + S_K(u_D, u_D) \right)^{1/2}, \quad (13) $$

where $S_K$ is given by (11). We note that $\|\cdot\|_{1,D}$ is genuinely a norm on $U_{0,D}$ since $\Omega \setminus \bar{\Gamma}$ is connected.

To define the norm on $M_D$, we recall that $\Gamma = \cup_{i \in I} \Gamma_i$ with each $\Gamma_i$ open connected subset of a hyperplane, and $\Gamma_i \cap \Gamma_j = \emptyset$ if $i \neq j$. We define $\Omega_i^+$ as the intersection of $\Omega$ with the half-plane defined by $\Gamma_i$ and $(\mathbf{n})|_{\Gamma_i}$; see Figure 2. The space $H^1(\Omega_i^+; \Gamma_i)$ is spanned by functions in $H^1(\Omega_i^+)$ that vanish on $\partial \Omega_i^+ \setminus \Gamma_i$.

<figure><img src="image_3.png" /><figcaption>Figure 2. Splitting of the fracture network and construction of sub-domains used to define the norm on $M_D$.</figcaption></figure>

**Definition 4.2** ($H_{0,j}^{-1/2}(\Gamma)$-like norm on $M_D$). The $H_{0,j}^{-1/2}(\Gamma)$-like norm on $M_D$ is defined by: for all $\lambda_D \in M_D$,

$$ \|\lambda_D\|_{-1/2, \Gamma} = \sum_{i \in I} \|\lambda_D\|_{-1/2, \Gamma_i} \quad \text{with} \quad \|\lambda_D\|_{-1/2, \Gamma_i} = \sup_{v_i \in H^1(\Omega_i^+; \Gamma_i)^d \setminus \{0\}} \frac{\int_{\Gamma_i} \lambda_D \cdot v_i}{\|v_i\|_{H^1(\Omega_i^+)}}. \quad (14) $$

**Remark 4.3** (Norm on the fracture network). The definition of this $H_{0,j}^{-1/2}(\Gamma)$-like norm is non-standard due to localisation to each planar component of $\Gamma$. This is however required to properly take into account the possible complex topology of the fracture network (triple – or more – intersections of fractures, etc.) and the fact that the inf-sup condition below is based on the jump of the functions, not their traces. See Remark 5.6 for more insight on this.

In case of a simple network in which no more than two planar fractures intersect at a given location, since these intersections (and the fracture tips) have a zero 2-capacity it can be checked that the norm (14) is equivalent to a more standard $H_0^{-1/2}(\Gamma)$-norm on the network.

**Proposition 4.4** (Existence and uniqueness result). There exists a unique solution $(u_D, \lambda_D) \in U_{0,D} \times C_D$ to (10).
Remark 3.2. Thanks to the fracture face-wise Lagrange multiplier, the variational inequality (11b) together with $\lambda_D \in C_D$ can be equivalently replaced by the following non linear equations (see [2]), in which $[s]_{\mathbb{R}^+} = \max(0, s)$, and $[\xi]_g$ is the projection of $\xi \in \mathbb{R}^d$ on the ball of center 0 and radius g:

$$ \left\{ \begin{aligned} \lambda_{D,n} &= [\lambda_{D,n} + \beta_{D,n} \llbracket \mathbf{u}_D \rrbracket]_{D,n}]_{\mathbb{R}^+} \\ \lambda_{D,\tau} &= [\lambda_{D,\tau} + \beta_{D,\tau} \llbracket \mathbf{u}_D \rrbracket]_{D,\tau}]_g \end{aligned} \right. \quad (13) $$

(these equations can easily be expressed locally to each fracture face). Here, $\beta_{D,n} > 0$, $\beta_{D,\tau} > 0$ are given face-wise constant functions along $\Gamma$.

Remark 3.3 (Virtual element interpretation). It is shown in [18] that the scheme (11) can be re-interpreted in a Virtual Element presentation.

Remark 3.4 (Incompressible limit). Note that an enrichment of the space of displacements with face bubbles in the normal direction on all faces of the mesh should provide a robust discretisation in the incompressible limit $\nu \to \frac{1}{2}$. This enrichment goes beyond the present work and remains to be investigated.

# 4 Main results

The following proposition states the well-posedness of the scheme (11).

**Proposition 4.1** (Existence and uniqueness of the solution to the scheme). There exists a unique solution $(\mathbf{u}_D, \lambda_D) \in U_{0,D} \times C_D$ to (11).

*Proof.* The abstract framework of [25, Theorem 3.8] on saddle-point problems gives the existence of a solution. To prove the uniqueness, let us define the discrete energy inner product $\langle \cdot, \cdot \rangle_{e,D}$ such that, for $\mathbf{u}_D, \mathbf{v}_D \in U_{0,D}$,

$$ \langle \mathbf{u}_D, \mathbf{v}_D \rangle_{e,D} = \int_{\Omega} \sigma^D(\mathbf{u}_D) : \varepsilon^D(\mathbf{v}_D) + S_{E,\nu,D}(\mathbf{u}_D, \mathbf{v}_D), \quad (14) $$

which is indeed an inner product owing to the discrete Korn inequality (45) below. Assume that (11) admits two solutions denoted by $(\mathbf{u}_D^{(1)}, \lambda_D^{(1)})$ and $(\mathbf{u}_D^{(2)}, \lambda_D^{(2)})$. Subtracting the variational equalities (11a) obtained with both solutions gives, for all $\mathbf{v}_D \in U_{0,D}$,

$$ \langle \mathbf{u}_D^{(1)} - \mathbf{u}_D^{(2)}, \mathbf{v}_D \rangle_{e,D} = \int_{\Gamma} (\lambda_D^{(2)} - \lambda_D^{(1)}) \cdot \llbracket \mathbf{v}_D \rrbracket_D. \quad (15) $$

Setting $\mathbf{v}_D = \mathbf{u}_D^{(1)} - \mathbf{u}_D^{(2)}$ in (15), we obtain

$$ \langle \mathbf{v}_D, \mathbf{v}_D \rangle_{e,D} = \int_{\Gamma} (\lambda_D^{(2)} - \lambda_D^{(1)}) \cdot \llbracket \mathbf{u}_D^{(1)} \rrbracket_D + \int_{\Gamma} (\lambda_D^{(1)} - \lambda_D^{(2)}) \cdot \llbracket \mathbf{u}_D^{(2)} \rrbracket_D. $$

Setting $\mu_D = \lambda_D^{(2)}$ in (11b) for the first solution and $\mu_D = \lambda_D^{(1)}$ in (11b) for the second solution, we see that the right-hand-side of the above equation is non-positive, which proves that $\mathbf{v}_D = \mathbf{0}$, that is, $\mathbf{u}_D^{(1)} = \mathbf{u}_D^{(2)}$. Then, (15) yields $\int_{\Gamma} (\lambda_D^{(2)} - \lambda_D^{(1)}) \cdot \llbracket \mathbf{v}_D \rrbracket_D = 0$ for all $\mathbf{v}_D \in U_{0,D}$, which implies $\lambda_D^{(1)} = \lambda_D^{(2)}$ by the discrete inf-sup property (23) below. $\square$
Proof. [25, Theorem 3.9] gives the existence of a solution. The uniqueness of $\mathbf{u}_D$ derives from the discrete Korn inequality (42) below and the uniqueness of $\lambda_D$ from the discrete inf-sup property (20) below. $\square$

To state the error estimates, we introduce the following notations:

• $\pi_{\bar{\Gamma}}^0$ is the orthogonal projection on $\mathbb{P}^0(\mathcal{F}_\Gamma)$.

• The (primal) consistency error is: for $\mathbf{v}_D \in \mathbf{U}_{0, D}$,

$$C_D(\mathbf{u}, \mathbf{v}_D) = \left( \|\nabla \mathbf{u} - \nabla^\mathcal{D} \mathbf{v}_D\|_{L^2(\Omega \setminus \bar{\Gamma})}^2 + S_D(\mathbf{v}_D, \mathbf{v}_D) \right)^{1/2}. \quad (15)$$

• Letting

$$\mathbf{W} = \{\sigma \in H_{\text{div}}(\Omega \setminus \bar{\Gamma}; \mathcal{S}^d(\mathbb{R})), \gamma_{\mathbf{n}}^{+}\sigma + \gamma_{\mathbf{n}}^{-}\sigma = \mathbf{0}, \gamma_{\mathbf{n}}^{+}\sigma \in L^2(\Gamma)^d\}, \quad (16)$$

(where $\mathcal{S}^d(\mathbb{R})$ is the space of symmetric $d \times d$ matrices with real coefficients), the adjoint consistency error (or limit-conformity measure) is defined, for $\sigma \in \mathbf{W}$, by

$$\mathcal{W}_D(\sigma) = \sup_{\mathbf{v}_D \in \mathbf{U}_{0,D}} \frac{w_D(\sigma, \mathbf{v}_D)}{\|\mathbf{v}_D\|_{1,D}}, \quad (17)$$

$$\text{where } w_D(\sigma, \mathbf{v}_D) = - \int_{\Omega} \sigma : \mathcal{C}_D(\mathbf{v}_D) + \int_{\Gamma} \gamma_{\mathbf{n}}^{+} \sigma \cdot [[\mathbf{v}_D]]_D - \int_{\Omega} \tilde{\Pi}^{\mathcal{D}} \mathbf{v}_D \cdot \text{div} \sigma.$$

• The discrete $H_0^{1/2}(\Gamma)$-norm is defined on $L^2(\Gamma)^d$ by: for $\mu \in L^2(\Gamma)^d$,

$$\|\mu\|_{1/2, \mathcal{D}} = \left( \sum_{\sigma \in \mathcal{F}_\Gamma} h_{\sigma}^{-1} \|\mu\|_{L^2(\sigma)}^2 \right)^{1/2}. \quad (18)$$

**Theorem 4.5 (Abstract error estimate).** For $(\mathbf{u}_D, \lambda_D)$ solution of (10) and $(\mathbf{u}, \lambda)$ solution of (2) with $\lambda \in L^2(\Gamma)^d$, we have the following abstract error estimate

$$\|\nabla^{\mathcal{D}} \mathbf{u}_D - \nabla \mathbf{u}\|_{L^2(\Omega \setminus \bar{\Gamma})} + \|\lambda_D - \lambda\|_{-1/2, \Gamma} \lesssim \|\lambda - \pi_{\bar{\Gamma}}^0 \lambda\|_{-1/2, \Gamma} + \mathcal{W}_D(\sigma(\mathbf{u})) \\ + \inf_{\mathbf{v}_D \in \mathbf{U}_{0, \mathcal{D}}} \left\{ \begin{aligned} & \|\lambda - \pi_{\bar{\Gamma}}^0 \lambda\|_{L^2(\Gamma)} \|[[\mathbf{v}_D]]_D - [[\mathbf{u}]\|_{L^2(\Gamma)}\|^{1/2} \\ & + \|[[\mathbf{v}_D]]_D - \pi_{\bar{\Gamma}}^0 [[\mathbf{u}]\|_{1/2, \mathcal{D}} + C_D(\mathbf{u}, \mathbf{v}_D)\|\end{aligned} \right\}. \quad (19)$$

Proof. See Section 5.3. $\square$

In the following theorem, we denote by $H^2(\mathcal{M})$ (resp. $H^1(\mathcal{F}_\Gamma)$) the space of functions defined on $\Omega$ that are $H^2$ on each $K \in \mathcal{M}$ (resp. defined on $\Gamma$ and $H^1$ on each $\sigma \in \mathcal{F}_\Gamma$). These spaces are endowed with their usual broken semi-norms.

**Theorem 4.6 (Error estimate).** Let $(\mathbf{u}, \lambda)$ be the solution to (2) and assume that $\mathbf{u} \in H^2(\mathcal{M})$ and $\lambda \in H^1(\mathcal{F}_\Gamma)$. Then the solution $(\mathbf{u}_D, \lambda_D)$ of (10) satisfies the following error estimate:

$$\|\nabla^{\mathcal{D}} \mathbf{u}_D - \nabla \mathbf{u}\|_{L^2(\Omega \setminus \bar{\Gamma})} + \|\lambda_D - \lambda\|_{-1/2, \Gamma} \lesssim h_D (|\lambda|_{H^1(\mathcal{F}_\Gamma)} + |\mathbf{u}|_{H^2(\mathcal{M})} + \|[[\mathbf{u}]\|_{H^1(\mathcal{F}_\Gamma)}).$$

Proof. See Section 5.4. $\square$
To carry out the convergence analysis of the scheme, we need to define norms on the spaces of unknowns: an $H^1$-like norm on the space of displacements and an $H^{-1/2}$-like norm on the space of Lagrange multipliers.

**Definition 4.2** (Discrete $H^1$-like semi-norm on $U_D$). The semi-norm $\|\cdot\|_{1,D}$ on $U_D$ is defined by: for all $u_D \in U_D$,

$$ \|u_D\|_{1,D} := \left( \sum_{K \in \mathcal{M}} \|u_D\|_{1,K}^2 \right)^{1/2} \quad \text{with} \quad \|u_D\|_{1,K} = \left( \|\nabla^K u_D\|_{L^2(K)}^2 + S_K(u_D, u_D) \right)^{1/2}, \quad (16) $$

where $S_K$ is given by (12). We note that $\|\cdot\|_{1,D}$ is actually a norm on $U_{0,D}$ by (1).

To define the norm on $M_D$, we recall that $\bar{\Gamma} = \cup_{i \in I} \Gamma_i$ with each $\Gamma_i$ an open simply connected subset of an hyperplane, and $\Gamma_i \cap \Gamma_j = \emptyset$ if $i \neq j$. We define $\Omega_i^+$ as the intersection of $\Omega$ with the half-space on the + side of $\Gamma_i$; see Figure 2. The space $H_0^1(\Omega_i^+; \Gamma_i)$ is spanned by vector-valued functions in $H^1(\Omega_i^+)$ that vanish on $\partial \Omega_i^+ \setminus \Gamma_i$.

<figure><img src="image_3.png" /><figcaption>Figure 2. Splitting of the fracture network and construction of sub-domains used to define the norm on $M_D$.</figcaption></figure>

**Definition 4.3** ($H^{-1/2}$-like norm on $L^2(\Gamma)^d$). The norm $\|\cdot\|_{-1/2,\Gamma}$ on $L^2(\Gamma)^d$ is defined by: for all $\lambda \in L^2(\Gamma)^d$,

$$ \|\lambda\|_{-1/2, \Gamma} = \sum_{i \in I} \|\lambda\|_{-1/2, \Gamma_i} \quad \text{with} \quad \|\lambda\|_{-1/2, \Gamma_i} = \sup_{v_i \in H_0^1(\Omega_i^+; \Gamma_i) \setminus \{0\}} \frac{\int_{\Gamma_i} \lambda \cdot v_i}{\|v_i\|_{H^1(\Omega_i^+) \setminus \{0\}}} \quad (17) $$

**Remark 4.4** (Norm on the fracture network). The definition of this $H^{-1/2}$-like norm is non-standard due to localisation to each planar component of $\Gamma$. This is however required to properly take into account the possible complex topology of the fracture network (triple – or more – intersections of fractures, etc.) and the fact that the inf-sup condition below is based on the jump of the functions, not their traces. See Remark 5.6 for more insight on this.

In the case of a simple network in which no more than two planar fractures intersect at a given location, since these intersections (and the fracture tips) have a zero 2-capacity in $\Omega$ it can be checked that the norm (17) is equivalent to a more standard $H^{-1/2}$-norm on the network.

To state the error estimates, we introduce the following notations:

* $\pi_{\bar{\Gamma}}^0$ is the orthogonal projection on $\mathbb{P}^0(\bar{\mathcal{F}}_\Gamma)^d$.
# 5 Proof of the error estimate

## 5.1 Inf-sup condition

We prove in this section the following result, which establishes that the bubble degree of freedom in (3) plays its role controlling the Lagrange multiplier through the jump of the displacements.

**Theorem 5.1 (Discrete inf-sup condition).** It holds

$$ \sup_{\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}} \setminus \{0\}} \frac{\int_{\Gamma} \lambda_{\mathcal{D}} \cdot [[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}}}{\|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}} \geq \|\lambda_{\mathcal{D}}\|_{-1/2,\Gamma} \quad \forall \lambda_{\mathcal{D}} \in \mathbf{M}_{\mathcal{D}}. \qquad (20) $$

### 5.1.1 Fracture-compatible averaged interpolator

To prove the inf-sup condition, we need a Clément-like interpolator that takes into account the fracture network. Due to the design of the $H_0^{1/2}(\Gamma)$-like norm, we will use this interpolator considering only one planar component at a time. In the following, we therefore fix $i \in I$ and, in the constructions below, the sets $\mathcal{K}s(i)$ are considered only in respect to $\Gamma_i$. Thus, for all $K \in \mathcal{M}$ and $s \in \mathcal{V}_K$, $\mathcal{K}s(i)$ is the set of all cells $L \in \mathcal{M}_s$ that lie on the same side of $\Gamma_i$ as $K$. In particular, $\mathcal{K}s \subset \mathcal{K}s(i)$ with equality whenever $s$ is internal to $\Gamma_i$.

For each $K \in \mathcal{M}$ and $s \in \mathcal{V}_K$, we take an open set $U_{\mathcal{K}s(i)} \subset \cup_{K \in \mathcal{K}s(i)} K$ and a function $\omega_{\mathcal{K}s(i)} \in L^\infty(U_{\mathcal{K}s(i)})$ such that

$$ |U_{\mathcal{K}s(i)}| \geq \max_{K \in \mathcal{K}s(i)} |K|, \quad |\omega_{\mathcal{K}s(i)}| \leq 1, \qquad (21a) $$

$$ \frac{1}{|U_{\mathcal{K}s(i)}|} \int_{U_{\mathcal{K}s(i)}} \omega_{\mathcal{K}s(i)} = 1, \quad \frac{1}{|U_{\mathcal{K}s(i)}|} \int_{U_{\mathcal{K}s(i)}} \mathbf{x}\omega_{\mathcal{K}s(i)} = \mathbf{x}_s. \qquad (21b) $$

If $s \notin \Gamma$, $U_{\mathcal{K}s(i)}$ can be taken as a ball centered at $s$, and $\omega_{\mathcal{K}s(i)} = 1$. Figure 3 illustrates possible choices for $U_{\mathcal{K}s(i)}$ depending on the nature of $s$, and Appendix A presents an explicit way to construct $(U_{\mathcal{K}s(i)}, \omega_{\mathcal{K}s(i)})$ in the generic case.

<figure><img src="image_4.png" /><figcaption>Figure 3. Domains for local fracture-compatible averages.</figcaption></figure>

The space $H_0^1(\Omega \setminus \Gamma_i)$ is the subspace of $H^1(\Omega \setminus \Gamma_i)$ spanned by functions that vanish on $\partial\Omega$, but not necessarily on $\Gamma_i$. For $\mathbf{v} \in H_0^1(\Omega \setminus \Gamma_i)^d$, the averaged interpolate $\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v} \in \mathbf{U}_{0,\mathcal{D}}$ of $\mathbf{v}$ is defined by

$$ (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s} = 0 \qquad \text{for all } s \in \mathcal{V}^{\text{ext}}, \qquad (22a) $$
• The (primal) consistency error is: for $\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}$,

$$C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}}) = \left( \|\nabla \mathbf{u} - \nabla^{\mathcal{D}} \mathbf{v}_{\mathcal{D}}\|_{L^2(\Omega \setminus \bar{\Gamma})}^2 + S_{\mathcal{D}}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) \right)^{1/2}. \quad (18)$$

• Letting

$$\Sigma = \{\sigma \in H_{\text{div}}(\Omega \setminus \bar{\Gamma}; S^d(\mathbb{R})) : \gamma_{\mathbf{n}}^{+}\sigma + \gamma_{\mathbf{n}}^{-}\sigma = \mathbf{0} \text{ on } \Gamma, \gamma_{\mathbf{n}}^{+}\sigma \in L^2(\Gamma)\}, \quad (19)$$

where $S^d(\mathbb{R})$ is the space of symmetric $d \times d$ matrices with real coefficients, the adjoint consistency error (or limit-conformity measure) is defined, for $\sigma \in \Sigma$, by

$$\mathcal{W}_{\mathcal{D}}(\sigma) = \sup_{\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}} \frac{w_{\mathcal{D}}(\sigma, \mathbf{v}_{\mathcal{D}})}{\|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}}, \quad (20)$$

$$\text{where } w_{\mathcal{D}}(\sigma, \mathbf{v}_{\mathcal{D}}) = - \int_{\Omega} \sigma : \mathcal{C}^{\mathcal{D}}(\mathbf{v}_{\mathcal{D}}) + \int_{\Gamma} \gamma_{\mathbf{n}}^{+}\sigma \cdot [[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}} - \int_{\Omega} \tilde{\Pi}^{\mathcal{D}}\mathbf{v}_{\mathcal{D}} \cdot \mathbf{div}\sigma.$$

• The discrete $H^{1/2}$-norm is defined on $L^2(\Gamma)$ by: for $\mu \in L^2(\Gamma)$,

$$\|\mu\|_{1/2, \mathcal{D}} = \left( \sum_{\sigma \in \mathcal{F}_{\Gamma}} h_{\sigma}^{-1} \|\mu\|_{L^2(\sigma)}^2 \right)^{1/2}. \quad (21)$$

**Theorem 4.5 (Abstract error estimate).** For $(\mathbf{u}_{\mathcal{D}}, \lambda_{\mathcal{D}})$ solution of (11) and $(\mathbf{u}, \lambda)$ solution of (3) with $\lambda \in L^2(\Gamma)$, we have the following abstract error estimate:

$$\|\nabla^{\mathcal{D}}\mathbf{u}_{\mathcal{D}} - \nabla\mathbf{u}\|_{L^2(\Omega \setminus \bar{\Gamma})} + \|\lambda_{\mathcal{D}} - \lambda\|_{-1/2, \Gamma} \lesssim \|\lambda - \pi_{\mathcal{F}_{\Gamma}}^{0}\lambda\|_{-1/2, \Gamma} + \mathcal{W}_{\mathcal{D}}(\sigma(\mathbf{u})) \\ + \inf_{\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}} \left\{ \left( \|\lambda - \pi_{\mathcal{F}_{\Gamma}}^{0}\lambda\|_{L^2(\Gamma)} \|[[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}} - [[\mathbf{u}]\|_{L^2(\Gamma)} \right)^{1/2} \right. \\ \left. + \|[[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}} - \pi_{\mathcal{F}_{\Gamma}}^{0}\|[[\mathbf{u}]\|_{1/2, \mathcal{D}} + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}}) \right\}. \quad (22)$$

*Proof.* See Section 5.3. □

In the following theorem, we denote by $H^2(\mathcal{M})$ (resp. $H^1(\mathcal{F}_{\Gamma})$) the space of vector-valued functions defined on $\Omega$ that are $H^2$ on each $K \in \mathcal{M}$ (resp. defined on $\Gamma$ and $H^1$ on each $\sigma \in \mathcal{F}_{\Gamma}$). These spaces are endowed with their usual broken semi-norms.

**Theorem 4.6 (Error estimate).** Let $(\mathbf{u}, \lambda)$ be the solution to (3) and assume that $\mathbf{u} \in H^2(\mathcal{M})$ and $\lambda \in H^1(\mathcal{F}_{\Gamma})$. Then the solution $(\mathbf{u}_{\mathcal{D}}, \lambda_{\mathcal{D}})$ of (11) satisfies the following error estimate:

$$\|\nabla^{\mathcal{D}}\mathbf{u}_{\mathcal{D}} - \nabla\mathbf{u}\|_{L^2(\Omega \setminus \bar{\Gamma})} + \|\lambda_{\mathcal{D}} - \lambda\|_{-1/2, \Gamma} \lesssim h_{\mathcal{D}}(|\lambda|_{H^1(\mathcal{F}_{\Gamma})} + |\mathbf{u}|_{H^2(\mathcal{M})} + |[[\mathbf{u}]|_{H^1(\mathcal{F}_{\Gamma})}]).$$

*Proof.* See Section 5.4. □
$$ (\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s} = \frac{1}{|U_{\mathcal{K}s(i)}|} \int_{U_{\mathcal{K}s(i)}} \bar{\omega}_{\mathcal{K}s(i)} \mathbf{v} \quad \text{for all } K \in \mathcal{M}, s \in \mathcal{V}_K \setminus \mathcal{V}^{\text{ext}}, \qquad (22b) $$

$$ (\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{K\sigma} = \frac{1}{|\sigma|} \int_{\sigma} (\gamma^{K\sigma} \mathbf{v} - \Pi^{K\sigma} (\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v})) \quad \text{for all } K \in \mathcal{M}, \sigma \in \mathcal{F}_{\Gamma_i, K}^{+}, \qquad (22c) $$

$$ (\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{K\sigma} = 0 \quad \text{for all } K \in \mathcal{M}, \sigma \in \mathcal{F}_{\Gamma, K}^{+} \setminus \mathcal{F}_{\Gamma, K}^{+}. \qquad (22d) $$

**Remark 5.2** (Averaged interpolator). The usage of averages in the definition of the nodal values is mandatory since we need well-defined and stable interpolations of functions with only $H^1$-regularity – see Proposition 5.5. However, the averages must not be done across $\Gamma_i$ (hence the condition $U_{\mathcal{K}s(i)} \subset \cup_{K \in \mathcal{K}s(i)} K$), and must have $\mathbf{x}_s$ as center of mass (hence the integral conditions on $\bar{\omega}_{\mathcal{K}s(i)}$), to ensure that the interpolator is exact on linear functions – see the proof of Proposition 5.5.

If $s \in \Gamma$ is not in the relative interior of $\Gamma_i$, then for all $K, L \in \mathcal{M}_s$ we have $\mathcal{K}s(i) = \mathcal{L}s(i)$ and thus $(\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s} = (\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{L}s}$. In other words, this $\Gamma_i$-adapted interpolator $\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a}$ produces single nodal values on $\Gamma_j$ for $j \neq i$, and possibly multiple nodal values only on vertices in the interior of $\Gamma_i$. This is coherent with the continuity properties of functions in $H_0^1(\Omega \setminus \Gamma_i)$ that it interpolates.

Finally, we note that a zero value is imposed for boundary nodes, so that $\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v} \in \mathbf{U}_{0,\mathcal{D}}$ (computing boundary nodal values by averaging would not ensure that they vanish), and that the “bubble” value is set at zero for faces not on $\Gamma_i$.

The main properties of this interpolator are its behaviour with respect to the jump, and its $H^1$-stability.

**Proposition 5.3** (Jump of the averaged interpolate). For all $\mathbf{v} \in H_0^1(\Omega \setminus \Gamma_i)^d$, it holds

$$ \forall K \in \mathcal{M}, \forall \sigma \in \mathcal{F}_K, \quad \int_{\sigma} [\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v}]_{\sigma} = \begin{cases} \int_{\sigma} (\gamma^{K\sigma} \mathbf{v} - \Pi^{L\sigma} \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v}) & \text{if } \sigma \in \mathcal{F}_{\Gamma_i, K}^{+}, \\ 0 & \text{if } \sigma \in \mathcal{F}_{\Gamma, K}^{+} \setminus \mathcal{F}_{\Gamma_i, K}^{+}. \end{cases} \qquad (23) $$

*Proof.* If $\sigma \in \mathcal{F}_{\Gamma_i, K}^{+}$, then the equality directly follows from the definition (6) of the jump on $\sigma$ and from (22c). If $\sigma \in \mathcal{F}_{\Gamma, K}^{+} \setminus \mathcal{F}_{\Gamma_i, K}^{+}$ then any $s \in \mathcal{V}_{\sigma}$ is not internal to $\Gamma_i$ and thus, as noted in Remark 5.2, for all $K, L \in \mathcal{M}_s$ we have $(\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s} = (\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{L}s}$. Hence, $\Pi^{K\sigma} \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v} = \Pi^{L\sigma} \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v}$. The definition (6) of the jump on $\sigma$, together with the fact that $(\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{K\sigma} = 0$ by (22d) yields the relation in (23). $\square$

Establishing the $H^1$-stability of $\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}}^{i,a}$ requires the following lemma.

**Lemma 5.4** (DOF-based bound on the discrete norm). Let $K \in \mathcal{M}$. Recalling the definition (13) of $\|\cdot\|_{1,K}$, we have, for all $\mathbf{w}_K = ((\mathbf{w}_{\mathcal{K}s})_{s \in \mathcal{V}_K}, (\mathbf{w}_{K\sigma})_{\sigma \in \mathcal{F}_{\Gamma, K}^{+}})$,

$$ \|\mathbf{w}_K\|_{1,K} \lesssim h_K^{-1} |K|^{1/2} \left( \max_{s \in \mathcal{V}_K} |\mathbf{w}_{\mathcal{K}s}| + \max_{\sigma \in \mathcal{F}_{\Gamma, K}^{+}} |\mathbf{w}_{K\sigma}| \right). \qquad (24) $$

*Proof.* Set

$$ |\mathbf{w}_K|_{\infty,K} = \max_{s \in \mathcal{V}_K} |\mathbf{w}_{\mathcal{K}s}| + \max_{\sigma \in \mathcal{F}_{\Gamma, K}^{+}} |\mathbf{w}_{K\sigma}|. $$

The definition (7) of $\nabla^K$ and the mesh regularity assumption (which ensures that $|\sigma|/|K| \lesssim h_K^{-1}$) shows that

$$ \|\nabla^K \mathbf{w}_K\|_{L^\infty(K)} \lesssim h_K^{-1} |\mathbf{w}_K|_{\infty,K} \qquad (25) $$
# 5 Proof of the error estimate

## 5.1 Inf-sup condition

We prove in this section the following result, which establishes that the bubble degree of freedom in (4) plays its role controlling the Lagrange multiplier through the jump of the displacements.

**Theorem 5.1** (Discrete inf-sup condition). It holds

$$ \sup_{\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}} \setminus \{\mathbf{0}\}} \frac{\int_{\Gamma} \lambda_{\mathcal{D}} \cdot [[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}}}{\|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}} \geq \|\lambda_{\mathcal{D}}\|_{-1/2,\Gamma} \quad \forall \lambda_{\mathcal{D}} \in \mathbf{M}_{\mathcal{D}}. \qquad (23) $$

### 5.1.1 Fracture-compatible averaged interpolator

To prove the inf-sup condition, we need a Clément-like interpolator that takes into account the fracture network. Due to the design of the $H^{-1/2}$-like norm on $M_{\mathcal{D}}$, we will use this interpolator considering only one fracture at a time. In the following, we therefore fix $i \in I$ and, in the constructions below, the sets $\mathcal{K}_s(i)$ are considered only in respect to $\Gamma_i$. Thus, for all $K \in \mathcal{M}$ and $s \in \mathcal{V}_K$, $\mathcal{K}_s(i)$ is the set of all cells $L \in \mathcal{M}_s$ that lie on the same side of $\Gamma_i$ as $K$. In particular, $\mathcal{K}_s \subset \mathcal{K}_s(i)$ with equality whenever $s$ is internal to $\Gamma_i$.

For each $K \in \mathcal{M}$ and $s \in \mathcal{V}_K$, we take an open set $U_{\mathcal{K}_s(i)} \subset \bigcup_{K \in \mathcal{K}_s(i)} K$ and a function $\varpi_{\mathcal{K}_s(i)} \in L^{\infty}(U_{\mathcal{K}_s(i)})$ such that

$$ |U_{\mathcal{K}_s(i)}| \geq \max_{K \in \mathcal{K}_s(i)} |K|, \quad |\varpi_{\mathcal{K}_s(i)}| \leq 1, \qquad (24a) $$

$$ \frac{1}{|U_{\mathcal{K}_s(i)}|} \int_{U_{\mathcal{K}_s(i)}} \varpi_{\mathcal{K}_s(i)} = 1, \quad \frac{1}{|U_{\mathcal{K}_s(i)}|} \int_{U_{\mathcal{K}_s(i)}} \mathbf{x}\varpi_{\mathcal{K}_s(i)} = \mathbf{x}_s. \qquad (24b) $$

If $s \notin \Gamma$, $U_{\mathcal{K}_s(i)}$ can be taken as a ball centered at $s$, and $\varpi_{\mathcal{K}_s(i)} = 1$. Figure 3 illustrates possible choices for $U_{\mathcal{K}_s(i)}$ depending on the nature of $s$, and Appendix A presents an explicit way to construct $(U_{\mathcal{K}_s(i)}, \varpi_{\mathcal{K}_s(i)})$ in the generic case.

<figure><img src="image_4.png" /><figcaption>Figure 3. Domains for local fracture-compatible averages.</figcaption></figure>

The space $H_0^1(\Omega \setminus \Gamma_i)$ is the subspace of $H^1(\Omega \setminus \Gamma_i)$ spanned by vector-valued functions that vanish on $\partial\Omega$, but not necessarily on $\Gamma_i$. For $\mathbf{v} \in H_0^1(\Omega \setminus \Gamma_i)$, the averaged interpolate $\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v} \in U_{0,\mathcal{D}}$ of $\mathbf{v}$ is defined by

$$ (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}_s} = \mathbf{0} \qquad \text{for all } s \in \mathcal{V}^{\text{bd}}, \qquad (25a) $$
(remember that the weights $\omega_s^\sigma$ are all nonnegative and sum up to 1). Plugged into (8) this shows that

$$ \|\Pi^K \mathbf{w}_K\|_{L^\infty(K)} \lesssim |\mathbf{w}_K|_{\infty,K}. \qquad (26) $$

The bound (24) follows by using (25) and the Cauchy–Schwarz inequality to get $\|\nabla^K \mathbf{w}_D\|_{L^2(K)} \lesssim h_K^{-1}|K|^{1/2}|\mathbf{w}_K|_{\infty,K}$, and by using (11) and (26) in the definition (13) of $\|\mathbf{w}_K\|_{1,K}$, recalling that $h_K^{d-2} \lesssim h_K^{-2}|K|$ and $\mathrm{Card}(\mathcal{V}_K) + \mathrm{Card}(\mathcal{F}_K) \lesssim 1$ by mesh regularity. $\square$

**Proposition 5.5** (Stability of the averaged interpolator). It holds

$$ \|\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}\|_{1,\mathcal{D}} \lesssim \|\nabla \mathbf{v}\|_{L^2(\Omega)} \quad \forall \mathbf{v} \in H_0^1(\Omega \setminus \Gamma_i)^d. \qquad (27) $$

*Proof.* Let $\mathbf{v} \in H_0^1(\Omega \setminus \Gamma_i)^d$. The proof is split in several steps. The idea is to first obtain a bound on a cell $K$, by locally approximating $\mathbf{v}$ by a linear function $\mathbf{q}$ and using the linear exactness of $\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a}$ resulting from the properties of $\omega_{\mathcal{K}_s}$. This allows us to bound the interpolates of $\mathbf{q}$ and of $\mathbf{v} - \mathbf{q}$, and to conclude.

**Step 1: local linear approximation $\mathbf{q}$**.

Let $K \in \mathcal{M}$ and set

$$ \mathcal{N}(K) = \bigcup_{s \in \mathcal{V}_K} \bigcup_{L \in \mathcal{K}_s(i)} L $$

the patch around $K$ made of all the cells in $\mathcal{K}_s(i)$ for each $s$ vertex of $K$. By mesh regularity, it can be checked that $\mathcal{N}(K)$ is connected by star-shaped sets as per [16, Definition 1.4], and by definition of $\mathcal{K}_s(i)$ we have $\mathbf{v} \in H^1(\mathcal{N}(K))^d$. As a consequence, the $L^2(\mathcal{N}(K))^d$-projection $\mathbf{q}$ of $\mathbf{v}$ on $\mathbb{P}^1(\mathcal{N}(K))^d$ enjoys the following approximation and stability properties:

$$ \|\mathbf{q} - \mathbf{v}\|_{L^2(\mathcal{N}(K))} \lesssim h_K \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}, \qquad (28) $$

$$ \|\mathbf{q} - \mathbf{v}\|_{L^2(\sigma)} \lesssim h_K^{1/2} \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))} \quad \forall \sigma \in \mathcal{F}_K, \qquad (29) $$

$$ \|\nabla \mathbf{q}\|_{L^2(\mathcal{N}(K))} \lesssim \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}. \qquad (30) $$

The relation (28) comes from [16, Theorem 1.45, Eq. (1.74)], accounting for the fact that the diameter of $\mathcal{N}(K)$ is $\lesssim h_K$ by mesh regularity; the trace approximation property (29) can be established following the same arguments as for Eq. (1.75) in this reference (which establishes the trace approximation property for $\sigma$ on the boundary of $\mathcal{N}(K)$); the boundedness (30) follows from [16, Remark 1.47].

**Step 2: bound on $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}$**.

In the following, by abuse of notation we denote by $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}$ the (non-averaged) local interpolate of $\mathbf{q}$ on $K$, that is, the vector obtained gathering only the degrees of freedom (9) attached to the chosen cell $K$ (for $s \in \mathcal{V}_K$ and $\sigma \in \mathcal{F}_{\Gamma,K}^+$). We also note that the DOFs lying on the boundary of $\Omega$ do not necessarily vanish for this local interpolate.

For all $s \in \mathcal{V}_K$ we have $(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}_s} = \mathbf{q}(\mathbf{x}_s)$ and thus, by choice (4) of the face weights and since $\mathbf{q}$ is linear, it holds $\nabla^{K\sigma} \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} = \nabla_\sigma \mathbf{q}$ for all $\sigma \in \mathcal{F}_K$ ($\nabla_\sigma$ being the tangential gradient), see [17, Lemma 14.8] for details. From (5), we easily infer

$$ \Pi^{K\sigma} \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} = \mathbf{q}|_{\sigma}. \qquad (31) $$
$$ (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s} = \frac{1}{|U_{\mathcal{K}s(i)}|} \int_{U_{\mathcal{K}s(i)}} \bar{\omega}_{\mathcal{K}s(i)} \mathbf{v} \quad \text{for all } K \in \mathcal{M}, s \in \mathcal{V}_K \setminus \mathcal{V}^{\text{bd}}, \qquad (25b) $$

$$ (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{K\sigma} = \frac{1}{|\sigma|} \int_{\sigma} (\gamma^{K\sigma} \mathbf{v} - \Pi^{K\sigma} (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})) \quad \text{for all } K \in \mathcal{M}, \sigma \in \mathcal{F}_{\Gamma_i, K}^{+}, \qquad (25c) $$

$$ (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{K\sigma} = \mathbf{0} \quad \text{for all } K \in \mathcal{M}, \sigma \in \mathcal{F}_{\Gamma, K}^{+} \setminus \mathcal{F}_{\Gamma_i, K}^{+}, \qquad (25d) $$

where, as expected, $\mathcal{F}_{\Gamma_i, K}^{+}$ is the set of faces on $\Gamma_i$ that belong to $\mathcal{F}_{\Gamma, K}^{+}$.

**Remark 5.2 (Averaged interpolator).** The usage of averages in the definition of the nodal values is mandatory since we need well-defined and stable interpolations of functions with only $H^1$-regularity – see Proposition 5.5. However, the averages must not be done across $\Gamma_i$ (hence the condition $U_{\mathcal{K}s(i)} \subset \cup_{K \in \mathcal{K}s(i)} K$), and must have $\mathbf{x}_s$ as center of mass (hence the integral conditions on $\bar{\omega}_{\mathcal{K}s(i)}$), to ensure that the interpolator is exact on linear functions – see the proof of Proposition 5.5.

If $s \in \Gamma$ is not in the relative interior of $\Gamma_i$, then for all $K, L \in \mathcal{M}_s$ we have $\mathcal{K}s(i) = \mathcal{L}s(i)$ and thus $(\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s} = (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{L}s}$. In other words, this $\Gamma_i$-adapted interpolator $\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a}$ produces single nodal values on $\Gamma_j$ for $j \neq i$, and possibly multiple nodal values only at vertices in the interior of $\Gamma_i$. This is coherent with the continuity properties of functions in $H_0^1(\Omega \setminus \Gamma_i)$ that it interpolates.

Finally, we note that a zero value is imposed for boundary nodes, so that $\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v} \in U_{0,\mathcal{D}}$ (computing boundary nodal values by averaging would not ensure that they vanish), and that the “bubble” value is set at zero for faces not on $\Gamma_i$.

The main properties of this interpolator are its behaviour with respect to the jump, and its $H^1$-stability.

**Proposition 5.3 (Jump of the averaged interpolate).** For all $\mathbf{v} \in H_0^1(\Omega \setminus \Gamma_i)$, it holds

$$ \forall K \in \mathcal{M}, \forall \sigma \in \mathcal{F}_K, \quad \int_{\sigma} \| \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v} \|_{\sigma} = \begin{cases} \int_{\sigma} (\gamma^{K\sigma} \mathbf{v} - \Pi^{L\sigma} \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}) & \text{if } \sigma \in \mathcal{F}_{\Gamma_i, K}^{+}, \\ \mathbf{0} & \text{if } \sigma \in \mathcal{F}_{\Gamma, K}^{+} \setminus \mathcal{F}_{\Gamma_i, K}^{+}. \end{cases} \qquad (26) $$

*Proof.* If $\sigma \in \mathcal{F}_{\Gamma_i, K}^{+}$, then the equality directly follows from the definition (7) of the jump on $\sigma$ and from (25c). If $\sigma \in \mathcal{F}_{\Gamma, K}^{+} \setminus \mathcal{F}_{\Gamma_i, K}^{+}$ then any $s \in \mathcal{V}_{\sigma}$ is not internal to $\Gamma_i$ and thus, as noted in Remark 5.2, for all $K, L \in \mathcal{M}_s$ we have $(\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s} = (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{L}s}$. Hence, $\Pi^{K\sigma} \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v} = \Pi^{L\sigma} \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}$. The definition (7) of the jump on $\sigma$, together with the fact that $(\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{K\sigma} = \mathbf{0}$ by (25d), yields the relation in (26). $\square$

Establishing the $H^1$-stability of $\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a}$ requires the following lemma.

**Lemma 5.4 (DOF-based bound on the discrete norm).** Let $K \in \mathcal{M}$. Recalling the definition (16) of $\|\cdot\|_{1,K}$, we have, for all $\mathbf{w}_K = ((\mathbf{w}_{\mathcal{K}s})_{s \in \mathcal{V}_K}, (\mathbf{w}_{K\sigma})_{\sigma \in \mathcal{F}_{\Gamma, K}^{+}})$,

$$ \|\mathbf{w}_K\|_{1,K} \lesssim h_K^{-1} |K|^{1/2} \left( \max_{s \in \mathcal{V}_K} |\mathbf{w}_{\mathcal{K}s}| + \max_{\sigma \in \mathcal{F}_{\Gamma, K}^{+}} |\mathbf{w}_{K\sigma}| \right). \qquad (27) $$

*Proof.* Set

$$ |\mathbf{w}_K|_{\infty,K} = \max_{s \in \mathcal{V}_K} |\mathbf{w}_{\mathcal{K}s}| + \max_{\sigma \in \mathcal{F}_{\Gamma, K}^{+}} |\mathbf{w}_{K\sigma}|. $$

The definition (8) of $\nabla^K$ and the mesh regularity assumption (which ensures that $|\sigma|/|K| \lesssim h_K^{-1}$) shows that

$$ \|\nabla^K \mathbf{w}_K\|_{L^\infty(K)} \lesssim h_K^{-1} |\mathbf{w}_K|_{\infty,K} \qquad (28) $$
Plugging this into the second relation of (9) yields

$$ (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{K\sigma} = 0 \quad \forall \sigma \in \mathcal{F}_{\Gamma,K}^{+}. \qquad (32) $$

Recalling (7) and using again [17, Lemma 14.8] and the fact that $\mathbf{q}$ is linear, we deduce that

$$ \nabla^K \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} = \nabla \mathbf{q}. \qquad (33) $$

The definition (8) of $\Pi^K$ then shows that $\Pi^K \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} = \mathbf{q}$. Together with (32) and the definition (11) of $S_K$, this gives

$$ S_K(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}, \mathbf{w}_{\mathcal{D}}) = 0 \quad \forall \mathbf{w}_{\mathcal{D}} \in \mathbf{U}_{\mathcal{D}}. \qquad (34) $$

Finally, recalling (33) and (30), as well as the definition (13) of $\|\cdot\|_{1,K}$, we obtain

$$ \|\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}\|_{1,K} \le \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}. \qquad (35) $$

Step 3: bound on $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}$.

We bound each type of degree of freedom (boundary vertex, internal vertex, fracture face), then use (24) to conclude.

If $s \in \mathcal{V}_K^{\text{ext}}$ then $(\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s} = 0$ and $(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}s} = \mathbf{q}(\mathbf{x}_s)$. We can then take $\sigma \in \mathcal{F}^{\text{ext}}$ that contains $s$ and write

$$ \begin{aligned} |(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}s} - (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s}| &\le |\sigma|^{-1/2} \|\mathbf{q}\|_{L^2(\sigma)} \le |\sigma|^{-1/2} h_K^{1/2} \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))} \\ &\le |K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}, \end{aligned} \qquad (36) $$

where we have used the inverse Lebesgue inequality [16, Lemma 1.25] for the first inequality, the bound (29) together with $\mathbf{v}|_{\sigma} = 0$ (since $\sigma \subset \partial\Omega$) for the second inequality, and the mesh regularity assumption to conclude.

If $s \in \mathcal{V}_K^{\text{int}}$ then (21b) and the fact that $\mathbf{q}$ is linear ensure that

$$ (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}s} = \mathbf{q}(\mathbf{x}_s) = \frac{1}{|U_{\mathcal{K}s(i)}|} \int_{U_{\mathcal{K}s(i)}} \omega_{\mathcal{K}s(i)} \mathbf{q}. $$

Hence,

$$ \begin{aligned} |(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}s} - (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s}| &= \left| \frac{1}{|U_{\mathcal{K}s(i)}|} \int_{U_{\mathcal{K}s(i)}} \omega_{\mathcal{K}s(i)} (\mathbf{q} - \mathbf{v}) \right| \\ &\le |U_{\mathcal{K}s(i)}|^{-1/2} \|\mathbf{q} - \mathbf{v}\|_{L^2(U_{\mathcal{K}s(i)})} \\ &\le |K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}, \end{aligned} \qquad (37) $$

where the first inequality follows from the bound $|\omega_{\mathcal{K}s(i)}| \le 1$ and the Cauchy-Schwarz inequality, while the conclusion is obtained using $|U_{\mathcal{K}s(i)}| \ge |K|$ (see (21a)), $U_{\mathcal{K}s(i)} \subset \mathcal{N}(K)$ and (28).

Finally, we consider the face degrees of freedom. If $\sigma \in \mathcal{F}_{\Gamma,K}^{+} \setminus \mathcal{F}_{\Gamma_i,K}^{+}$ then by (22d) and (32) we have $(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{K\sigma} = (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{K\sigma} = 0$. If $\sigma \in \mathcal{F}_{\Gamma_i,K}^{+}$ then (9) and (22c) yield

$$ (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{K\sigma} - (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{K\sigma} = \frac{1}{|\sigma|} \int_{\sigma} \gamma^{K\sigma}(\mathbf{q} - \mathbf{v}) - \Pi^{K\sigma} (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}). \qquad (38) $$
(remember that the weights $\omega_s^\sigma$ are all nonnegative and sum up to 1). Plugged into (9) this shows that

$$ \|\Pi^K \mathbf{w}_K\|_{L^\infty(K)} \lesssim |\mathbf{w}_K|_{\infty,K}. \qquad (29) $$

The bound (27) follows by using (28) and the Cauchy–Schwarz inequality to get $\|\nabla^K \mathbf{w}_D\|_{L^2(K)} \lesssim h_K^{-1}|K|^{1/2}|\mathbf{w}_K|_{\infty,K}$, and by using (12) and (29) in the definition (16) of $\|\mathbf{w}_K\|_{1,K}$, recalling that $h_K^{d-2} \lesssim h_K^{-2}|K|$ and $\mathrm{Card}(\mathcal{V}_K) + \mathrm{Card}(\mathcal{F}_K) \lesssim 1$ by mesh regularity. $\square$

**Proposition 5.5** (Stability of the averaged interpolator). It holds

$$ \|\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}\|_{1,\mathcal{D}} \lesssim \|\nabla \mathbf{v}\|_{L^2(\Omega)} \quad \forall \mathbf{v} \in \mathbf{H}_0^1(\Omega \setminus \Gamma_i). \qquad (30) $$

*Proof.* Let $\mathbf{v} \in \mathbf{H}_0^1(\Omega \setminus \Gamma_i)$. The proof is split in several steps. The idea is to first obtain a bound on a cell $K$, by locally approximating $\mathbf{v}$ by a linear function $\mathbf{q}$ and using the linear exactness of $\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a}$ resulting from the properties of $\omega_{\mathcal{K}_s}$. This allows us to bound the interpolates of $\mathbf{q}$ and of $\mathbf{v} - \mathbf{q}$, and to conclude.

*Step 1: local linear approximation q.*

Let $K \in \mathcal{M}$ and set

$$ \mathcal{N}(K) = \bigcup_{s \in \mathcal{V}_K} \bigcup_{L \in \mathcal{K}_s(i)} L $$

the patch around $K$ made of all the cells in $\mathcal{K}_s(i)$ for each $s$ vertex of $K$. By mesh regularity, it can be checked that $\mathcal{N}(K)$ is connected by star-shaped sets as per [16, Definition 1.4], and by definition of $\mathcal{K}_s(i)$ we have $\mathbf{v} \in \mathbf{H}^1(\mathcal{N}(K))$. As a consequence, the $L^2(\mathcal{N}(K))$-projection $\mathbf{q}$ of $\mathbf{v}$ on $\mathbb{P}^1(\mathcal{N}(K))^d$ enjoys the following approximation and stability properties:

$$ \|\mathbf{q} - \mathbf{v}\|_{L^2(\mathcal{N}(K))} \lesssim h_K \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}, \qquad (31) $$

$$ \|\mathbf{q} - \mathbf{v}\|_{L^2(\sigma)} \lesssim h_K^{1/2} \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))} \quad \forall \sigma \in \mathcal{F}_K, \qquad (32) $$

$$ \|\nabla \mathbf{q}\|_{L^2(\mathcal{N}(K))} \lesssim \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}. \qquad (33) $$

The relation (31) comes from [16, Theorem 1.45, Eq. (1.74)], accounting for the fact that the diameter of $\mathcal{N}(K)$ is $\lesssim h_K$ by mesh regularity; the trace approximation property (32) can be established following the same arguments as for Eq. (1.75) in this reference (which establishes the trace approximation property for $\sigma$ on the boundary of $\mathcal{N}(K)$); the boundedness (33) follows from [16, Remark 1.47].

*Step 2: bound on $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}$.*

In the following, by abuse of notation we denote by $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}$ the (non-averaged) local interpolate of $\mathbf{q}$ on $K$, that is, the vector obtained gathering only the degrees of freedom (10) attached to the chosen cell $K$ (for $s \in \mathcal{V}_K$ and $\sigma \in \mathcal{F}_{\Gamma,K}^+$). We also note that the DOFs lying on the boundary of $\Omega$ do not necessarily vanish for this local interpolate.

For all $s \in \mathcal{V}_K$ we have $(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}_s} = \mathbf{q}(\mathbf{x}_s)$ and thus, by choice (5) of the face weights and since $\mathbf{q}$ is linear, it holds $\nabla^{K\sigma} \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} = \nabla_\sigma \mathbf{q}$ for all $\sigma \in \mathcal{F}_K$ ($\nabla_\sigma$ being the tangential gradient), see [17, Lemma 14.8] for details. From (6), we easily infer

$$ \Pi^{K\sigma} \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} = \mathbf{q}|_{\sigma}. \qquad (34) $$

Plugging this into the second relation of (10) yields

$$ (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{K\sigma} = \mathbf{0} \quad \forall \sigma \in \mathcal{F}_{\Gamma,K}^{+}. \qquad (35) $$
By (36) and (37), all vertex degrees of freedom of $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}$ are $O(|K|^{-1/2}h_K \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))})$. The operator $\Pi^{K\sigma}$ only depends on these degrees of freedom so, by the same arguments that led to (26), we find that $\|\Pi^{K\sigma}(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}|_{\sigma} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})\|_{L^{\infty}(\sigma)} \lesssim |K|^{-1/2}h_K \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}$. Plugging this into (38) and using the Cauchy–Schwarz inequality, we infer that

$$ \begin{aligned} |(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{K\sigma} - (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{K\sigma}| & \lesssim |\sigma|^{-1/2} \|\mathbf{q} - \mathbf{v}\|_{L^2(\sigma)} + |K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))} \\ & \lesssim |K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}, \end{aligned} \quad (39) $$

where the mesh regularity and (29) were used to conclude.

Gathering (36), (37) and (39) and invoking the bound (24) with $\mathbf{w}_K$ gathering the degrees of freedom on $K$ of $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}$ yields

$$ \|\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}\|_{1,K} \lesssim \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}. \quad (40) $$

Step 4: conclusion.

Combining (35) and (40), squaring, summing over $K \in \mathcal{M}$ and gathering the integrals by cells we infer

$$ \|\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}\|_{1,\mathcal{D}}^2 \lesssim \sum_{K \in \mathcal{M}} \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}^2 = \sum_{K' \in \mathcal{M}} \|\nabla \mathbf{v}\|_{L^2(K')}^2 \mathrm{Card}\{K \in \mathcal{M} : K' \in \mathcal{N}(K)\}. $$

The proof of (27) is complete by using the mesh regularity properties to see that $\mathrm{Card}\{K \in \mathcal{M} : K' \in \mathcal{N}(K)\} \lesssim 1$. $\square$

### 5.1.2 Proof of the inf-sup condition

We are now ready to prove Theorem 5.1. Let $\lambda_{\mathcal{D}} \in \mathcal{M}_{\mathcal{D}}$. Recalling the definition (14) of the $H_{0,j}^{-1/2}(\Gamma)$-like norm, the inf-sup condition (20) holds provided we can prove that, for all $i \in I$ and all $\mathbf{v}_i \in H_0^1(\Omega_i^+; \Gamma_i)^d$, there exists $\mathbf{v}_{\mathcal{D}} \in U_{0,\mathcal{D}}$ such that

$$ \frac{\int_{\Gamma} \lambda_{\mathcal{D}} \cdot [[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}}}{\|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}} \geq \frac{\int_{\Gamma_i} \lambda_{\mathcal{D}} \cdot \mathbf{v}_i}{\|\mathbf{v}_i\|_{H^1(\Omega_i^+)}}, \quad (41) $$

Since $\mathbf{v}_i$ vanishes on $\partial \Omega_i^+ \setminus \Gamma_i$, its extension $\tilde{\mathbf{v}}_i$ by 0 to $\Omega$ belongs to $H_0^1(\Omega \setminus \Gamma_i)^d$ and satisfies $\|\tilde{\mathbf{v}}_i\|_{H^1(\Omega \setminus \Gamma_i)} = \|\mathbf{v}_i\|_{H_0^1(\Omega_i^+)}$. Let us consider $\mathbf{v}_{\mathcal{D}} = \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i$ with the interpolator defined by (22); by the construction done in Appendix A, we can define this interpolator in such a way that, for all $s \in \bar{\Gamma}_i$ and $L$ on the negative side of $\Gamma_i$, $U_{\mathcal{L}s(i)}$ is fully contained in $\Omega \setminus \Omega_i^+$ (since $U_{\mathcal{L}s(i)}$ can be fully contained in any chosen cell that contains $s$). This ensures that $(\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i)_{\mathcal{L}s(i)} = 0$ for all such $s$, and thus that $\Pi^{L\sigma} \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i = 0$ for all such $L$. Hence, by (23) and since $\lambda_{\mathcal{D}}$ is piecewise constant on $\mathcal{F}_{\Gamma}$,

$$ \int_{\Gamma} \lambda_{\mathcal{D}} \cdot [[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}} = \sum_{\sigma \in \mathcal{F}_{\Gamma}} \lambda_{\sigma} \cdot \int_{\sigma} [[\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i]]_{\sigma} = \sum_{\sigma \in \mathcal{F}_{\Gamma_i}} \lambda_{\sigma} \cdot \int_{\sigma} \gamma^{K\sigma} \tilde{\mathbf{v}}_i = \int_{\Gamma_i} \lambda_{\mathcal{D}} \cdot \mathbf{v}_i, $$

where the conclusion follows from the fact that $\gamma^{K\sigma}$ is the trace on the positive side, where $\tilde{\mathbf{v}}_i = \mathbf{v}_i$. Dividing throughout by $\|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}} = \|\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i\|_{1,\mathcal{D}} \lesssim \|\tilde{\mathbf{v}}_i\|_{H^1(\Omega \setminus \Gamma_i)} = \|\mathbf{v}_i\|_{H_0^1(\Omega_i^+)}$ (see (27)) we infer that (41) holds.
Recalling (8) and using again [17, Lemma 14.8] and the fact that $\mathbf{q}$ is linear, we deduce that

$$ \nabla^K \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} = \nabla \mathbf{q}. \qquad (36) $$

The definition (9) of $\Pi^K$ then shows that $\Pi^K \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} = \mathbf{q}$. Together with (35) and the definition (12) of $S_K$, this gives

$$ S_K(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}, \mathbf{w}_{\mathcal{D}}) = 0 \quad \forall \mathbf{w}_{\mathcal{D}} \in \mathbf{U}_{\mathcal{D}}. \qquad (37) $$

Finally, recalling (36) and (33), as well as the definition (16) of $\|\cdot\|_{1,K}$, we obtain

$$ \| \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} \|_{1,K} \lesssim \| \nabla \mathbf{v} \|_{L^2(N(K))}. \qquad (38) $$

Step 3: bound on $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}$.

We bound each type of degree of freedom (boundary vertex, internal vertex, fracture face), then use (27) to conclude.

If $s \in \mathcal{V}_K^{\text{bd}}$ then $(\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s} = \mathbf{0}$ and $(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}s} = \mathbf{q}(\mathbf{x}_s)$. We can then take $\sigma \in \mathcal{F}^{\text{bd}}$ that contains $s$ and write

$$ \begin{aligned} |(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}s} - (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s}| &\lesssim |\sigma|^{-1/2} \|\mathbf{q}\|_{L^2(\sigma)} \lesssim |\sigma|^{-1/2} h_K^{1/2} \|\nabla \mathbf{v}\|_{L^2(N(K))} \\ &\lesssim |K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(N(K))}, \end{aligned} \qquad (39) $$

where we have used the inverse Lebesgue inequality [16, Lemma 1.25] for the first inequality, the bound (32) together with $\mathbf{v}|_{\sigma} = \mathbf{0}$ (since $\sigma \subset \partial\Omega$) for the second inequality, and the mesh regularity assumption to conclude.

If $s \in \mathcal{V}_K^{\text{in}}$ then (24b) and the fact that $\mathbf{q}$ is linear ensure that

$$ (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}s} = \mathbf{q}(\mathbf{x}_s) = \frac{1}{|U_{\mathcal{K}s(i)}|} \int_{U_{\mathcal{K}s(i)}} \overline{\omega_{\mathcal{K}s(i)} \mathbf{q}}. $$

Hence,

$$ \begin{aligned} |(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}s} - (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}s}| &= \left| \frac{1}{|U_{\mathcal{K}s(i)}|} \int_{U_{\mathcal{K}s(i)}} \overline{\omega_{\mathcal{K}s(i)}(\mathbf{q}-\mathbf{v})} \right| \\ &\lesssim |U_{\mathcal{K}s(i)}|^{-1/2} \|\mathbf{q}-\mathbf{v}\|_{L^2(U_{\mathcal{K}s(i)})} \\ &\lesssim |K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(N(K))}, \end{aligned} \qquad (40) $$

where the first inequality follows from the bound $|\overline{\omega_{\mathcal{K}s(i)}}| \lesssim 1$ and the Cauchy–Schwarz inequality, while the conclusion is obtained using $|U_{\mathcal{K}s(i)}| \gtrsim |K|$ (see (24a)), $U_{\mathcal{K}s(i)} \subset \mathcal{N}(K)$ and (31).

Finally, we consider the face degrees of freedom. If $\sigma \in \mathcal{F}_{\Gamma, K}^{+} \setminus \mathcal{F}_{\Gamma_i, K}^{+}$ then by (25d) and (35) we have $(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}\sigma} = (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}\sigma} = \mathbf{0}$. If $\sigma \in \mathcal{F}_{\Gamma_i, K}^{+}$ then (10) and (25c) yield

$$ (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}\sigma} - (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}\sigma} = \frac{1}{|\sigma|} \int_{\sigma} \gamma^{\mathcal{K}\sigma}(\mathbf{q}-\mathbf{v}) - \Pi^{\mathcal{K}\sigma} (\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}). \qquad (41) $$

By (39) and (40), all vertex degrees of freedom of $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}$ are $O(|K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(N(K))})$. The operator $\Pi^{\mathcal{K}\sigma}$ only depends on these degrees of freedom so, by the same arguments that led to (29), we find that $\|\Pi^{\mathcal{K}\sigma}(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}|_{\sigma} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})\|_{L^{\infty}(\sigma)} \lesssim |K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(N(K))}$. Plugging this into (41) and using the Cauchy–Schwarz inequality, we infer that

$$ |(\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q})_{\mathcal{K}\sigma} - (\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v})_{\mathcal{K}\sigma}| \lesssim |\sigma|^{-1/2} \|\mathbf{q} - \mathbf{v}\|_{L^2(\sigma)} + |K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(N(K))} $$
Remark 5.6 (Jump-based norm on the fracture network). An apparent simpler – and perhaps more natural – choice, that would prevent us from having to cut the fracture network into planar components $(\Gamma_i)_{i \in I}$, would be to define the $H_{0,i}^{1/2}(\Omega \setminus \Gamma)$-norm using test functions $\mathbf{v} \in H_0^1(\Omega \setminus \Gamma)^d$ and replacing their trace on the positive side of the fractures by their jump across the fractures. This choice is fully valid if we consider discretisations that have two bubble degrees of freedom, one on each side of the fracture, because in that case we have

$$ \int_{\sigma} \left[ \left\| \mathcal{I}_{U_{0,\mathcal{D}}}^{a} \mathbf{v} \right\|_{\sigma} \right] = \int_{\sigma} \left[ \left\| \mathbf{v} \right\| \right] \quad \forall \mathbf{v} \in H_{0}^{1}(\Omega \setminus \Gamma)^{d}. $$

However, in the case of a single bubble degree of freedom located only on one side of the fracture (which has some practical interest in case of non-matching interfaces), the best relation that can be established between $\left[ \left\| \mathcal{I}_{U_{0,\mathcal{D}}}^{a} \mathbf{v} \right\|_{\sigma} \right]$ and $\mathbf{v}$ is (23) (dropping the index $i$ since we are discussing the option to handle the whole network at once), in which $\Pi^{L\sigma} \mathcal{I}_{U_{0,\mathcal{D}}}^{a} \mathbf{v}$ is not the trace of $\mathbf{v}$ on the negative side, but some approximation thereof. As a consequence, when bubble degrees of freedom are only taken on the positive side of the network, we have to ensure that $\Pi^{L\sigma} \mathcal{I}_{U_{0,\mathcal{D}}}^{a} \mathbf{v} = 0$ on the negative side; this is achieved in the proof above by working on each individual planar component of the network and using the zero extension outside $\Omega_i^+$ (which provides zero values on the negative side of $\Gamma$). For simple networks we might be able to create a $\mathbf{v}$ that is zero on the whole negative side, but for complex network with multiple intersection lines, this does not look feasible, and explains why we resort to the decomposed norm (14).

## 5.2 Discrete Korn inequality

**Theorem 5.7 (Discrete Korn inequality).** It holds

$$ \|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}^2 \lesssim \|\varepsilon_{\mathcal{D}}(\mathbf{v}_{\mathcal{D}})\|_{L^2(\Omega \setminus \Gamma)}^2 + S_{\mathcal{D}}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) \quad \forall \mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}. \quad (42) $$

*Proof.* We begin by substituting the quantities $\nabla^{\mathcal{D}}\mathbf{v}_{\mathcal{D}}$ in $\|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}$ and $\varepsilon_{\mathcal{D}}(\mathbf{v}_{\mathcal{D}})$ with their equivalents $\nabla_{\mathcal{M}}(\Pi^{\mathcal{D}}\mathbf{v}_{\mathcal{D}})$ and $\varepsilon_{\mathcal{M}}(\Pi^{\mathcal{D}}\mathbf{v}_{\mathcal{D}})$ (where $\nabla_{\mathcal{M}}$ and $\varepsilon_{\mathcal{M}}$ are the broken gradient and strain on $\mathcal{M}$) to expose the function $\Pi^{\mathcal{D}}\mathbf{v}_{\mathcal{D}} \in \mathbb{P}^1(\mathcal{M})^d$, and enable the usage of [16, Lemma 7.23] with adjustments to account for fractures.

Specifically, the node-averaging operator $\mathcal{I}_{\text{av},h}^1$ from [16, Section 7.3.2] can be adapted to avoid crossing fractures, by defining the node-averaged value at a given $s$ on the side $K$ of the fracture as follows

$$ (\mathcal{I}_{\text{av},h}^1 \Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}})_{\mathcal{K}_s} = \frac{1}{\#\mathcal{K}_s} \sum_{L \in \mathcal{K}_s} \Pi^L \mathbf{v}_{\mathcal{D}}(\mathbf{x}_s). $$

It results that the jumps $(\mathcal{I}_{\text{av},h}^1 \Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}})_{\mathcal{K}_s} - \Pi^L \mathbf{v}_{\mathcal{D}}(\mathbf{x}_s)$ can be controlled using only non-fracture faces, instead of all mesh faces as in [16, Section 7.3.2]. Applying the techniques there, we obtain the following modified version of [16, Eq. (7.66)]:

$$ \|\nabla_{\mathcal{M}}(\Pi^{\mathcal{D}}\mathbf{v}_{\mathcal{D}})\|_{L^2(\Omega \setminus \Gamma)}^2 \lesssim \|\varepsilon_{\mathcal{M}}(\Pi^{\mathcal{D}}\mathbf{v}_{\mathcal{D}})\|_{L^2(\Omega \setminus \Gamma)}^2 + \sum_{\sigma \in \mathcal{F} \setminus \mathcal{F}_{\Gamma}} h_{\sigma}^{-1} \|\left[\Pi^{\mathcal{D}}\mathbf{v}_{\mathcal{D}}\right]_{|\sigma}\|_{L^2(\sigma)}^2. \quad (43) $$

To conclude, we need to bound the jump terms in the right-hand side. We consider a face $K|L = \sigma \in \mathcal{F}^{\text{int}} \setminus \mathcal{F}_{\Gamma}$. For $s \in \mathcal{V}_{\sigma}$, noticing that $\mathbf{v}_{\mathcal{K}_s} = \mathbf{v}_{\mathcal{L}_s}$, we write

$$ \begin{aligned} \left|\left[\Pi^{\mathcal{D}}\mathbf{v}_{\mathcal{D}}\right]_{|\sigma}(\mathbf{x}_s)\right|^2 &= |(\Pi^K \mathbf{v}_{\mathcal{D}})(\mathbf{x}_s) - (\Pi^L \mathbf{v}_{\mathcal{D}})(\mathbf{x}_s)|^2 \\ &\le 2|\mathbf{v}_{\mathcal{K}_s} - (\Pi^K \mathbf{v}_{\mathcal{D}})(\mathbf{x}_s)|^2 + 2|\mathbf{v}_{\mathcal{L}_s} - (\Pi^L \mathbf{v}_{\mathcal{D}})(\mathbf{x}_s)|^2 \end{aligned} $$
$$ \lesssim |K|^{-1/2} h_K \|\nabla \mathbf{v}\|_{L^2(\mathcal{N}(K))}, \quad (42) $$

where the mesh regularity and (32) were used to conclude.

Gathering (39), (40) and (42) and invoking the bound (27) with $\mathbf{w}_K$ gathering the degrees of freedom on $K$ of $\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v}$ yields

$$ \| \mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q} - \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v} \|_{1,K} \lesssim \| \nabla \mathbf{v} \|_{L^2(\mathcal{N}(K))}. \quad (43) $$

Step 4: conclusion.

Combining (38) and (43), squaring, summing over $K \in \mathcal{M}$ and gathering the integrals by cells we infer

$$ \| \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \mathbf{v} \|_{1,\mathcal{D}}^2 \lesssim \sum_{K \in \mathcal{M}} \| \nabla \mathbf{v} \|_{L^2(\mathcal{N}(K))}^2 = \sum_{K' \in \mathcal{M}} \| \nabla \mathbf{v} \|_{L^2(K')}^2 \mathrm{Card}\{K \in \mathcal{M} : K' \in \mathcal{N}(K)\}. $$

The proof of (30) is complete by using the mesh regularity properties to see that $\mathrm{Card}\{K \in \mathcal{M} : K' \in \mathcal{N}(K)\} \lesssim 1$. $\square$

### 5.1.2 Proof of the inf-sup condition

We are now ready to prove Theorem 5.1. Let $\lambda_{\mathcal{D}} \in \mathcal{M}_{\mathcal{D}}$. Recalling the definition (17) of the $H^{-1/2}$-like norm, the inf-sup condition (23) holds provided we can prove that, for all $i \in I$ and all $\mathbf{v}_i \in H_0^1(\Omega_i^+; \Gamma_i)$, there exists $\mathbf{v}_{\mathcal{D}} \in U_{0,\mathcal{D}}$ such that

$$ \frac{\int_{\Gamma} \lambda_{\mathcal{D}} \cdot [[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}}}{\|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}} \geq \frac{\int_{\Gamma_i} \lambda_{\mathcal{D}} \cdot \mathbf{v}_i}{\|\mathbf{v}_i\|_{H^1(\Omega_i^+)}}, \quad (44) $$

Since $\mathbf{v}_i$ vanishes on $\partial \Omega_i^+ \setminus \Gamma_i$, its extension $\tilde{\mathbf{v}}_i$ by $\mathbf{0}$ to $\Omega$ belongs to $H_0^1(\Omega \setminus \Gamma_i)$ and satisfies $\|\tilde{\mathbf{v}}_i\|_{H^1(\Omega \setminus \Gamma_i)} = \|\mathbf{v}_i\|_{H_0^1(\Omega_i^+)}$. Let us consider $\mathbf{v}_{\mathcal{D}} = \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i$ with the interpolator defined by (25); by the construction done in Appendix A, we can define this interpolator in such a way that, for all $s \in \bar{\Gamma}_i$ and $L$ on the negative side of $\Gamma_i$, $U_{\mathcal{L}_s(i)}$ is fully contained in $\Omega \setminus \Omega_i^+$ (since $U_{\mathcal{L}_s(i)}$ can be fully contained in any chosen cell that contains $s$). This ensures that $(\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i)_{\mathcal{L}_s(i)} = \mathbf{0}$ for all such $s$, and thus that $\Pi^{\mathcal{L}\sigma} \mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i = \mathbf{0}$ for all such $L$. Hence, by (26) and since $\lambda_{\mathcal{D}}$ is piecewise constant on $\mathcal{F}_{\Gamma}$,

$$ \int_{\Gamma} \lambda_{\mathcal{D}} \cdot [[\mathbf{v}_{\mathcal{D}}]]_{\mathcal{D}} = \sum_{\sigma \in \mathcal{F}_{\Gamma}} \lambda_{\sigma} \cdot \int_{\sigma} [[\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i]]_{\sigma} = \sum_{\sigma \in \mathcal{F}_{\Gamma_i}} \lambda_{\sigma} \cdot \int_{\sigma} \gamma^{\mathcal{K}\sigma} \tilde{\mathbf{v}}_i = \int_{\Gamma_i} \lambda_{\mathcal{D}} \cdot \mathbf{v}_i, $$

where the conclusion follows from the fact that $\gamma^{\mathcal{K}\sigma}$ is the trace on the positive side, where $\tilde{\mathbf{v}}_i = \mathbf{v}_i$. Dividing throughout by $\|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}} = \|\mathcal{I}_{U_{0,\mathcal{D}}}^{i,a} \tilde{\mathbf{v}}_i\|_{1,\mathcal{D}} \lesssim \|\tilde{\mathbf{v}}_i\|_{H^1(\Omega \setminus \Gamma_i)} = \|\mathbf{v}_i\|_{H^1(\Omega_i^+)}$ (see (30)) we infer that (44) holds.

**Remark 5.6** (Jump-based norm on the fracture network). An apparent simpler – and perhaps more natural – choice, that would prevent us from having to cut the fracture network into planar components $(\Gamma_i)_{i \in I}$, would be to define the $H^{-1/2}$-like norm on $\mathcal{M}_{\mathcal{D}}$ using test functions $\mathbf{v} \in U_0$ and to replace their trace on the positive side of the fractures by their jump across the fractures. This choice is fully valid if we consider discretisations that have two bubble degrees of freedom, one on each side of the fracture, because in that case we have

$$ \int_{\sigma} [[\mathcal{I}_{U_{0,\mathcal{D}}}^{a} \mathbf{v}]]_{\sigma} = \int_{\sigma} [[\mathbf{v}]] \quad \forall \mathbf{v} \in U_0. $$
$$ \lesssim h_K^{2-d} (S_K(\mathbf{v}_D, \mathbf{v}_D) + S_L(\mathbf{v}_D, \mathbf{v}_D)), \quad (44) $$

where the conclusion follows from the definition (11) of the local stabilisation terms and the mesh regularity (which ensures $h_K \approx h_L$). Further, the fact that $\|\Pi^D \mathbf{v}_D\|_{\sigma} \in \mathbb{P}^1(\sigma)^d$ along with (44) yields

$$ \begin{aligned} \|\|\Pi^D \mathbf{v}_D\|\|_{L^2(\sigma)}\|_{\sigma}^2 &\le |\sigma| \|\|\Pi^D \mathbf{v}_D\|\|_{\sigma}\|_{\infty}^2 \lesssim h_K^{d-1} \max_{s \in \mathcal{V}_{\sigma}} \|\|\Pi^D \mathbf{v}_D\|\|_{\sigma}(\mathbf{x}_s)\|^2 \\ &\lesssim h_K (S_K(\mathbf{v}_D, \mathbf{v}_D) + S_L(\mathbf{v}_D, \mathbf{v}_D)). \end{aligned} \quad (45) $$

In the second inequality, we have used the fact that, by mesh regularity property, any $\mathbf{x} \in \sigma$ can be written $\mathbf{x} = \sum_{s \in \mathcal{V}_{\sigma}} \rho_s \mathbf{x}_s$ with $\sum_{s \in \mathcal{V}_{\sigma}} \rho_s = 1$ and $\sum_{s \in \mathcal{V}_{\sigma}} |\rho_s| \le 1$, and that $\|\Pi^D \mathbf{v}_D\|(\mathbf{x}) = \sum_{s \in \mathcal{V}_{\sigma}} \rho_s \|\Pi^D \mathbf{v}_D\|(\mathbf{x}_s)$. We deduce from (45) that

$$ h_{\sigma}^{-1} \|\|\Pi^D \mathbf{v}_D\|\|_{\sigma} \|_{L^2(\sigma)}^2 \lesssim S_K(\mathbf{v}_D, \mathbf{v}_D) + S_L(\mathbf{v}_D, \mathbf{v}_D). \quad (46) $$

Consider now a face $\sigma \in \mathcal{F}^{\text{ext}} \cap \mathcal{F}_K$, $K \in \mathcal{M}$. For $s \in \mathcal{V}_{\sigma}$, noticing that $\mathbf{v}_{\mathcal{K}s} = 0$, we obtain

$$ \|\|\Pi^D \mathbf{v}_D\|\|_{\sigma}(\mathbf{x}_s)\|^2 = |(\Pi^K \mathbf{v}_D)(\mathbf{x}_s) - \mathbf{v}_{\mathcal{K}s}|^2 \le h_K^{2-d} S_K(\mathbf{v}_D, \mathbf{v}_D). $$

Then, following the same arguments that lead to (46), we get

$$ h_{\sigma}^{-1} \|\|\Pi^D \mathbf{v}_D\|\|_{\sigma} \|_{L^2(\sigma)}^2 \lesssim S_K(\mathbf{v}_D, \mathbf{v}_D). \quad (47) $$

Summing (46) over the faces $\sigma \in \mathcal{F}^{\text{int}} \setminus \mathcal{F}_{\Gamma}$ and (47) over the faces $\sigma \in \mathcal{F}^{\text{ext}}$, we deduce

$$ \sum_{\sigma \in \mathcal{F} \setminus \mathcal{F}_{\Gamma}} h_{\sigma}^{-1} \|\|\Pi^D \mathbf{v}_D\|\|_{\sigma} \|_{L^2(\sigma)}^2 \lesssim \sum_{K \in \mathcal{M}} S_K(\mathbf{v}_D, \mathbf{v}_D). \quad (48) $$

The proof of (42) is concluded by noticing that

$$ \|\mathbf{v}_D\|_{1,\mathcal{D}}^2 = \|\nabla_{\mathcal{M}}(\Pi^D \mathbf{v}_D)\|_{L^2(\Omega \setminus \bar{\Gamma})}^2 + \sum_{K \in \mathcal{M}} S_K(\mathbf{v}_D, \mathbf{v}_D), $$

and by using (43) and (48). □

## 5.3 Proof of the abstract error estimate

We prove here Theorem 4.5. To shorten the notations, let us define the discrete energy inner product $\langle \cdot, \cdot \rangle_{e,\mathcal{D}}$ such that, for $\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{\mathcal{D}}$,

$$ \langle \mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}} \rangle_{e, \mathcal{D}} = \int_{\Omega} \varphi_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}) : \varphi_{\mathcal{D}}(\mathbf{v}_{\mathcal{D}}) + S_{\mu, \lambda, \mathcal{D}}(\mathbf{u}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}), $$

and denote by $\|\cdot\|_{e,\mathcal{D}}$ its associated norm. From the discrete Korn inequality (42), we deduce the following bound for all $\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}$:

$$ \|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}} \lesssim \|\mathbf{v}_{\mathcal{D}}\|_{e,\mathcal{D}}. \quad (49) $$

Since $\mathbf{u}$ is a weak solution to the contact problem and by the regularity assumed on $\lambda$ in the theorem, we have $-\text{div}(\sigma(\mathbf{u})) = \mathbf{f} \in L^2(\Omega)$ and $\lambda = -\gamma_{\mathbf{n}}^{+} \sigma(\mathbf{u}) \in L^2(\Gamma)$. Hence, $\sigma(\mathbf{u}) \in \mathbf{W}$ defined by (16) and, using the definition (17) of $w_{\mathcal{D}}$, we have for all $\mathbf{w}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}$

$$ \int_{\Omega} \varphi(\mathbf{u}) : \varphi_{\mathcal{D}}(\mathbf{w}_{\mathcal{D}}) - \int_{\Omega} \mathbf{f} \cdot \tilde{\Pi}^{\mathcal{D}} \mathbf{w}_{\mathcal{D}} = -w_{\mathcal{D}}(\sigma(\mathbf{u}), \mathbf{w}_{\mathcal{D}}) - \int_{\Gamma} \lambda \cdot \|\mathbf{w}_{\mathcal{D}}\|_{\mathcal{D}}. \quad (50) $$
However, in the case of a single bubble degree of freedom located only on one side of the fracture (which has some practical interest in case of non-matching interfaces), the best relation that can be established between $\llbracket I_{U_{0,D}}^a \mathbf{v} \rrbracket_\sigma$ and $\mathbf{v}$ is (26) (dropping the index $i$ since we are discussing the option to handle the whole network at once), in which $\Pi^{L\sigma} I_{U_{0,D}}^a \mathbf{v}$ is not the trace of $\mathbf{v}$ on the negative side, but some approximation thereof. As a consequence, when bubble degrees of freedom are only taken on the positive side of the network, we have to ensure that $\Pi^{L\sigma} I_{U_{0,D}}^a \mathbf{v} = \mathbf{0}$ on the negative side; this is achieved in the proof above by working on each individual planar component of the network and using the zero extension outside $\Omega_i^+$ (which provides zero values on the negative side of $\Gamma_i$). For simple networks we might be able to create a $\mathbf{v}$ that is zero on the whole negative side, but for complex network with multiple intersection lines, this does not look feasible, and explains why we resort to the decomposed norm (17).

## 5.2 Discrete Korn inequality

**Theorem 5.7 (Discrete Korn inequality).** It holds

$$ \|\mathbf{v}_D\|_{1,D}^2 \lesssim \|\mathfrak{e}^{\mathcal{D}}(\mathbf{v}_D)\|_{L^2(\Omega\setminus\Gamma)}^2 + S_D(\mathbf{v}_D, \mathbf{v}_D) \quad \forall \mathbf{v}_D \in \mathbf{U}_{0,D}. \quad (45) $$

*Proof.* We begin by substituting the quantities $\nabla^{\mathcal{D}}\mathbf{v}_D$ in $\|\mathbf{v}_D\|_{1,D}$ and $\mathfrak{e}^{\mathcal{D}}(\mathbf{v}_D)$ with their equivalents $\nabla_M(\Pi^{\mathcal{D}}\mathbf{v}_D)$ and $\mathfrak{e}^M(\Pi^{\mathcal{D}}\mathbf{v}_D)$ (where $\nabla_M$ and $\mathfrak{e}^M$ are the broken gradient and strain on $M$) to expose the function $\Pi^{\mathcal{D}}\mathbf{v}_D \in \mathbb{P}^1(M)^d$, and enable the usage of [16, Lemma 7.23] with adjustments to account for fractures.

Specifically, the node-averaging operator $I_{av,h}^1$ from [16, Section 7.3.2] can be adapted to avoid crossing fractures, by defining the node-averaged value at a given $s$ on the side $K$ of the fracture as follows

$$ (I_{av,h}^1 \Pi^{\mathcal{D}} \mathbf{v}_D)_{\mathcal{K}_s} = \frac{1}{\#\mathcal{K}_s} \sum_{L \in \mathcal{K}_s} \Pi^L \mathbf{v}_D(\mathbf{x}_s). $$

It results that, for all $L \in \mathcal{K}_s$, the jump $(I_{av,h}^1 \Pi^{\mathcal{D}} \mathbf{v}_D)_{\mathcal{K}_s} - \Pi^L \mathbf{v}_D(\mathbf{x}_s)$ can be controlled using only non-fracture faces, instead of all mesh faces as in [16, Section 7.3.2]. Applying the techniques there, we obtain the following modified version of [16, Eq. (7.66)]:

$$ \|\nabla_M(\Pi^{\mathcal{D}}\mathbf{v}_D)\|_{L^2(\Omega\setminus\Gamma)}^2 \lesssim \|\mathfrak{e}^M(\Pi^{\mathcal{D}}\mathbf{v}_D)\|_{L^2(\Omega\setminus\Gamma)}^2 + \sum_{\sigma \in \mathcal{F}\setminus\mathcal{F}_\Gamma} h_\sigma^{-1} \|\llbracket\Pi^{\mathcal{D}}\mathbf{v}_D\rrbracket\|_{\sigma} \|\Pi^L(\mathbf{v}_D)\|_{L^2(\sigma)}^2. \quad (46) $$

We now bound the jump terms above. We consider a face $K|L = \sigma \in \mathcal{F}^{\text{in}} \setminus \mathcal{F}_\Gamma$. For $s \in \mathcal{V}_\sigma$, noticing that $\mathbf{v}_{\mathcal{K}_s} = \mathbf{v}_{\mathcal{L}_s}$, we write

$$ \begin{aligned} |\llbracket \Pi^{\mathcal{D}} \mathbf{v}_D \rrbracket|_{\sigma}(\mathbf{x}_s)|^2 &= |(\Pi^K \mathbf{v}_D)(\mathbf{x}_s) - (\Pi^L \mathbf{v}_D)(\mathbf{x}_s)|^2 \\ &\le 2|\mathbf{v}_{\mathcal{K}_s} - (\Pi^K \mathbf{v}_D)(\mathbf{x}_s)|^2 + 2|\mathbf{v}_{\mathcal{L}_s} - (\Pi^L \mathbf{v}_D)(\mathbf{x}_s)|^2 \\ &\lesssim h_K^{2-d} (S_K(\mathbf{v}_D, \mathbf{v}_D) + S_L(\mathbf{v}_D, \mathbf{v}_D)), \end{aligned} \quad (47) $$

where the conclusion follows from the definition (12) of the local stabilisation terms and the mesh regularity (which ensures $h_K \approx h_L$). Further, the fact that $\llbracket \Pi^{\mathcal{D}} \mathbf{v}_D \rrbracket|_{\sigma} \in \mathbb{P}^1(\sigma)^d$ along with (47) yields

$$ \begin{aligned} \|\llbracket \Pi^{\mathcal{D}} \mathbf{v}_D \rrbracket|_{\sigma} \|_{L^2(\sigma)}^2 &\le |\sigma| \|\llbracket \Pi^{\mathcal{D}} \mathbf{v}_D \rrbracket|_{\sigma} \|_{\infty}^2 \lesssim h_K^{d-1} \max_{s \in \mathcal{V}_\sigma} |\llbracket \Pi^{\mathcal{D}} \mathbf{v}_D \rrbracket|_{\sigma}(\mathbf{x}_s)|^2 \\ &\lesssim h_K(S_K(\mathbf{v}_D, \mathbf{v}_D) + S_L(\mathbf{v}_D, \mathbf{v}_D)). \end{aligned} \quad (48) $$
Subtracting (10a) (with $v_D = w_D$) from (50), we obtain

$$ \int_{\Omega} (\sigma(\mathbf{u}) - \sigma_D(\mathbf{u}_D)) : \mathcal{E}_D(\mathbf{w}_D) - S_{\mu, \lambda, D}(\mathbf{u}_D, \mathbf{w}_D) = -w_D(\sigma(\mathbf{u}), \mathbf{w}_D) + \int_{\Gamma} (\lambda_D - \lambda) \cdot [[\mathbf{w}_D]]_D. \quad (51) $$

Take $v_D \in \mathbf{U}_{0,D}$ and set $w_D = v_D - u_D$ in (51) to get

$$ \begin{aligned} \|v_D - u_D\|_{e, D}^2 = & -w_D(\sigma(\mathbf{u}), v_D - u_D) - \int_{\Omega} (\sigma(\mathbf{u}) - \sigma_D(v_D)) : \mathcal{E}_D(v_D - u_D) \\ & + S_{\mu, \lambda, D}(v_D, v_D - u_D) + \int_{\Gamma} (\lambda_D - \lambda) \cdot ([[v_D]]_D - [[\mathbf{u}]]) \\ & + \int_{\Gamma} (\lambda_D - \lambda) \cdot ([[\mathbf{u}]] - [[u_D]]_D). \end{aligned} \quad (52) $$

Recalling that $\lambda_D \in C_D \subset C_f$, we deduce from (2b) that $\int_{\Gamma} (\lambda_D - \lambda) \cdot [[\mathbf{u}]] \le 0$. Furthermore, as $\pi_{\mathcal{F}_\Gamma}^0 \lambda \in C_D$ (since g is assumed to be piecewise constant on $\mathcal{F}_\Gamma$), we obtain from the fact that $[[u_D]]_D$ is piecewise constant on $\mathcal{F}_\Gamma$ and (10b) that

$$ \int_{\Gamma} (\lambda - \lambda_D) \cdot [[u_D]]_D = \int_{\Gamma} (\pi_{\mathcal{F}_\Gamma}^0 \lambda - \lambda_D) \cdot [[u_D]]_D \le 0, $$

and consequently

$$ \int_{\Gamma} (\lambda_D - \lambda) \cdot ([[\mathbf{u}]] - [[u_D]]_D) \le 0. $$

Plugging this relation into (52), using the norm estimate (49), the definitions (15) of $C_D$ and (17) of $W_D$, a Cauchy–Schwarz inequality on $S_{\mu, \lambda, D}$ and Young's inequality, we infer

$$ \|v_D - u_D\|_{1, D}^2 \lesssim W_D(\sigma(\mathbf{u}))^2 + C_D(\mathbf{u}, v_D)^2 + \int_{\Gamma} (\lambda_D - \lambda) \cdot ([[v_D]]_D - [[\mathbf{u}]]). \quad (53) $$

We now return to (51), which shows that, for all $w_D, v_D \in \mathbf{U}_{0,D}$,

$$ \begin{aligned} \int_{\Gamma} (\lambda_D - \lambda) \cdot [[w_D]]_D = & w_D(\sigma(\mathbf{u}), w_D) + \int_{\Omega} (\sigma(\mathbf{u}) - \sigma_D(\mathbf{u}_D)) : \mathcal{E}_D(w_D) \\ & - S_{\mu, \lambda, D}(v_D, w_D) - S_{\mu, \lambda, D}(u_D - v_D, w_D). \end{aligned} \quad (54) $$

Setting $\mu_D = \pi_{\mathcal{F}_\Gamma}^0 \lambda$ and noticing that $\int_{\Gamma} \lambda \cdot [[w_D]]_D = \int_{\Gamma} \mu_D \cdot [[w_D]]_D$ since $[[w_D]]_D$ is piecewise constant on $\mathcal{F}_\Gamma$, we deduce from (54) that

$$ \begin{aligned} \int_{\Gamma} (\lambda_D - \mu_D) \cdot [[w_D]]_D = & w_D(\sigma(\mathbf{u}), w_D) + \int_{\Omega} (\sigma(\mathbf{u}) - \sigma_D(\mathbf{u}_D)) : \mathcal{E}_D(w_D) \\ & - S_{\mu, \lambda, D}(v_D, w_D) - S_{\mu, \lambda, D}(u_D - v_D, w_D). \end{aligned} $$

The discrete inf-sup condition (20) then leads to

$$ \begin{aligned} \|\lambda_D - \mu_D\|_{-1/2, \Gamma} &\lesssim \sup_{\mathbf{w}_D \in \mathbf{U}_{0,D}} \frac{\int_{\Gamma} (\lambda_D - \mu_D) \cdot [[\mathbf{w}_D]]_D}{\|\mathbf{w}_D\|_{1, D}} \\ &\lesssim W_D(\sigma(\mathbf{u})) + \|\sigma(\mathbf{u}) - \sigma_D(\mathbf{v}_D)\|_{L^2(\Omega \setminus \Gamma)} + \|\sigma_D(\mathbf{u}_D - \mathbf{v}_D)\|_{L^2(\Omega \setminus \Gamma)} \\ &\quad + S_{\mu, \lambda, D}(\mathbf{v}_D, \mathbf{v}_D)^{1/2} + S_{\mu, \lambda, D}(\mathbf{u}_D - \mathbf{v}_D, \mathbf{u}_D - \mathbf{v}_D)^{1/2}, \end{aligned} $$
In the second inequality, we have used the fact that, by mesh regularity property, any $\mathbf{x} \in \sigma$ can be written $\mathbf{x} = \sum_{s \in \mathcal{V}_{\sigma}} \rho_s \mathbf{x}_s$ with $\sum_{s \in \mathcal{V}_{\sigma}} \rho_s = 1$ and $\sum_{s \in \mathcal{V}_{\sigma}} |\rho_s| \lesssim 1$, and that $\llbracket \Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}} \rrbracket(\mathbf{x}) = \sum_{s \in \mathcal{V}_{\sigma}} \rho_s \llbracket \Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}} \rrbracket(\mathbf{x}_s)$. We deduce from (48) that

$$h_{\sigma}^{-1} \| \llbracket \Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}} \rrbracket |_{\sigma} \|_{L^2(\sigma)}^2 \lesssim S_K(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) + S_L(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}). \quad (49)$$

Consider now a face $\sigma \in \mathcal{F}^{\text{bd}} \cap \mathcal{F}_K$, $K \in \mathcal{M}$. For $s \in \mathcal{V}_{\sigma}$, noticing that $\mathbf{v}_{\mathcal{K}_s} = 0$, we obtain

$$|\llbracket \Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}} \rrbracket|_{\sigma}(\mathbf{x}_s)|^2 = |(\Pi^K \mathbf{v}_{\mathcal{D}})(\mathbf{x}_s) - \mathbf{v}_{\mathcal{K}_s}|^2 \le h_K^{2-d} S_K(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}).$$

Then, following the same arguments that lead to (49), we get

$$h_{\sigma}^{-1} \| \llbracket \Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}} \rrbracket |_{\sigma} \|_{L^2(\sigma)}^2 \lesssim S_K(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}). \quad (50)$$

Summing (49) over the faces $\sigma \in \mathcal{F}^{\text{in}} \setminus \mathcal{F}_{\Gamma}$ and (50) over the faces $\sigma \in \mathcal{F}^{\text{bd}}$, we deduce

$$\sum_{\sigma \in \mathcal{F} \setminus \mathcal{F}_{\Gamma}} h_{\sigma}^{-1} \| \llbracket \Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}} \rrbracket |_{\sigma} \|_{L^2(\sigma)}^2 \lesssim \sum_{K \in \mathcal{M}} S_K(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}). \quad (51)$$

The proof of (45) is concluded by noticing that

$$\| \mathbf{v}_{\mathcal{D}} \|_{1, \mathcal{D}}^2 = \| \nabla_{\mathcal{M}} (\Pi^{\mathcal{D}} \mathbf{v}_{\mathcal{D}}) \|_{L^2(\Omega \setminus \Gamma)}^2 + \sum_{K \in \mathcal{M}} S_K(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}),$$

and by using (46) and (51).

□

## 5.3 Proof of the abstract error estimate

We prove here Theorem 4.5. To shorten the notations, let us recall the definition (14) of the discrete energy inner product $\langle \cdot, \cdot \rangle_{e, \mathcal{D}}$ and denote by $\| \cdot \|_{e, \mathcal{D}}$ its associated norm. From the discrete Korn inequality (45), we deduce the following bound for all $\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0, \mathcal{D}}$:

$$\| \mathbf{v}_{\mathcal{D}} \|_{1, \mathcal{D}} \lesssim \| \mathbf{v}_{\mathcal{D}} \|_{e, \mathcal{D}}. \quad (52)$$

Since $\mathbf{u}$ is a weak solution of the contact problem and by the regularity assumed on $\lambda$ in the theorem, we have $-\text{div}(\sigma(\mathbf{u})) = \mathbf{f} \in L^2(\Omega)$ and $\lambda = -\gamma_{\mathbf{n}}^{\dagger} \sigma(\mathbf{u}) \in L^2(\Gamma)$. Hence, $\sigma(\mathbf{u}) \in \Sigma$ (see (19)) and, using the definition (20) of $w_{\mathcal{D}}$, we have for all $\mathbf{w}_{\mathcal{D}} \in \mathbf{U}_{0, \mathcal{D}}$

$$\int_{\Omega} \sigma(\mathbf{u}) : \varepsilon^{\mathcal{D}}(\mathbf{w}_{\mathcal{D}}) - \int_{\Omega} \mathbf{f} \cdot \tilde{\Pi}^{\mathcal{D}} \mathbf{w}_{\mathcal{D}} = -w_{\mathcal{D}}(\sigma(\mathbf{u}), \mathbf{w}_{\mathcal{D}}) - \int_{\Gamma} \lambda \cdot \llbracket \mathbf{w}_{\mathcal{D}} \rrbracket_{\mathcal{D}}. \quad (53)$$

Subtracting (11a) (with $\mathbf{v}_{\mathcal{D}} = \mathbf{w}_{\mathcal{D}}$) from (53), we obtain

$$\int_{\Omega} (\sigma(\mathbf{u}) - \sigma^{\mathcal{D}}(\mathbf{u}_{\mathcal{D}})) : \varepsilon^{\mathcal{D}}(\mathbf{w}_{\mathcal{D}}) - S_{E, \nu, \mathcal{D}}(\mathbf{u}_{\mathcal{D}}, \mathbf{w}_{\mathcal{D}}) = -w_{\mathcal{D}}(\sigma(\mathbf{u}), \mathbf{w}_{\mathcal{D}}) + \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot \llbracket \mathbf{w}_{\mathcal{D}} \rrbracket_{\mathcal{D}}. \quad (54)$$

Take $\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0, \mathcal{D}}$ and set $\mathbf{w}_{\mathcal{D}} = \mathbf{v}_{\mathcal{D}} - \mathbf{u}_{\mathcal{D}}$ in (54) to get

$$\| \mathbf{v}_{\mathcal{D}} - \mathbf{u}_{\mathcal{D}} \|_{e, \mathcal{D}}^2 = - w_{\mathcal{D}}(\sigma(\mathbf{u}), \mathbf{v}_{\mathcal{D}} - \mathbf{u}_{\mathcal{D}}) - \int_{\Omega} (\sigma(\mathbf{u}) - \sigma^{\mathcal{D}}(\mathbf{v}_{\mathcal{D}})) : \varepsilon^{\mathcal{D}}(\mathbf{v}_{\mathcal{D}} - \mathbf{u}_{\mathcal{D}}) \\ + S_{E, \nu, \mathcal{D}}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}} - \mathbf{u}_{\mathcal{D}}) + \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket_{\mathcal{D}})$$
from which we infer

$$ \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma} \lesssim \|\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}\|_{1, \mathcal{D}} + \mathcal{W}_{\mathcal{D}}(\sigma(\mathbf{u})) + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}}). \quad (55) $$

Combining (53) and (55) yields

$$ \begin{aligned} \|\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}\|_{1, \mathcal{D}}^2 + \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma}^2 &\lesssim \mathcal{W}_{\mathcal{D}}(\sigma(\mathbf{u}))^2 + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}})^2 \\ &\quad + \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket). \end{aligned} \quad (56) $$

Let us introduce the discrete $H_0^{1/2}(\Gamma)$-norm, dual of the discrete $H_0^{1/2}(\Gamma)$-norm (18) and defined for all $\mu \in L^2(\Gamma)^d$ by

$$ \|\mu\|_{-1/2, \mathcal{D}} = \left( \sum_{\sigma \in \mathcal{F}_{\Gamma}} h_{\sigma} \|\mu\|_{L^2(\sigma)}^2 \right)^{1/2}. \quad (57) $$

Using a weighted Cauchy–Schwarz inequality, we note that

$$ \int_{\Gamma} \mu \cdot \xi \le \|\mu\|_{-1/2, \mathcal{D}} \|\xi\|_{1/2, \mathcal{D}} \quad \forall \mu, \xi \in L^2(\Gamma)^d. \quad (58) $$

Then, the contact term in the right-hand side of (56) can be estimated by introducing $\mu_{\mathcal{D}} = \pi_{\mathcal{F}_{\Gamma}}^0 \lambda$ and writing

$$ \begin{aligned} \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket) &= \int_{\Gamma} (\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \pi_{\mathcal{F}_{\Gamma}}^0 \llbracket \mathbf{u} \rrbracket) + \int_{\Gamma} (\mu_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket) \\ &\le \frac{c}{2} \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \mathcal{D}}^2 + \frac{1}{2c} \|\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \pi_{\mathcal{F}_{\Gamma}}^0 \llbracket \mathbf{u} \rrbracket\|_{1/2, \mathcal{D}}^2 \\ &\quad + \|\mu_{\mathcal{D}} - \lambda\|_{L^2(\Gamma)} \|\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket\|_{L^2(\Gamma)}, \end{aligned} $$

where the introduction of the projector $\pi_{\mathcal{F}_{\Gamma}}^0$ in front of $\llbracket \mathbf{u} \rrbracket$ in the first line is justified by the fact that $\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}$ is piecewise constant on $\mathcal{F}_{\Gamma}$, and the conclusion follows from Cauchy–Schwarz inequalities (including (58)) and a Young inequality (with an arbitrary $c > 0$). We then apply Lemma B.2 in the appendix to write $\|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \mathcal{D}}^2 \lesssim \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma}^2$ and plug the resulting inequality in (56). Adjusting $c$ to absorb $\|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma}^2$ in the left-hand side, we infer

$$ \begin{aligned} \|\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}\|_{1, \mathcal{D}}^2 + \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma}^2 &\lesssim \|\mu_{\mathcal{D}} - \lambda\|_{L^2(\Gamma)} \|\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket\|_{L^2(\Gamma)} \\ &\quad + \|\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \pi_{\mathcal{F}_{\Gamma}}^0 \llbracket \mathbf{u} \rrbracket\|_{1/2, \mathcal{D}}^2 + \mathcal{W}_{\mathcal{D}}(\sigma(\mathbf{u}))^2 + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}})^2. \end{aligned} \quad (59) $$

Applying triangular inequalities we can write

$$ \begin{aligned} &\|\nabla^{\mathcal{D}} \mathbf{u}_{\mathcal{D}} - \nabla \mathbf{u}\|_{L^2(\Omega \setminus \Gamma)} + \|\lambda_{\mathcal{D}} - \lambda\|_{-1/2, \Gamma} \\ &\lesssim \|\nabla^{\mathcal{D}} \mathbf{u}_{\mathcal{D}} - \nabla^{\mathcal{D}} \mathbf{v}_{\mathcal{D}}\|_{L^2(\Omega \setminus \Gamma)} + \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma} + \|\nabla^{\mathcal{D}} \mathbf{v}_{\mathcal{D}} - \nabla \mathbf{u}\|_{L^2(\Omega \setminus \Gamma)} + \|\mu_{\mathcal{D}} - \lambda\|_{-1/2, \Gamma} \\ &\lesssim \|\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}\|_{1, \mathcal{D}} + \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma} + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}}) + \|\mu_{\mathcal{D}} - \lambda\|_{-1/2, \Gamma}. \end{aligned} \quad (60) $$

Taking the square root of (59), plugging the resulting estimate into (60), recalling that $\mu_{\mathcal{D}} = \pi_{\mathcal{F}_{\Gamma}}^0 \lambda$, then taking the infimum over $\mathbf{v}_{\mathcal{D}}$ concludes the proof of (19).
$$ + \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{u} \rrbracket - \llbracket \mathbf{u}_{\mathcal{D}} \rrbracket_{\mathcal{D}}). \qquad (55) $$

Recalling that $\lambda_{\mathcal{D}} \in \mathbf{C}_{\mathcal{D}} \subset C_f$, we deduce from (3b) that $\int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot \llbracket \mathbf{u} \rrbracket \le 0$. Furthermore, as $\pi_{\mathcal{F}_{\Gamma}}^0 \lambda \in \mathbf{C}_{\mathcal{D}}$ (since g is assumed to be piecewise constant on $\mathcal{F}_{\Gamma}$), we obtain from the fact that $\llbracket \mathbf{u}_{\mathcal{D}} \rrbracket_{\mathcal{D}}$ is piecewise constant on $\mathcal{F}_{\Gamma}$ and (11b) that

$$ \int_{\Gamma} (\lambda - \lambda_{\mathcal{D}}) \cdot \llbracket \mathbf{u}_{\mathcal{D}} \rrbracket_{\mathcal{D}} = \int_{\Gamma} (\pi_{\mathcal{F}_{\Gamma}}^0 \lambda - \lambda_{\mathcal{D}}) \cdot \llbracket \mathbf{u}_{\mathcal{D}} \rrbracket_{\mathcal{D}} \le 0, $$

and consequently

$$ \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{u} \rrbracket - \llbracket \mathbf{u}_{\mathcal{D}} \rrbracket_{\mathcal{D}}) \le 0. $$

Plugging this relation into (55), using the norm estimate (52), the definitions (18) of $C_{\mathcal{D}}$ and (20) of $W_{\mathcal{D}}$, a Cauchy–Schwarz inequality on $S_{E,\nu,\mathcal{D}}$ and the Young inequality, we infer

$$ \|\mathbf{v}_{\mathcal{D}} - \mathbf{u}_{\mathcal{D}}\|_{1,\mathcal{D}}^2 \lesssim W_{\mathcal{D}}(\sigma(\mathbf{u}))^2 + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}})^2 + \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket). \quad (56) $$

We now return to (54), which shows that, for all $\mathbf{w}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}$,

$$ \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot \llbracket \mathbf{w}_{\mathcal{D}} \rrbracket_{\mathcal{D}} = w_{\mathcal{D}}(\sigma(\mathbf{u}), \mathbf{w}_{\mathcal{D}}) + \int_{\Omega} (\sigma(\mathbf{u}) - \sigma^{\mathcal{D}}(\mathbf{u}_{\mathcal{D}})) : \varepsilon^{\mathcal{D}}(\mathbf{w}_{\mathcal{D}}) \\ - S_{E,\nu,\mathcal{D}}(\mathbf{v}_{\mathcal{D}}, \mathbf{w}_{\mathcal{D}}) - S_{E,\nu,\mathcal{D}}(\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}, \mathbf{w}_{\mathcal{D}}). \quad (57) $$

Setting $\mu_{\mathcal{D}} = \pi_{\mathcal{F}_{\Gamma}}^0 \lambda$ and noticing that $\int_{\Gamma} \lambda \cdot \llbracket \mathbf{w}_{\mathcal{D}} \rrbracket_{\mathcal{D}} = \int_{\Gamma} \mu_{\mathcal{D}} \cdot \llbracket \mathbf{w}_{\mathcal{D}} \rrbracket_{\mathcal{D}}$ since $\llbracket \mathbf{w}_{\mathcal{D}} \rrbracket_{\mathcal{D}}$ is piecewise constant on $\mathcal{F}_{\Gamma}$, we deduce from (57) that

$$ \int_{\Gamma} (\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}) \cdot \llbracket \mathbf{w}_{\mathcal{D}} \rrbracket_{\mathcal{D}} = w_{\mathcal{D}}(\sigma(\mathbf{u}), \mathbf{w}_{\mathcal{D}}) + \int_{\Omega} (\sigma(\mathbf{u}) - \sigma^{\mathcal{D}}(\mathbf{u}_{\mathcal{D}})) : \varepsilon^{\mathcal{D}}(\mathbf{w}_{\mathcal{D}}) \\ - S_{E,\nu,\mathcal{D}}(\mathbf{v}_{\mathcal{D}}, \mathbf{w}_{\mathcal{D}}) - S_{E,\nu,\mathcal{D}}(\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}, \mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}). $$

The discrete inf-sup condition (23) then leads to

$$ \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma} \lesssim \sup_{\mathbf{w}_{\mathcal{D}} \in \mathbf{U}_{0, \mathcal{D}}} \frac{\int_{\Gamma} (\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}) \cdot \llbracket \mathbf{w}_{\mathcal{D}} \rrbracket_{\mathcal{D}}}{\|\mathbf{w}_{\mathcal{D}}\|_{1, \mathcal{D}}} \\ \lesssim W_{\mathcal{D}}(\sigma(\mathbf{u})) + \|\sigma(\mathbf{u}) - \sigma^{\mathcal{D}}(\mathbf{v}_{\mathcal{D}})\|_{L^2(\Omega \setminus \bar{\Gamma})} + \|\sigma^{\mathcal{D}}(\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}})\|_{L^2(\Omega \setminus \bar{\Gamma})} \\ + S_{E, \nu, \mathcal{D}}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}})^{1/2} + S_{E, \nu, \mathcal{D}}(\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}, \mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}})^{1/2}, $$

from which we infer

$$ \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma} \lesssim \|\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}\|_{1, \mathcal{D}} + W_{\mathcal{D}}(\sigma(\mathbf{u})) + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}}). \quad (58) $$

Combining (56) and (58) yields

$$ \|\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}\|_{1, \mathcal{D}}^2 + \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma}^2 \lesssim W_{\mathcal{D}}(\sigma(\mathbf{u}))^2 + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}})^2 \\ + \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket). \quad (59) $$

Let us introduce the discrete $H^{-1/2}$-norm, dual of the discrete $H^{1/2}$-norm (21) and defined for all $\mu \in L^2(\Gamma)$ by

$$ \|\mu\|_{-1/2, \mathcal{D}} = \left( \sum_{\sigma \in \mathcal{F}_{\Gamma}} h_{\sigma} \|\mu\|_{L^2(\sigma)}^2 \right)^{1/2}. \quad (60) $$
## 5.4 Proof of the error estimate

Theorem 4.6 directly follows from the abstract error estimate (19), Lemma B.1 (with $s = 1$) in the appendix, and Lemmas 5.8, 5.9 and 5.10 below.

**Lemma 5.8 (Consistency of the gradient reconstruction).** If $\mathbf{u} \in \mathbf{U}_0 \cap H^2(\mathcal{M})^d$ then, recalling the definition (15) of $C_{\mathcal{D}}$, it holds

$$C_{\mathcal{D}}(\mathbf{u}, \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u}) \lesssim h_{\mathcal{D}} |\mathbf{u}|_{H^2(\mathcal{M})} \quad (61)$$

*Proof.* We first note that the regularity assumption ensures that $\mathbf{u} \in C_0^0(\bar{\Omega} \setminus \Gamma)$: the $H^2$ regularity ensures the continuity of $\mathbf{u}$ on each $\bar{K}$, and the values on each side of each $\sigma \notin \mathcal{F}_\Gamma$ match by the $H^1$ regularity. Hence, $\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u}$ is well-defined.

Let $K \in \mathcal{M}$ and $\mathbf{q}$ be the $L^2(K)$-orthogonal projection of $\mathbf{u}$ on $\mathbb{P}^1(K)^d$. By the approximation properties of the polynomial projector [16, Theorem 1.45], we have

$$|\mathbf{u} - \mathbf{q}|_{H^s(K)} \lesssim h_K^s |\mathbf{u}|_{H^2(K)}, \quad \forall s \in \{0, 1, 2\}. \quad (62)$$

Applying the bound [16, Eq. (5.110)] to $\mathbf{u} - \mathbf{q}$ yields

$$\begin{aligned} \max_{\bar{K}} |\mathbf{u} - \mathbf{q}| & \lesssim |K|^{-1/2} \left( \|\mathbf{u} - \mathbf{q}\|_{L^2(K)} + h_K |\mathbf{u} - \mathbf{q}|_{H^1(K)} + h_K^2 |\mathbf{u} - \mathbf{q}|_{H^2(K)} \right) \\ & \lesssim |K|^{-1/2} h_K^2 |\mathbf{u}|_{H^2(K)}. \end{aligned} \quad (63)$$

where the conclusion follows from (62). Plugging this estimate into the bound (24) for $w_K$ gathering the DOFs on $K$ of $\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} (\mathbf{u} - \mathbf{q})$ and recalling the definition (13) of the local norm, we infer

$$(\|\nabla^K \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} (\mathbf{u} - \mathbf{q})\|_{L^2(K)}^2 + S_K (\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u} - \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{q}, \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u} - \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{q}))^{1/2} \lesssim h_K |\mathbf{u}|_{H^2(K)}.$$

The linear exactness properties (33) and (34) allow us to manipulate the left-hand side to obtain

$$(\|\nabla^K \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u} - \nabla \mathbf{q}\|_{L^2(K)}^2 + S_K (\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u}, \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u}))^{1/2} \lesssim h_K |\mathbf{u}|_{H^2(K)}.$$

Recalling the definition (15) of $C_{\mathcal{D}}(\mathbf{u}, \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u})$, the estimate (61) follows by introducing $\pm \nabla \mathbf{u}$ in the left-hand side above, using a triangle inequality, invoking again the approximation property (62) (with $s = 1$), squaring, summing over $K$, using the bound $h_K \le h_{\mathcal{D}}$ and taking the square root. $\square$

**Lemma 5.9 (Consistency of the jump reconstructions).** If $\mathbf{u} \in H^2(\mathcal{M})^d$ then

$$\|\|\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u}\|_{\mathcal{D}} - \|\mathbf{u}\|_{L^2(\Gamma)}\| \lesssim h_{\mathcal{D}} \left( \|\|\mathbf{u}\|_{H^1(\mathcal{F}_\Gamma)} + h_{\mathcal{D}}^{1/2} |\mathbf{u}|_{H^2(\mathcal{M})} \right) \quad (64)$$

and

$$\|\|\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u}\|_{\mathcal{D}} - \pi_{\mathcal{F}_\Gamma}^0 \|\mathbf{u}\|_{1/2, \mathcal{D}} \lesssim h_{\mathcal{D}} |\mathbf{u}|_{H^2(\mathcal{M})} \quad (65)$$

*Proof.* These approximation properties are, similarly to (61), a consequence of the linear exactness of the jump reconstruction (upon projection on piecewise constant functions) and of a bound on this operator.
Using a weighted Cauchy–Schwarz inequality, we note that

$$ \int_{\Gamma} \mu \cdot \xi \leq \|\mu\|_{-1/2, \mathcal{D}} \|\xi\|_{1/2, \mathcal{D}} \quad \forall \mu, \xi \in L^2(\Gamma). \qquad (61) $$

Then, the contact term in the right-hand side of (59) can be estimated by introducing $\mu_{\mathcal{D}} = \pi_{\mathcal{F}\Gamma}^0 \lambda$ and writing

$$ \begin{aligned} \int_{\Gamma} (\lambda_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket) &= \int_{\Gamma} (\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \pi_{\mathcal{F}\Gamma}^0 \llbracket \mathbf{u} \rrbracket) + \int_{\Gamma} (\mu_{\mathcal{D}} - \lambda) \cdot (\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket) \\ &\leq \frac{c}{2} \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \mathcal{D}}^2 + \frac{1}{2c} \|\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \pi_{\mathcal{F}\Gamma}^0 \llbracket \mathbf{u} \rrbracket\|_{1/2, \mathcal{D}}^2 \\ &\quad + \|\mu_{\mathcal{D}} - \lambda\|_{L^2(\Gamma)} \|\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket\|_{L^2(\Gamma)}, \end{aligned} $$

where the introduction of the projector $\pi_{\mathcal{F}\Gamma}^0$ in front of $\llbracket \mathbf{u} \rrbracket$ in the first line is justified by the fact that $\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}$ is piecewise constant on $\mathcal{F}\Gamma$, and the conclusion follows from Cauchy–Schwarz inequalities (including (61)) and a Young inequality (with an arbitrary $c > 0$). We then apply Lemma B.2 in the appendix to write $\|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \mathcal{D}}^2 \lesssim \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma}^2$ and plug the resulting inequality in (59). Adjusting $c$ to absorb $\|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma}^2$ in the left-hand side, we infer

$$ \begin{aligned} \|\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}\|_{1, \mathcal{D}}^2 + \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma}^2 &\lesssim \|\mu_{\mathcal{D}} - \lambda\|_{L^2(\Gamma)} \|\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \llbracket \mathbf{u} \rrbracket\|_{L^2(\Gamma)} \\ &\quad + \|\llbracket \mathbf{v}_{\mathcal{D}} \rrbracket_{\mathcal{D}} - \pi_{\mathcal{F}\Gamma}^0 \llbracket \mathbf{u} \rrbracket\|_{1/2, \mathcal{D}}^2 + \mathcal{W}_{\mathcal{D}}(\sigma(\mathbf{u}))^2 + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}})^2. \end{aligned} \qquad (62) $$

Applying triangular inequalities we can write

$$ \begin{aligned} & \|\nabla^{\mathcal{D}} \mathbf{u}_{\mathcal{D}} - \nabla \mathbf{u}\|_{L^2(\Omega \setminus \Gamma)} + \|\lambda_{\mathcal{D}} - \lambda\|_{-1/2, \Gamma} \\ &\lesssim \|\nabla^{\mathcal{D}} \mathbf{u}_{\mathcal{D}} - \nabla^{\mathcal{D}} \mathbf{v}_{\mathcal{D}}\|_{L^2(\Omega \setminus \Gamma)} + \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma} + \|\nabla^{\mathcal{D}} \mathbf{v}_{\mathcal{D}} - \nabla \mathbf{u}\|_{L^2(\Omega \setminus \Gamma)} + \|\mu_{\mathcal{D}} - \lambda\|_{-1/2, \Gamma} \\ &\lesssim \|\mathbf{u}_{\mathcal{D}} - \mathbf{v}_{\mathcal{D}}\|_{1, \mathcal{D}} + \|\lambda_{\mathcal{D}} - \mu_{\mathcal{D}}\|_{-1/2, \Gamma} + C_{\mathcal{D}}(\mathbf{u}, \mathbf{v}_{\mathcal{D}}) + \|\mu_{\mathcal{D}} - \lambda\|_{-1/2, \Gamma}. \end{aligned} \qquad (63) $$

Taking the square root of (62), plugging the resulting estimate into (63), recalling that $\mu_{\mathcal{D}} = \pi_{\mathcal{F}\Gamma}^0 \lambda$, then taking the infimum over $\mathbf{v}_{\mathcal{D}}$ concludes the proof of (22).

## 5.4 Proof of the error estimate

Theorem 4.6 directly follows from the abstract error estimate (22), Lemma B.1 (with $s = 1$) in the appendix, and Lemmas 5.8, 5.9 and 5.10 below.

**Lemma 5.8 (Consistency of the gradient reconstruction).** If $\mathbf{u} \in \mathbf{U}_0 \cap H^2(\mathcal{M})$ then, recalling the definition (18) of $C_{\mathcal{D}}$, it holds

$$ C_{\mathcal{D}}(\mathbf{u}, \mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u}) \lesssim h_{\mathcal{D}} |\mathbf{u}|_{H^2(\mathcal{M})}. \qquad (64) $$

*Proof.* We first note that the regularity assumption ensures that $\mathbf{u} \in C_0^0(\bar{\Omega} \setminus \Gamma)$: the $H^2$ regularity ensures the continuity of $\mathbf{u}$ on each $\bar{K}$, and the values on each side of each $\sigma \notin \mathcal{F}\Gamma$ match by the $H^1$ regularity. Hence, $\mathcal{I}_{\mathbf{U}_{0,\mathcal{D}}} \mathbf{u}$ is well-defined.

Let $K \in \mathcal{M}$ and $q$ be the $L^2(K)$-orthogonal projection of $\mathbf{u}$ on $\mathbb{P}^1(K)^d$. By the approximation properties of the polynomial projector [16, Theorem 1.45], we have

$$ |\mathbf{u} - \mathbf{q}|_{H^s(K)} \lesssim h_K^{2-s} |\mathbf{u}|_{H^2(K)}, \quad \forall s \in \{0, 1, 2\}. \qquad (65) $$
The definition (6) of the jump reconstruction directly gives

$$ \|[\mathbf{v}_D]\sigma\| \le \max_{s \in \mathcal{V}_\sigma} |\mathbf{v}_{\mathcal{K}s}| + \max_{s \in \mathcal{V}_\sigma} |\mathbf{v}_{\mathcal{L}s}| + |\mathbf{v}_{\mathcal{K}\sigma}| \quad \forall \mathbf{v}_D \in \mathbf{U}_D, \forall \sigma = K|L| \in \mathcal{F}_\Gamma. \quad (66) $$

Consider now, for such a $\sigma$, two linear functions $\mathbf{q}_K \in \mathbb{P}^1(K)^d$, resp. $\mathbf{q}_L \in \mathbb{P}^1(L)^d$, that satisfy (63) for $K$, resp. for $L$. Let $\mathbf{q}$ be the piecewise polynomial on $K \cup L$ defined by $\mathbf{q}_K$ and $\mathbf{q}_L$. Applying (63) on each side of $\sigma$ we have

$$ \|[\mathbf{u} - \mathbf{q}]\|_{L^\infty(\sigma)} \le |K|^{-1/2} h_K^2 |\mathbf{u}|_{H^2(K)} + |L|^{-1/2} h_L^2 |\mathbf{u}|_{H^2(L)}. \quad (67) $$

Moreover, even though $\mathbf{q}$ is not globally defined or in $\mathbf{U}_0$, its interpolate on the degrees of freedom attached to $\sigma$ is well defined. The definition (9) of $\mathcal{I}_{\mathbf{U}_{0,D}}$ easily shows that the degrees of freedom of $\mathcal{I}_{\mathbf{U}_{0,D}}(\mathbf{u} - \mathbf{q})$ on $\sigma$ are bounded by the maximum of $|\mathbf{u} - \mathbf{q}|$ on $K \cup L$. Hence, by (66) and (63),

$$ \|[\mathcal{I}_{\mathbf{U}_{0,D}}(\mathbf{u} - \mathbf{q})]\|_{\sigma} \le |K|^{-1/2} h_K^2 |\mathbf{u}|_{H^2(K)} + |L|^{-1/2} h_L^2 |\mathbf{u}|_{H^2(L)}. \quad (68) $$

The linear exactness (31) applied on each side of $\sigma$ together with (32) shows that $[\mathcal{I}_{\mathbf{U}_{0,D}}\mathbf{q}]_{\sigma} = \pi_{\sigma}^{0}[\mathbf{q}]$ (with $\pi_{\sigma}^{0}$ the $L^2$-projector on $\mathbb{P}^0(\sigma)^d$). We can therefore write

$$ \begin{aligned} \|[\mathcal{I}_{\mathbf{U}_{0,D}} \mathbf{u}]_{\sigma} - [\mathbf{u}]\|_{L^2(\sigma)} &\le \|[\mathcal{I}_{\mathbf{U}_{0,D}}(\mathbf{u} - \mathbf{q})]\|_{L^2(\sigma)} + \|\pi_{\sigma}^{0}[\mathbf{q} - \mathbf{u}]\|_{L^2(\sigma)} + \|\pi_{\sigma}^{0}[\mathbf{u}] - [\mathbf{u}]\|_{L^2(\sigma)} \\ &\le |\sigma|^{1/2} \left( |K|^{-1/2} h_K^2 |\mathbf{u}|_{H^2(K)} + |L|^{-1/2} h_L^2 |\mathbf{u}|_{H^2(L)} \right) \\ &\quad + \|\pi_{\sigma}^{0}[\mathbf{u}] - [\mathbf{u}]\|_{L^2(\sigma)} \\ &\le h_D^{3/2} \left( |\mathbf{u}|_{H^2(K)} + |\mathbf{u}|_{H^2(L)} \right) + h_{\sigma} \|[\mathbf{u}]\|_{H^1(\sigma)}, \end{aligned} \quad (69) $$

where we have introduced $\pm\pi_{\sigma}^{0}[\mathbf{q} - \mathbf{u}] = \pm([\mathcal{I}_{\mathbf{U}_{0,D}}\mathbf{q}]_{\sigma} - \pi_{\sigma}^{0}[\mathbf{u}])$ and used a triangle inequality in the first line, invoked (68) and used the $L^2(\sigma)$-boundedness of $\pi_{\sigma}^{0}$ together with (67) in the second inequality, and used the mesh regularity property in the conclusion to write $|\sigma|^{1/2}h_K^{1/2} \le |K|^{1/2}$ and $|\sigma|^{1/2}h_L^{1/2} \le |L|^{1/2}$, together with the approximation properties of $\pi_{\sigma}^{0}$ [16, Theorem 1.45]. Squaring, summing over $\sigma$ and taking the square root yields (64).

To prove the second estimate, we first notice that, by definition (6) of $[\cdot]_{\sigma}$ and (9) of the interpolator,

$$ [\mathcal{I}_{\mathbf{U}_{0,D}} \mathbf{u}]_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} (\gamma^{K\sigma} \mathbf{u} - \Pi^{L\sigma} (\mathcal{I}_{\mathbf{U}_{0,D}} \mathbf{u})) = \pi_{\sigma}^{0} (\gamma^{K\sigma} \mathbf{u} - \Pi^{L\sigma} (\mathcal{I}_{\mathbf{U}_{0,D}} \mathbf{u})). $$

Hence,

$$ \begin{aligned} \|[\mathcal{I}_{\mathbf{U}_{0,D}} \mathbf{u}]_{\sigma} - \pi_{\sigma}^{0} [\mathbf{u}]\|_{L^2(\sigma)} &= \|\pi_{\sigma}^{0} (\gamma^{K\sigma} \mathbf{u} - \Pi^{L\sigma} (\mathcal{I}_{\mathbf{U}_{0,D}} \mathbf{u}) - [\mathbf{u}]\|_{L^2(\sigma)} \\ &\le \|\gamma^{L\sigma} \mathbf{u} - \Pi^{L\sigma} (\mathcal{I}_{\mathbf{U}_{0,D}} \mathbf{u})\|_{L^2(\sigma)}. \end{aligned} $$

Introduce as above a linear approximation $\mathbf{q}_L$ of $\mathbf{u}$ in $L$, and use the linear exactness (31) (with $L$ instead of $K$) of $\Pi^{L\sigma}$ to deduce

$$ \|[\mathcal{I}_{\mathbf{U}_{0,D}} \mathbf{u}]_{\sigma} - \pi_{\sigma}^{0} [\mathbf{u}]\|_{L^2(\sigma)} \le \|\mathbf{u}|_L - \mathbf{q}_L\|_{L^2(\sigma)} + \|\Pi^{L\sigma} \mathcal{I}_{\mathbf{U}_{0,D}} (\mathbf{u}|_L - \mathbf{q}_L)\|_{L^2(\sigma)}. $$

In a similar way as (26) we can show that $\Pi^{L\sigma}\mathbf{w}_D$ is bounded above, on $\sigma$, by the maximum of the absolute values of the degrees of freedom associated with $\sigma$. By definition of $\mathcal{I}_{\mathbf{U}_{0,D}}$ and (63) (with $L$ instead of $K$) and using similar arguments as for the first terms in (69), we infer

$$ \|[\mathcal{I}_{\mathbf{U}_{0,D}} \mathbf{u}]_{\sigma} - \pi_{\sigma}^{0} [\mathbf{u}]\|_{L^2(\sigma)} \le h_L^{3/2} |\mathbf{u}|_{H^2(L)}. $$

Squaring, multiplying by $h_{\sigma}^{-1} \le h_L^{-1}$, summing over $\sigma$ and taking the square root yields (65). $\square$
Applying the bound [16, Eq. (5.110)] to $\mathbf{u}-\mathbf{q}$ yields

$$ \begin{aligned} \max_{\bar{K}} |\mathbf{u}-\mathbf{q}| &\lesssim |K|^{-1/2} (\|\mathbf{u}-\mathbf{q}\|_{L^2(K)} + h_K |\mathbf{u}-\mathbf{q}|_{H^1(K)} + h_K^2 |\mathbf{u}-\mathbf{q}|_{H^2(K)}) \\ &\lesssim |K|^{-1/2} h_K^2 |\mathbf{u}|_{H^2(K)}, \end{aligned} \quad (66) $$

where the conclusion follows from (65). Plugging this estimate into the bound (27) for $w_K$ gathering the DOFs on $K$ of $\mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u}-\mathbf{q})$ and recalling the definition (16) of the local norm, we infer

$$ (\|\nabla^K \mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u}-\mathbf{q})\|_{L^2(K)}^2 + S_K(\mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u}-\mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{q}), \mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u}-\mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{q})))^{1/2} \lesssim h_K |\mathbf{u}|_{H^2(K)}. $$

The linear exactness properties (36) and (37) allow us to manipulate the left-hand side to obtain

$$ (\|\nabla^K \mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u}-\nabla\mathbf{q})\|_{L^2(K)}^2 + S_K(\mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u}, \mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u})))^{1/2} \lesssim h_K |\mathbf{u}|_{H^2(K)}. $$

Recalling the definition (18) of $C_{\mathcal{D}}(\mathbf{u}, \mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u}))$, the estimate (64) follows by introducing $\pm\nabla\mathbf{u}$ in the left-hand side above, using a triangle inequality, invoking again the approximation property (65) (with $s=1$), squaring, summing over $K$, using the bound $h_K \le h_{\mathcal{D}}$ and taking the square root. $\square$

**Lemma 5.9 (Consistency of the jump reconstructions).** If $\mathbf{u} \in H^2(\mathcal{M})$ then

$$ \|\|\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{u}\|_{\mathcal{D}} - \|\|\mathbf{u}\|\|_{L^2(\Gamma)} \lesssim h_{\mathcal{D}} \left( \|\|\mathbf{u}\|\|_{H^1(\mathcal{T})} + h_{\mathcal{D}}^{1/2} |\mathbf{u}|_{H^2(\mathcal{M})} \right) \quad (67) $$

and

$$ \|\|\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{u}\|_{\mathcal{D}} - \pi_{\mathcal{T}}^0 \|\|\mathbf{u}\|\|_{1/2,\mathcal{D}} \lesssim h_{\mathcal{D}} |\mathbf{u}|_{H^2(\mathcal{M})}. \quad (68) $$

*Proof.* These approximation properties are, similarly to (64), a consequence of the linear exactness of the jump reconstruction (upon projection on piecewise constant functions) and of a bound on this operator.

The definition (7) of the jump reconstruction directly gives

$$ \|\|\mathbf{v}_{\mathcal{D}}\|\|_{\sigma} \lesssim \max_{s \in \mathcal{V}_{\sigma}} |\mathbf{v}_{\mathcal{K}_s}| + \max_{s \in \mathcal{V}_{\sigma}} |\mathbf{v}_{\mathcal{L}_s}| + |\mathbf{v}_{K\sigma}| \quad \forall \mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{\mathcal{D}}, \forall \sigma = K|L \in \mathcal{T}. \quad (69) $$

Consider now, for such a $\sigma$, two linear functions $\mathbf{q}_K \in \mathbb{P}^1(K)^d$, resp. $\mathbf{q}_L \in \mathbb{P}^1(L)^d$, that satisfy (66) for $K$, resp. for $L$. Let $\mathbf{q}$ be the piecewise polynomial on $K \cup L$ defined by $\mathbf{q}_K$ and $\mathbf{q}_L$. Applying (66) on each side of $\sigma$ we have

$$ \|\|\mathbf{u}-\mathbf{q}\|\|_{L^{\infty}(\sigma)} \lesssim |K|^{-1/2} h_K^2 |\mathbf{u}|_{H^2(K)} + |L|^{-1/2} h_L^2 |\mathbf{u}|_{H^2(L)}. \quad (70) $$

Moreover, even though $\mathbf{q}$ is not globally defined or in $U_0$, its interpolate on the degrees of freedom attached to $\sigma$ is well defined. The definition (10) of $\mathcal{I}_{U_{0,\mathcal{D}}}$ easily shows that the degrees of freedom of $\mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u}-\mathbf{q})$ on $\sigma$ are bounded by the maximum of $|\mathbf{u}-\mathbf{q}|$ on $K \cup L$. Hence, by (69) and (66),

$$ \|\|\mathcal{I}_{U_{0,\mathcal{D}}}(\mathbf{u}-\mathbf{q})\|\|_{\sigma} \lesssim |K|^{-1/2} h_K^2 |\mathbf{u}|_{H^2(K)} + |L|^{-1/2} h_L^2 |\mathbf{u}|_{H^2(L)}. \quad (71) $$

The linear exactness (34) applied on each side of $\sigma$ together with (35) shows that $\|\|\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{q}\|\|_{\sigma} = \pi_{\sigma}^0 \|\|\mathbf{q}\|\|$ (with $\pi_{\sigma}^0$ the $L^2$-projector on $\mathbb{P}^0(\sigma)^d$). We can therefore write

$$ \|\|\mathcal{I}_{U_{0,\mathcal{D}}} \mathbf{u}\|\|_{\sigma} - \|\|\mathbf{u}\|\|_{L^2(\sigma)} \lesssim \|\|\mathcal{I}_{U_{0,\mathcal{D}}} (\mathbf{u}-\mathbf{q})\|\|_{\sigma} \|_{L^2(\sigma)} + \|\pi_{\sigma}^0 \|\|\mathbf{q}-\mathbf{u}\|\|_{L^2(\sigma)} + \|\pi_{\sigma}^0 \|\|\mathbf{u}\| - \|\|\mathbf{u}\|\|_{L^2(\sigma)} $$
**Lemma 5.10 (Adjoint consistency).** *Recalling the definition (17) of the adjoint consistency error* $W_D$, we have for all $\nu \in W$ such that $\nu \in H^1(M)^{d \times d}$ the estimate

$$ W_D(\nu) \lesssim h_D |\nu|_{H^1(M)}. \qquad (70) $$

*Proof.* For each $\nu \in W \cap H^1(M)^{d \times d}$, the definitions (17) of $w_D$ and (6) of $[\cdot]_\sigma$ yield

$$ \begin{aligned} w_D(\nu, \mathbf{v}_D) &= \sum_{K \in \mathcal{M}} \left( -\left( \int_K \nu \right) : \mathfrak{e}_K(\mathbf{v}_D) - \sum_{\sigma \in \mathcal{F}_K} \bar{\mathbf{v}}_K \cdot \int_\sigma (\nu|_K \mathbf{n}_{K\sigma}) \right) \\ &\quad + \sum_{\sigma=K|L \in \mathcal{F}_\Gamma} \left( \int_\sigma \nu|_K \mathbf{n}_{K\sigma} \right) \cdot \left( \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma} + \mathbf{v}_{K\sigma} \right) \\ &= \sum_{K \in \mathcal{M}} \left( -|K| \nu_K : \mathfrak{e}_K(\mathbf{v}_D) - \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_K \cdot \tau_{K\sigma} \right) \\ &\quad + \sum_{\sigma=K|L \in \mathcal{F}_\Gamma} |\sigma| \tau_{K\sigma} \cdot \left( \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma} + \mathbf{v}_{K\sigma} \right), \end{aligned} $$

with

$$ \nu_K = \frac{1}{|K|} \int_K \nu \quad \text{and} \quad \tau_{K\sigma} = \frac{1}{|\sigma|} \int_\sigma (\nu|_K \mathbf{n}_{K\sigma}). $$

Recalling the definition (7) of $\nabla^K$ and noticing that $\nu_K : \mathfrak{e}_K(\mathbf{v}_D) = \nu_K : \nabla^K \mathbf{v}_D$ since $\nu_K$ is symmetric, we infer

$$ \begin{aligned} w_D(\nu, \mathbf{v}_D) &= \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} \left( -|\sigma| \bar{\mathbf{v}}_{K\sigma} \cdot (\nu_K \mathbf{n}_{K\sigma}) - |\sigma| \bar{\mathbf{v}}_K \cdot \tau_{K\sigma} \right) \\ &\quad - \sum_{\sigma=K|L \in \mathcal{F}_\Gamma} |\sigma| \mathbf{v}_{K\sigma} \cdot (\nu_K \mathbf{n}_{K\sigma}) + \sum_{\sigma=K|L \in \mathcal{F}_\Gamma} |\sigma| \tau_{K\sigma} \cdot (\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma} + \mathbf{v}_{K\sigma}). \end{aligned} \qquad (71) $$

By the normal continuity property embedded in the space $W$ we have $\tau_{K\sigma} = -\tau_{L\sigma}$ for all $\sigma = K|L \in \mathcal{F}^{\text{int}}$. Since $\bar{\mathbf{v}}_{K\sigma} = \bar{\mathbf{v}}_{L\sigma}$ whenever $\sigma \notin \mathcal{F}_\Gamma$, we infer

$$ \sum_{\sigma=K|L \in \mathcal{F}_\Gamma} |\sigma| \tau_{K\sigma} \cdot (\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma}) = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_{K\sigma} \cdot \tau_{K\sigma}. $$

Moreover, as $\sum_{\sigma \in \mathcal{F}_K} |\sigma| \mathbf{n}_{K\sigma} = 0$ for all $K \in \mathcal{M}$,

$$ \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_K \cdot (\nu_K \mathbf{n}_{K\sigma}) = \sum_{K \in \mathcal{M}} \bar{\mathbf{v}}_K \cdot (\nu_K \sum_{\sigma \in \mathcal{F}_K} |\sigma| \mathbf{n}_{K\sigma}) = 0. $$

Plugging these relations into (71) leads to

$$ \begin{aligned} w_D(\nu, \mathbf{v}_D) &= \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_{K\sigma} \cdot (-\nu_K \mathbf{n}_{K\sigma} + \tau_{K\sigma}) + \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_K \cdot (\nu_K \mathbf{n}_{K\sigma} - \tau_{K\sigma}) \\ &\quad + \sum_{\sigma=K|L \in \mathcal{F}_\Gamma} |\sigma| \tau_{K\sigma} \cdot \mathbf{v}_{K\sigma} - \sum_{\sigma=K|L \in \mathcal{F}_\Gamma} |\sigma| (\nu_K \mathbf{n}_{K\sigma}) \cdot \mathbf{v}_{K\sigma} \\ &= \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| (\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K) \cdot (\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}) + \sum_{\sigma=K|L \in \mathcal{F}_\Gamma} |\sigma| \mathbf{v}_{K\sigma} \cdot (\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}). \end{aligned} $$

Using the Cauchy–Schwarz inequality and invoking Lemma 5.11, we infer

$$ w_D(\nu, \mathbf{v}_D) \le \left( \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} \frac{|\sigma|}{h_K} |\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| h_K |\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2 \right)^{1/2} $$
$$
\begin{align}
&\lesssim |\sigma|^{1/2} \left( |K|^{-1/2} h_K^2 |\mathbf{u}|_{H^2(K)} + |L|^{-1/2} h_L^2 |\mathbf{u}|_{H^2(L)} \right) \nonumber \\
&\qquad + \|\pi_\sigma^0 [[\mathbf{u}] - [[\mathbf{u}]]]\|_{L^2(\sigma)} \nonumber \\
&\lesssim h_D^{3/2} \left( |\mathbf{u}|_{H^2(K)} + |\mathbf{u}|_{H^2(L)} \right) + h_\sigma |[[\mathbf{u}]]|_{H^1(\sigma)}, \tag{72}
\end{align}
$$

where we have introduced $\pm \pi_{\sigma}^{0} [[\mathbf{q}] - \mathbf{u}] = \pm ( [[\mathcal{I}_{U_{0,D}} \mathbf{q}}]]_{\sigma} - \pi_{\sigma}^{0} [[\mathbf{u}] ) )$ and used a triangle inequality in the first line, invoked (71) and used the $L^2(\sigma)$-boundedness of $\pi_{\sigma}^{0}$ together with (70) in the second inequality, and used the mesh regularity property in the conclusion to write $|\sigma|^{1/2} h_{K}^{1/2} \lesssim |K|^{1/2}$ and $|\sigma|^{1/2} h_{L}^{1/2} \lesssim |L|^{1/2}$, together with the approximation properties of $\pi_{\sigma}^{0}$ [16, Theorem 1.45]. Squaring, summing over $\sigma$ and taking the square root yields (67).

To prove the second estimate, we first notice that, by definition (7) of $[[\cdot]]_{\sigma}$ and (10) of the interpolator,

$$
[[\mathcal{I}_{U_{0,D}} \mathbf{u}]_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} (\gamma^{K\sigma} \mathbf{u} - \Pi^{L\sigma} (\mathcal{I}_{U_{0,D}} \mathbf{u})) = \pi_{\sigma}^{0} (\gamma^{K\sigma} \mathbf{u} - \Pi^{L\sigma} (\mathcal{I}_{U_{0,D}} \mathbf{u})).
$$

Hence,

$$
\begin{align*}
\| [[\mathcal{I}_{U_{0,D}} \mathbf{u}]_{\sigma} - \pi_{\sigma}^{0} [[\mathbf{u}]]]_{L^2(\sigma)} &= \| \pi_{\sigma}^{0} (\gamma^{K\sigma} \mathbf{u} - \Pi^{L\sigma} (\mathcal{I}_{U_{0,D}} \mathbf{u}) - [[\mathbf{u}]]) \|_{L^2(\sigma)} \\
&\le \| \gamma^{L\sigma} \mathbf{u} - \Pi^{L\sigma} (\mathcal{I}_{U_{0,D}} \mathbf{u}) \|_{L^2(\sigma)}.
\end{align*}
$$

Introduce as above a linear approximation $\mathbf{q}_L$ of $\mathbf{u}$ in $L$, and use the linear exactness (34) (with $L$ instead of $K$) of $\Pi^{L\sigma}$ to deduce

$$
\| [[ \mathcal{I}_{U_{0,D}} \mathbf{u} ] ]_{\sigma} - \pi_{\sigma}^{0} [[ \mathbf{u} ] ] \|_{L^2(\sigma)} \le \| \mathbf{u}|_L - \mathbf{q}_L \|_{L^2(\sigma)} + \| \Pi^{L\sigma} \mathcal{I}_{U_{0,D}} ( \mathbf{u}|_L - \mathbf{q}_L ) \|_{L^2(\sigma)}.
$$

In a similar way as (29) we can show that $\Pi^{L\sigma} \mathbf{w}_D$ is bounded above, on $\sigma$, by the maximum of the absolute values of the degrees of freedom associated with $\sigma$. By definition of $\mathcal{I}_{U_{0,D}}$ and (66) (with $L$ instead of $K$) and using similar arguments as for the first terms in (72), we infer

$$
\| [[ \mathcal{I}_{U_{0,D}} \mathbf{u} ] ]_{\sigma} - \pi_{\sigma}^{0} [[ \mathbf{u} ] ] \|_{L^2(\sigma)} \lesssim h_{L}^{3/2} |\mathbf{u}|_{H^2(L)}.
$$

Squaring, multiplying by $h_{\sigma}^{-1} \lesssim h_{L}^{-1}$, summing over $\sigma$ and taking the square root yields (68). $\square$

**Lemma 5.10 (Adjoint consistency).** *Recalling the definition (20) of the adjoint consistency error* $\mathcal{W}_D$, *it holds, for all* $\tau \in \Sigma$ *such that* $\tau \in H^1(\mathcal{M})^{d \times d}$,

$$
\mathcal{W}_D(\tau) \lesssim h_D |\tau|_{H^1(\mathcal{M})}. \tag{73}
$$

*Proof.* Let $\tau \in \Sigma \cap H^1(\mathcal{M})^{d \times d}$. The definitions (20) of $w_D$ and (7) of $[[\cdot]]_\sigma$ yield

$$
\begin{align*}
w_D(\tau, \mathbf{v}_D) &= \sum_{K \in \mathcal{M}} \left( -\left( \int_K \tau \right) : \mathcal{E}^K(\mathbf{v}_D) - \sum_{\sigma \in \mathcal{F}_K} \bar{\mathbf{v}}_K \cdot \int_{\sigma} (\tau|_K \mathbf{n}_{K\sigma}) \right) \\
&\quad + \sum_{\sigma=K|L \in \mathcal{F}_T} \left( \int_{\sigma} \tau|_K \mathbf{n}_{K\sigma} \right) \cdot \left( \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma} + \mathbf{v}_{K\sigma} \right) \\
&= \sum_{K \in \mathcal{M}} \left( -|K|_{\tau K} : \mathcal{E}^K(\mathbf{v}_D) - \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_K \cdot \boldsymbol{\tau}_{K\sigma} \right) \\
&\quad + \sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| \boldsymbol{\tau}_{K\sigma} \cdot \left( \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma} + \mathbf{v}_{K\sigma} \right),
\end{align*}
$$
$$ \begin{aligned} & + \left( \sum_{\sigma=K|L \in \mathcal{F}_T} \frac{|\sigma|}{h_K} |\mathbf{v}_{K\sigma}| \right)^{1/2} \left( \sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| h_K |\boldsymbol{\tau}_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2 \right)^{1/2} \\ & \lesssim \|\mathbf{v}_D\|_{1,D} \left( \sum_{K \in M} \sum_{\sigma \in \mathcal{F}_K} |\sigma| h_K |\boldsymbol{\tau}_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2 \right)^{1/2} . \end{aligned} \quad (72) $$

Note that the term $\sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| h_K |\boldsymbol{\tau}_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2$ has been included in the last factor in the right-hand side. By [17, Lemma B.6],

$$ |\sigma| h_K |\boldsymbol{\tau}_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2 \lesssim h_K^2 |\nu|^{1/2} H^1(K). $$

Plugging this into (72), dividing by $\|\mathbf{v}_D\|_{1,D}$ and taking the supremum over $\mathbf{v}_D$ concludes the proof. □

**Lemma 5.11.** For all $\mathbf{v}_D \in \mathbf{U}_{0,D}$, the following two inequalities hold:

$$ \left( \sum_{\sigma=K|L \in \mathcal{F}_T} \frac{|\sigma|}{h_K} |\mathbf{v}_{K\sigma}|^2 \right)^{1/2} \lesssim \|\mathbf{v}_D\|_{1,D}, \quad (73) $$

$$ \left( \sum_{K \in M} \sum_{\sigma \in \mathcal{F}_K} \frac{|\sigma|}{h_K} |\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K|^2 \right)^{1/2} \lesssim \|\mathbf{v}_D\|_{1,D}. \quad (74) $$

*Proof.* To prove (73), we simply write, by mesh regularity property, $\frac{|\sigma|}{h_K} \lesssim h_K^{d-2}$, so that, by definition (13) of the discrete norm,

$$ \sum_{\sigma=K|L \in \mathcal{F}_T} \frac{|\sigma|}{h_K} |\mathbf{v}_{K\sigma}|^2 \lesssim \sum_{\sigma=K|L \in \mathcal{F}_T} h_K^{d-2} |\mathbf{v}_{K\sigma}|^2 \lesssim S_D(\mathbf{v}_D, \mathbf{v}_D) \lesssim \|\mathbf{v}_D\|_{1,D}^2. $$

We now turn to (74). Let $K \in M$ and $\sigma \in \mathcal{F}_K$. By the choice (4) of the weights $(\omega_s^\sigma)_{s \in V_\sigma}$, the definition (5) of $\bar{\mathbf{v}}_{K\sigma}$ and the definition (8) of the linear function $\Pi^K \mathbf{v}_D$, we have

$$ \sum_{s \in V_\sigma} \omega_s^\sigma (\mathbf{v}_{Ks} - \Pi^K \mathbf{v}_D(\mathbf{x}_s)) = \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K - \nabla^K \mathbf{v}_D(\bar{\mathbf{x}}_\sigma - \bar{\mathbf{x}}_K). $$

The convexity of the function $x \to |x|^2$ then yields

$$ |\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K - \nabla^K \mathbf{v}_D(\bar{\mathbf{x}}_\sigma - \bar{\mathbf{x}}_K)|^2 \le \sum_{s \in V_\sigma} \omega_s^\sigma |\mathbf{v}_{Ks} - \Pi^K \mathbf{v}_D(\mathbf{x}_s)|^2. $$

Apply the inequality $\frac{1}{2}|\mathbf{a}|^2 \le |\mathbf{b}|^2 + |\mathbf{a} - \mathbf{b}|^2$ with $\mathbf{a} = \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K$ and $\mathbf{b} = \nabla^K \mathbf{v}_D(\bar{\mathbf{x}}_\sigma - \bar{\mathbf{x}}_K)$. Recalling the definition (11) of $S_K$, and using the facts that each $\omega_s^\sigma$ is in $[0, 1]$ and that the number of faces that meet at a vertex $s \in V_K$ is $\le 1$, we deduce that

$$ \frac{1}{2} |\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K|^2 \le h_K^2 |\nabla^K \mathbf{u}_D|^2 + h_K^{2-d} S_K(\mathbf{v}_D, \mathbf{v}_D). $$

The estimate (74) then follows by multiplying by $|\sigma|/h_K$, by noticing that $|\sigma| h_K \lesssim |K| \lesssim h_K^d$, and by summing over $K \in M$. □
with

$$ \nu_K = \frac{1}{|K|} \int_K \nu \quad \text{and} \quad \tau_{K\sigma} = \frac{1}{|\sigma|} \int_{\sigma} (\nu_K \mathbf{n}_{K\sigma}). $$

Recalling the definition (8) of $\nabla^K$ and noticing that $\nu_K : \mathcal{C}^K(\mathbf{v}_D) = \nu_K : \nabla^K \mathbf{v}_D$ since $\nu_K$ is symmetric, we infer

$$ w_D(\nu, \mathbf{v}_D) = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} \left( -|\sigma| \bar{\mathbf{v}}_{K\sigma} \cdot (\nu_K \mathbf{n}_{K\sigma}) - |\sigma| \bar{\mathbf{v}}_K \cdot \tau_{K\sigma} \right) \\ - \sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| \mathbf{v}_{K\sigma} \cdot (\nu_K \mathbf{n}_{K\sigma}) + \sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| \tau_{K\sigma} \cdot (\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma} + \mathbf{v}_{K\sigma}). \quad (74) $$

By the normal continuity property embedded in the space $\Sigma$ we have $\tau_{K\sigma} = -\tau_{L\sigma}$ for all $\sigma = K|L \in \mathcal{F}^{\text{in}}$. Since $\bar{\mathbf{v}}_{K\sigma} = \bar{\mathbf{v}}_{L\sigma}$ whenever $\sigma \notin \mathcal{F}_T$, we infer

$$ \sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| \tau_{K\sigma} \cdot (\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_{L\sigma}) = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_{K\sigma} \cdot \tau_{K\sigma}. $$

Moreover, as $\sum_{\sigma \in \mathcal{F}_K} |\sigma| \mathbf{n}_{K\sigma} = 0$ for all $K \in \mathcal{M}$,

$$ \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_K \cdot (\nu_K \mathbf{n}_{K\sigma}) = \sum_{K \in \mathcal{M}} \bar{\mathbf{v}}_K \cdot (\nu_K \sum_{\sigma \in \mathcal{F}_K} |\sigma| \mathbf{n}_{K\sigma}) = 0. $$

Plugging these relations into (74) leads to

$$ w_D(\nu, \mathbf{v}_D) = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_{K\sigma} \cdot (-\nu_K \mathbf{n}_{K\sigma} + \tau_{K\sigma}) + \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| \bar{\mathbf{v}}_K \cdot (\nu_K \mathbf{n}_{K\sigma} - \tau_{K\sigma}) \\ + \sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| \tau_{K\sigma} \cdot \mathbf{v}_{K\sigma} - \sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| (\nu_K \mathbf{n}_{K\sigma}) \cdot \mathbf{v}_{K\sigma} \\ = \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| (\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K) \cdot (\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}) + \sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| \mathbf{v}_{K\sigma} \cdot (\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}). $$

Using the Cauchy–Schwarz inequality and invoking Lemma 5.11 below, we infer

$$ w_D(\nu, \mathbf{v}_D) \le \left( \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} \frac{|\sigma|}{h_K} |\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K|^2 \right)^{1/2} \left( \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| h_K |\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2 \right)^{1/2} \\ + \left( \sum_{\sigma=K|L \in \mathcal{F}_T} \frac{|\sigma|}{h_K} |\mathbf{v}_{K\sigma}| \right)^{1/2} \left( \sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| h_K |\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2 \right)^{1/2} \\ \le \|\mathbf{v}_D\|_{1,\mathcal{D}} \left( \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} |\sigma| h_K |\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2 \right)^{1/2}. \quad (75) $$

Note that the term $\sum_{\sigma=K|L \in \mathcal{F}_T} |\sigma| h_K |\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2$ has been included in the last factor in the right-hand side. By [17, Lemma B.6],

$$ |\sigma| h_K |\tau_{K\sigma} - \nu_K \mathbf{n}_{K\sigma}|^2 \lesssim h_K^2 |\nu|_{H^1(K)}^2. $$

Plugging this into (75), dividing by $\|\mathbf{v}_D\|_{1,\mathcal{D}}$ and taking the supremum over $\mathbf{v}_D$ concludes the proof. $\square$
# 6 Numerical Experiments

## 6.1 Unbounded 2D domain with a single fracture under compression

This test case presented in [28, 21, 22, 18] for a Coulomb frictional contact model also applies here to the Tresca frictional contact model (1) since the normal traction of the analytical solution and the friction coefficient are constant along the fracture. It consists of a 2D unbounded domain containing a single fracture and subject to a compressive remote stress $\sigma = 100$ MPa. The fracture inclination with respect to the x-direction is $\psi = \pi/9$ and its length is $2\ell = 2$ m. The Coulomb-friction coefficient, Young's modulus and Poisson's ratio are set to $F = 1/\sqrt{3}$, $E = 25$ GPa and $\nu = 0.25$. The analytical solution in terms of the Lagrange multiplier $\lambda_n$ and of the jump of the tangential displacement field is given by:

$$ \lambda_n = \sigma \sin^2(\psi), \quad |\llbracket \mathbf{u} \rrbracket_{\tau}| = \frac{4(1-\nu)}{E} \sigma \sin(\psi) \left( \cos(\psi) - \frac{g}{\lambda_n} \sin(\psi) \right) \sqrt{\ell^2 - (\ell^2 - \tau^2)}, \quad (75) $$

where $0 \le \tau \le 2\ell$ is a curvilinear abscissa along the fracture. Note that, since $\lambda_n > 0$, we have $\llbracket \mathbf{u} \rrbracket_n = 0$ on the fracture. It results that the Tresca threshold is constant along the fracture and defined by $g = F\lambda_n$.

Boundary conditions are imposed on $\mathbf{u}$ at specific nodes of the mesh, as shown in Figure 4, to respect the symmetry of the expected solution. For this simulation, we sample a 320m x 320 m square, and carry out uniform refinements at each step in such a way to compute the solution on meshes containing 100, 200, 400, and 800 faces on the fracture (corresponding, respectively, to 12 468, 49 872, 199 488, and 797 952 triangular elements). The initial mesh is refined in a neighborhood of the fracture; starting from this mesh, we perform global uniform refinements at each step.

Figure 5 shows the comparison between the analytical and numerical Lagrange multipliers $\lambda_n$ and tangential displacement jump $\llbracket \mathbf{u} \rrbracket_{\tau}$ computed on the finest mesh. The Lagrange multiplier $\lambda_n$ presents some oscillations in a neighborhood of the fracture tips. As already explained in [21], this is due to the sliding of faces close to the fracture tips (in this test case, all fracture faces are in a contact-slip state). The discrete tangential displacement jump cannot be distinguished from the analytical solution on this fine mesh. Figure 6 displays the convergence of the tangential displacement jump and of the normal Lagrange multiplier as a function of the size of the largest fracture face denoted by $h$. Note that the $L^2$ error for the Lagrange multiplier is computed 5% away from each tip to circumvent the lack of convergence induced by the oscillations as in [21]. A first-order convergence for the displacement jump and a 1.5 convergence order for the Lagrange multiplier are observed. The former (low) rate is related to the low regularity of $\llbracket \mathbf{u} \rrbracket_{\tau}$ close to the tips (cf. the analytical expression (75)), the latter (higher than expected) rate is likely related to the fact that $\lambda_n$ is constant.
**Lemma 5.11.** For all $\mathbf{v}_{\mathcal{D}} \in \mathbf{U}_{0,\mathcal{D}}$, the following two inequalities hold:

$$ \left( \sum_{\sigma=K|L \in \mathcal{F}_T} \frac{|\sigma|}{h_K} |\mathbf{v}_{K\sigma}|^2 \right)^{1/2} \lesssim \|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}, \quad (76) $$

$$ \left( \sum_{K \in \mathcal{M}} \sum_{\sigma \in \mathcal{F}_K} \frac{|\sigma|}{h_K} |\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K|^2 \right)^{1/2} \lesssim \|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}. \quad (77) $$

*Proof.* To prove (76), we simply write $\frac{|\sigma|}{h_K} \lesssim h_K^{d-2}$ by mesh regularity property, so that, by definition (16) of the discrete norm,

$$ \sum_{\sigma=K|L \in \mathcal{F}_T} \frac{|\sigma|}{h_K} |\mathbf{v}_{K\sigma}|^2 \lesssim \sum_{\sigma=K|L \in \mathcal{F}_T} h_K^{d-2} |\mathbf{v}_{K\sigma}|^2 \lesssim S_{\mathcal{D}}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) \lesssim \|\mathbf{v}_{\mathcal{D}}\|_{1,\mathcal{D}}^2. $$

We now turn to (77). Let $K \in \mathcal{M}$ and $\sigma \in \mathcal{F}_K$. By the choice (5) of the weights $(\omega_s^\sigma)_{s \in \mathcal{V}_\sigma}$, the definition (6) of $\bar{\mathbf{v}}_{K\sigma}$ and the definition (9) of the linear function $\Pi^K \mathbf{v}_{\mathcal{D}}$, we have

$$ \sum_{s \in \mathcal{V}_\sigma} \omega_s^\sigma (\mathbf{v}_{Ks} - \Pi^K \mathbf{v}_\mathcal{D}(\mathbf{x}_s)) = \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K - \nabla^K \mathbf{v}_\mathcal{D}(\bar{\mathbf{x}}_\sigma - \bar{\mathbf{x}}_K). $$

The convexity of the function $x \to |x|^2$ then yields

$$ |\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K - \nabla^K \mathbf{v}_\mathcal{D}(\bar{\mathbf{x}}_\sigma - \bar{\mathbf{x}}_K)|^2 \le \sum_{s \in \mathcal{V}_\sigma} \omega_s^\sigma |\mathbf{v}_{Ks} - \Pi^K \mathbf{v}_\mathcal{D}(\mathbf{x}_s)|^2. $$

Apply the inequality $\frac{1}{2}|\mathbf{a}|^2 \le |\mathbf{b}|^2 + |\mathbf{a} - \mathbf{b}|^2$ with $\mathbf{a} = \bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K$ and $\mathbf{b} = \nabla^K \mathbf{v}_\mathcal{D}(\bar{\mathbf{x}}_\sigma - \bar{\mathbf{x}}_K)$. Recalling the definition (12) of $S_K$, and using the facts that each $\omega_s^\sigma$ is in $[0, 1]$ and that the number of faces that meet at a vertex $s \in \mathcal{V}_K$ is $\le 1$, we deduce that

$$ \frac{1}{2} |\bar{\mathbf{v}}_{K\sigma} - \bar{\mathbf{v}}_K|^2 \le h_K^2 |\nabla^K \mathbf{u}_\mathcal{D}|^2 + h_K^{2-d} S_K(\mathbf{v}_\mathcal{D}, \mathbf{v}_\mathcal{D}). $$

The estimate (77) then follows by multiplying by $|\sigma|/h_K$, by noticing that $|\sigma|h_K \lesssim |K| \lesssim h_K^d$, and by summing over $K \in \mathcal{M}$. $\square$

# 6 Numerical Experiments

## 6.1 Unbounded 2D domain with a single fracture under compression

This test case, presented in [18, 21, 22, 28] for a Coulomb frictional contact model, also applies here to the Tresca frictional contact model (2) since the normal traction of the analytical solution and the friction coefficient are constant along the fracture. It consists of a 2D unbounded domain containing a single fracture and subject to a compressive remote stress $\sigma = 100$ MPa. The fracture inclination with respect to the x-direction is $\psi = \pi/9$ and its length is $2\ell = 2$ m. The Coulomb-friction coefficient, the Young modulus and the Poisson ratio are set to $F = 1/\sqrt{3}$, $E = 25$ GPa and $\nu = 0.25$. The analytical solution in terms of the Lagrange multiplier $\lambda_n$ and of the jump of the tangential displacement field is given by:

$$ \lambda_n = \sigma \sin^2(\psi), \quad \|\|\mathbf{u}\|\|_\tau = \frac{4(1-\nu)}{E} \sigma \sin(\psi) \left( \cos(\psi) - \frac{g}{\lambda_n} \sin(\psi) \right) \sqrt{\ell^2 - (\ell^2 - \tau^2)}, \quad (78) $$
<figure><img src="image_5.png" /><figcaption>Figure 4. Unbounded domain containing a single fracture under uniform compression (a) and mesh including nodes for boundary conditions (♦: $u_x = 0$, ■: $u_y = 0$), for the example of Section 6.1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Comparison between the numerical and analytical solutions on the finest mesh (with 800 fracture faces), in terms of $\lambda_n$ (a) and $\|\mathbf{u}\|_{\tau}$ (b), example of Section 6.1.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Convergence of the relative $L^2$ error of $\|\mathbf{u}\|_{\tau} - \|\mathbf{u}_{\mathcal{D},\tau}\|_{\mathcal{D},\tau}$ and of $\lambda_n - \lambda_{\mathcal{D},n}$ away from the tip, as a function of the size of the largest fracture face denoted by $h$. Test case of Section 6.1.</figcaption></figure>
where $0 \le \tau \le 2\ell$ is a curvilinear abscissa along the fracture. Note that, since $\lambda_n > 0$, we have $[[u]]_n = 0$ on the fracture. It results that the Tresca threshold is constant along the fracture and defined by $g = F\lambda_n$.

Boundary conditions are imposed on **u** at specific nodes of the mesh, as shown in Figure 4, to respect the symmetry of the expected solution. For this simulation, we sample a 320m x 320 m square, and carry out uniform refinements at each step in such a way to compute the solution on meshes containing 100, 200, 400, and 800 faces on the fracture (corresponding, respectively, to 12 468, 49 872, 199 488, and 797 952 triangular elements). The initial mesh is refined in a neighbourhood of the fracture; starting from this mesh, we perform global uniform refinements at each step.

Figure 5 shows the comparison between the analytical and numerical Lagrange multipliers $\lambda_n$ and tangential displacement jump $[[u]]_\tau$ computed on the finest mesh. The Lagrange multiplier $\lambda_n$ presents some oscillations in a neighbourhood of the fracture tips. As already explained in [21], this is due to the sliding of faces close to the fracture tips (in this test case, all fracture faces are in a contact-slip state). The discrete tangential displacement jump cannot be distinguished from the analytical solution on this fine mesh. Figure 6 displays the convergence of the tangential displacement jump and of the normal Lagrange multiplier as a function of the size of the largest fracture face denoted by $h$. Note that the $L^2$ error for the Lagrange multiplier is computed 5% away from each tip to circumvent the lack of convergence induced by the oscillations as in [21]. A first-order convergence for the displacement jump and a 1.5 convergence order for the Lagrange multiplier are observed. The former (low) rate is related to the low regularity of $[[u]]_\tau$ close to the tips (cf. the analytical expression (78)), the latter (higher than expected) rate is likely related to the fact that $\lambda_n$ is constant.

<figure><img src="image_5.png" /><figcaption>Figure 4. (a) unbounded domain containing a single fracture under uniform compression, and (b) mesh including nodes for boundary conditions (♦: $u_x = 0$, ■: $u_y = 0$), for the example of Section 6.1.</figcaption></figure>