## 6.2 3D manufactured solution for the Tresca friction model

We consider the 3D domain $\Omega = (-1, 1)^3$ with the single non-immersed fracture $\Gamma = \{0\} \times (-1, 1)^2$. The Tresca coefficient $g$ is set to 1 and the Lamé coefficients are set to $\mu = \lambda = 1$. Note that this situation does not formally match the assumptions made in the analysis of the scheme: specifically, $\Omega \setminus \Gamma$ is not connected; however, this assumption was solely made to ensure that the Korn inequality is valid in the fractured domain, which is the case here since each connected component of $\Omega \setminus \Gamma$ sees $\partial \Omega$ (as a matter of fact, the connectedness assumption could be replaced by the assumption that each connected component of $\Omega \setminus \Gamma$ touches the boundary of $\Omega$ along some hypersurface).

The exact solution

$$ \mathbf{u}(x, y, z) = \begin{cases} \begin{pmatrix} h(x, y)P(z) - gy \\ P(z) \\ x^2P(z) \end{pmatrix} & \text{if } z \ge 0, \\ \begin{pmatrix} h(x, y)Q(z) - gy \\ 2Q(z) \\ x^2Q(z) \end{pmatrix} & \text{if } z < 0, x < 0, \\ \begin{pmatrix} h(x, y)Q(z) - gy \\ Q(z) \\ x^2Q(z) \end{pmatrix} & \text{if } z < 0, x \ge 0, \end{cases} \quad (76) $$

with $h(x, y) = -\sin(x)\cos(y)$, $P(z) = z^2$ and $Q(z) = z^2/4$, is designed to satisfy the Tresca frictional-contact conditions at the matrix fracture interface $\Gamma$. The right hand side $f = -\operatorname{div}\sigma(\mathbf{u})$ is deduced and the trace of $\mathbf{u}$ is imposed as Dirichlet boundary condition on $\partial\Omega$. Note that the fracture $\Gamma$ is in sticky-contact state for $z > 0$ ($\llbracket\mathbf{u}\rrbracket_n = 0$, $\llbracket\mathbf{u}\rrbracket_\tau = 0$) and slippy-contact for $z < 0$ ($\llbracket\mathbf{u}\rrbracket_n = 0$, $\llbracket\mathbf{u}\rrbracket_\tau > 0$). The convergence of the mixed $\mathbb{P}^1$-bubble VEM-$\mathbb{P}^0$ formulation is investigated on families of uniform Cartesian, tetrahedral, and hexahedral meshes. Starting from a uniform Cartesian mesh, a hexahedral mesh is generated by random perturbation of the nodes. This leads to non-planar faces (except on the fracture) which are dealt with either by cutting the faces into two triangles or by applying the modified gradient operator as described in remark 3.1. These two choices illustrated in Figure 7 are denoted respectively by Hexa-cut and Hexa-bary in the following.

<figure><img src="image_8.png" /><figcaption>Figure 7. Example of randomly perturbed Cartesian cell with non planar faces: hexahedral cell with a non planar face cut into two triangles (left, denoted by Hexa-cut) or four triangles using the isobarycenter of the face nodes (right, denoted by Hexa-bary). Note that in the Hexa-bary case, the displacement at the isobarycenter is eliminated by linear combination of the displacements at the face nodes as detailed in Remark 3.1.</figcaption></figure>

Figure 8 exhibits the relative $L^2$ norms of the errors $\mathbf{u} - \Pi_{\mathcal{D}}\mathbf{u}_{\mathcal{D}}$, $\llbracket\mathbf{u}\rrbracket - \llbracket\mathbf{u}_{\mathcal{D}}\rrbracket_{\mathcal{D}}$, $\nabla\mathbf{u} - \nabla_{\mathcal{D}}\mathbf{u}_{\mathcal{D}}$ and $\lambda_n - \lambda_{\mathcal{D},n}$ on the three families of refined meshes as a function of the cubic root of the number of
<figure><img src="image_6.png" /><figcaption>Figure 5. Comparison between the numerical and analytical solutions on the finest mesh (with 800 fracture faces), for (a) $\lambda_n$, and (b) $[\mathbf{u}]_{\tau}$; 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}_D]_{D,\tau}$ and of $\lambda_n - \lambda_{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>

## 6.2 3D manufactured solution for the Tresca friction model

We consider the 3D domain $\Omega = (-1, 1)^3$ with the single non-immersed fracture $\Gamma = \{0\} \times (-1, 1)^2$. The Tresca coefficient $g$ is set to 1, the Young modulus at $E = 2.5$ and the Poisson ratio at $\nu = 0.25$.
cells. It shows, as expected for such a smooth solution, a second-order convergence for **u** and [[**u**]] with all families of meshes. A first-order convergence, coherent with Theorem 4.6, is obtained for $\nabla u$ and $\lambda_n$ with both the hexahedral and tetrahedral families of meshes, while a second-order convergence for $\nabla u$ and a 1.5th-order convergence for $\lambda_n$ is observed with the family of Cartesian meshes (these improved rates being probably due to the symmetry and uniformity of the mesh).

<figure><img src="image_9.png" /><figcaption>Figure 8. Relative $L^2$ norms of the errors $\mathbf{u} - \Pi_{\mathcal{D}}\mathbf{u}_{\mathcal{D}}$, $[\mathbf{u}] - [\mathbf{u}_{\mathcal{D}}]_{\mathcal{D}}$, $\nabla \mathbf{u} - \nabla_{\mathcal{D}}\mathbf{u}_{\mathcal{D}}$ and $\lambda_n - \lambda_{\mathcal{D},n}$ as a function of the square root of the number of cells, using the families of Cartesian (a), tetrahedral (b), Hexa-cut (c) and Hexa-bary (d) meshes. Test case of Section 6.2.</figcaption></figure>

The discrete solution for the face-wise constant normal jump $[\mathbf{u}_{\mathcal{D}}]_{\mathcal{D},n}$ is essentially zero (machine precision of $10^{-16}$) for any mesh, as we are in a contact state on the fracture $\Gamma$. In contrast, the nodal normal jumps approach zero as the mesh is refined, as illustrated in Figure 9 for the family of Hexa-cut meshes. The nodal normal jumps in Figure 9 are plotted as a function of $z$ along the “broken” line corresponding to $x = y = 0$ before perturbation of the mesh. Figure 10 plots, for the Hexa-cut family of meshes, the face-wise constant non-zero tangential component jump $[\mathbf{u}_{\mathcal{D}}]_{\mathcal{D},y}$ on $\Gamma$ as well as the nodal tangential jumps as a function of $z$ as for $[\mathbf{u}_{\mathcal{D}}]_{\mathcal{D},n}$ in Figure 9. We recall that, on the fracture, the exact tangential jump $[\mathbf{u}]_y$ depends only on $z$ and is equal to $\min(z/2, 0)^2$.
The exact solution

$$ \mathbf{u}(x, y, z) = \begin{cases} \begin{pmatrix} h(x, y)P(z) - gy \\ P(z) \\ x^2P(z) \end{pmatrix} & \text{if } z \ge 0, \\ \begin{pmatrix} h(x, y)Q(z) - gy \\ 2Q(z) \\ x^2Q(z) \end{pmatrix} & \text{if } z < 0, x < 0, \\ \begin{pmatrix} h(x, y)Q(z) - gy \\ Q(z) \\ x^2Q(z) \end{pmatrix} & \text{if } z < 0, x \ge 0, \end{cases} \quad (79) $$

with $h(x, y) = -\sin(x)\cos(y)$, $P(z) = z^2$ and $Q(z) = z^2/4$, is designed to satisfy the Tresca frictional-contact conditions at the matrix fracture interface $\Gamma$. The right hand side $\mathbf{f} = -\mathbf{div}\sigma(\mathbf{u})$ is deduced and the trace of $\mathbf{u}$ is imposed as Dirichlet boundary condition on $\partial\Omega$. Note that the fracture $\Gamma$ is in sticky-contact state for $z > 0$ ($\llbracket\mathbf{u}\rrbracket_n = 0$, $\llbracket\mathbf{u}\rrbracket_\tau = 0$) and slippy-contact for $z < 0$ ($\llbracket\mathbf{u}\rrbracket_n = 0$, $|\llbracket\mathbf{u}\rrbracket_\tau| > 0$). The convergence of the polytopal scheme (11) is investigated on families of uniform Cartesian, tetrahedral, and hexahedral meshes. Starting from a uniform Cartesian mesh, an hexahedral mesh is generated by random perturbation of the nodes. This lead to non-planar faces (except on the fracture) which are handled either by cutting these faces into two triangles or by applying the modified gradient operator as described in Remark 3.1. These two choices, illustrated in Figure 7, are respectively denoted by Hexa-cut and Hexa-bary in the following.

<figure><img src="image_8.png" /><figcaption>Figure 7. Example of randomly perturbed Cartesian cell with non planar faces: hexahedral cell with a non planar face cut into two triangles (left, denoted by Hexa-cut) or four triangles using the isobarycenter of the face nodes (right, denoted by Hexa-bary). Note that in the Hexa-bary case, the displacement at the isobarycenter is eliminated by linear combination of the displacements at the face nodes as detailed in Remark 3.1.</figcaption></figure>

Figure 8 exhibits the relative $L^2$ norms of the errors $\mathbf{u} - \Pi_D\mathbf{u}_D$, $\llbracket\mathbf{u}\rrbracket - \llbracket\mathbf{u}_D\rrbracket_D$, $\nabla\mathbf{u} - \nabla_D\mathbf{u}_D$ and $\lambda_n - \lambda_{D,n}$ on the three families of refined meshes as a function of the cubic root of the number of cells. It shows, as expected for such a smooth solution, a second-order convergence for $\mathbf{u}$ and $\llbracket\mathbf{u}\rrbracket$ with all families of meshes. A first-order convergence, coherent with Theorem 4.6, is obtained for $\nabla\mathbf{u}$ and $\lambda_n$ with both the hexahedral and tetrahedral families of meshes, while a second-order convergence for $\nabla\mathbf{u}$ and a 1.5th-order convergence for $\lambda_n$ is observed with the family of Cartesian meshes (these improved rates being probably due to the symmetry and uniformity of this particular mesh).
<figure><img src="image_10.png" /><figcaption>Figure 9. Nodal normal jumps $\llbracket \mathbf{u}_D \rrbracket_{D,n}$ along the line $x = y = 0$ as a function of $z$ for the discrete solutions on the Hexa-cut meshes with $2^{3m}$ cells, $m = 3, 4, 5$ and for the continuous solution. Test case of Section 6.2.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. (a) Face-wise constant non-zero tangential component jump $\llbracket \mathbf{u}_D \rrbracket_{D,y}$ on $\Gamma$ obtained on the Hexa-cut mesh with $2^{3m}$ cells, $m = 5$. (b) Nodal tangential jumps along the line $x = y = 0$ as a function of $z$ both for the discrete solutions on the Hexa-cut meshes with $2^{3m}$ cells, $m = 3, 4, 5$ and for the continuous solution depending only on $z$. Test case of Section 6.2.</figcaption></figure>

# 7 Conclusions

We have presented in this work the convergence analysis of a fully discrete polytopal discretisation of a contact-mechanical model with Tresca friction at matrix-fracture interfaces. The analysis accounts for general elements and network of planar fractures including immersed and non-immersed fractures,
<figure><img src="image_9.png" /><figcaption>Figure 8. Relative $L^2$ norms of the errors $\mathbf{u} - \Pi_{\mathcal{D}}\mathbf{u}_{\mathcal{D}}$, $[\mathbf{u}] - [\mathbf{u}_{\mathcal{D}}]_{\mathcal{D}}$, $\nabla \mathbf{u} - \nabla_{\mathcal{D}}\mathbf{u}_{\mathcal{D}}$ and $\lambda_n - \lambda_{\mathcal{D},n}$ as a function of the cubic root of the number of cells, using the families of Cartesian (a), tetrahedral (b), Hexa-cut (c) and Hexa-bary (d) meshes. Test case of Section 6.2.</figcaption></figure>

The discrete solution for the face-wise constant normal jump $[\mathbf{u}_{\mathcal{D}}]_{\mathcal{D},n}$ is essentially zero (machine precision of $10^{-16}$) for any mesh, as we are in a contact state on the fracture $\Gamma$. In contrast, the nodal normal jumps approach zero as the mesh is refined, as illustrated in Figure 9 for the family of Hexa-cut meshes. In this figure, we plot the nodal normal jumps and the nodal tangential non-zero component jumps $[\mathbf{u}_{\mathcal{D}}]_{\mathcal{D},y}$ as functions of $z$ along the broken line corresponding to $x = y = 0$ before perturbation of the mesh. The tangential jumps are also presented in Figure 10 as face-wise constant values on $\Gamma$. We recall that, on the fracture, the exact tangential jump $[\mathbf{u}]_y$ depends only on $z$ and is equal to $\min(z/2, 0)^2$.
with intersections, corners and tips. The main new features are related to the polytopal bubble additional displacement degree of freedom, the proof of a discrete inf-sup condition using a specific $H_{0,j}^{-1/2}(\Gamma)$-norm to deal with networks of fractures and one-sided bubbles, and the proof of a discrete Korn inequality taking fracture networks into account. Numerical experiments based on two analytical solutions are presented and confirm the established first order error estimates on the $L^2(\Omega)$-norm of the displacement gradient and on the $H_{0,j}^{-1/2}(\Gamma)$-norm of the Lagrange multiplier.

## A Construction of $U_{\mathcal{K}_s(i)}$ and $\omega_{\mathcal{K}_s(i)}$ for the averaged interpolator

For each $K \in \mathcal{M}$ and $s \in \mathcal{V}_K$, we construct an explicit $U_{\mathcal{K}_s(i)}$ and $\omega_{\mathcal{K}_s(i)}$ that satisfies (21). The construction shows that $U_{\mathcal{K}_s(i)}$ can be entirely contained in any single element in $\mathcal{K}_s(i)$. For simplicity of presentation we assume here that the space dimension is $d = 3$, but the same construction can also be done in 2D.

Let $L \in \mathcal{K}_s(i)$. By mesh regularity there is a simplex $S \subset L$ that contains $s$ as one of its vertices, and that is shape-regular (with regularity factor bounded above by the global mesh regularity factor). We take $U_{\mathcal{K}_s(i)} = S$. We then build $\omega_{\mathcal{K}_s(i)}$ on the reference simplex

$$ \hat{S} = \text{co}\{(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)\} $$

and linearly transport it onto $S = U_{\mathcal{K}_s(i)}$; the relation (21b) are preserved by this linear transport, and the upper bound in (21a) is the same as the one on the reference simplex.

Without loss of generality we can therefore assume that $\mathbf{x}_s = (0, 0, 0)$. The simplex $\hat{S}$ has center of mass $\mathbf{x}_{\hat{S}} = (1/3, 1/3, 1/3)$, and the simplex $S_1 = \frac{1}{2}\hat{S}$ has center of mass $\mathbf{x}_{S_1} = (1/6, 1/6, 1/6)$. We therefore have

$$ \mathbf{x}_s = 2\mathbf{x}_{S_1} - \mathbf{x}_{\hat{S}}. \qquad (77) $$

Let us define $\omega_{\mathcal{K}_s(i)} = 16\mathbf{1}_{S_1} - \mathbf{1}_{\hat{S}}$, where $\mathbf{1}_A$ is the characteristic function of A. Then, since $|S_1| = |\hat{S}|/8$,

$$ \int_{\hat{S}} \omega_{\mathcal{K}_s(i)} = 16|S_1| - |\hat{S}| = |\hat{S}|, $$

which establishes the first relation in (21b). To prove the second, by definition of the centers of mass we write

$$ \int_{\hat{S}} \mathbf{x} \omega_{\mathcal{K}_s(i)} = \int_{S_1} 16\mathbf{x} - \int_{\hat{S}} \mathbf{x} = 16|S_1|\mathbf{x}_{S_1} - |\hat{S}|\mathbf{x}_{\hat{S}} = |\hat{S}|(2\mathbf{x}_{S_1} - \mathbf{x}_{\hat{S}}) = |\hat{S}|\mathbf{x}_s, $$

where the conclusion follows from (77).

## B Estimate on fracture norms

We establish here two estimates involving the $H_{0,j}^{-1/2}(\Gamma)$-like norms, that are used in the error estimates. The first one is an approximation property of the $L^2$-projector $\pi_{\mathcal{F}_T}^0$ on piecewise constant functions.

**Lemma B.1** (Approximation properties of $\pi_{\mathcal{F}_T}^0$ in $L^2$ and discrete $H_{0,j}^{-1/2}(\Gamma)$ norms). We have

$$ \|\lambda - \pi_{\mathcal{F}_T}^0 \lambda\|_{L^2(\Gamma)} \le h_D |\lambda|_{H^1(\mathcal{F}_T)} \quad \forall \lambda \in H^1(\mathcal{F}_T)^d, \qquad (78) $$

and, for all $s \in [0, 1]$,

$$ \|\lambda - \pi_{\mathcal{F}_T}^0 \lambda\|_{-1/2, \Gamma} \le h_D^{\frac{1}{2}+s} |\lambda|_{H^s(\mathcal{F}_T)} \quad \forall \lambda \in H^s(\mathcal{F}_T)^d. \qquad (79) $$
<figure><img src="image_10.png" /><figcaption>Figure 9. (a) Nodal normal jumps $\llbracket \mathbf{u}_D \rrbracket_{D,n}$ and (b) nodal tangential jumps $\llbracket \mathbf{u}_D \rrbracket_{D,y}$, along the line $x = y = 0$ as functions of $z$ for the discrete solutions on the Hexa-cut meshes with $2^{3m}$ cells, $m = 3, 4, 5$. Test case of Section 6.2.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. Face-wise constant non-zero tangential component jump $\llbracket \mathbf{u}_D \rrbracket_{D,y}$ on $\Gamma$ obtained on the Hexa-cut mesh with $2^{3m}$ cells, $m = 5$. Test case of Section 6.2.</figcaption></figure>

# 7 Conclusions

We presented in this work the convergence analysis of a fully discrete polytopal discretisation of a contact-mechanical model with Tresca friction at matrix-fracture interfaces, using a mixed model with the displacement and a Lagrange multiplier (representing the surface traction) as unknowns. The analysis accounts for general mesh elements and a network of planar fractures – including immersed and non-immersed fractures, with intersections, corners and tips. The main new features are related to the additional polytopal bubble displacement degree of freedom, the proof of a discrete inf-sup condition using a suitable dual norm on the Lagrange multiplier (to deal with networks of fractures and one-sided bubble degrees of freedom), and the proof of a discrete Korn inequality taking into account the fracture network. Numerical experiments based on two analytical solutions are presented
Proof. Let $\lambda \in H^s(\mathcal{F}_\Gamma)^d$. The approximation properties of $L^2$-orthogonal projectors on polynomial spaces give, for all $\sigma \in \mathcal{F}_\Gamma$ and denoting by $\pi_\sigma^0$ the $L^2$-orthogonal projector on $\mathbb{P}^0(\sigma)^d$,

$$ \|\lambda - \pi_{\sigma}^{0} \lambda\|_{L^2(\sigma)} \lesssim h_{\mathcal{D}}^{s} |\lambda|_{H^{s}(\sigma)} $$

(this estimate is established in [16, Theorem 1.45] for the scalar case and $s \in \{0, 1\}$, and can be obtained for vector-valued functions and $s \in (0, 1)$ working component by component and interpolating between $s = 0$ and $s = 1$). Squaring, summing over $\sigma \in \mathcal{F}_\Gamma$ and taking the square root, we infer

$$ \|\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda\|_{L^2(\Gamma)} \lesssim h_{\mathcal{D}}^s |\lambda|_{H^s(\mathcal{F}_\Gamma)} \quad (80) $$

The case $s = 1$ corresponds to (78).

We now turn to (79). Given the definition (14), we only need to bound the norm on $\Gamma_i$ for each $i \in I$. For all $\mathbf{v}_i \in H^1(\Omega_i^+; \Gamma_i) \setminus \{0\}$, by orthogonality property of $\pi_{\mathcal{F}_\Gamma}^0$ and a Cauchy–Schwarz inequality we have

$$ \begin{aligned} \int_{\Gamma_i} (\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda) \cdot \mathbf{v}_i &= \int_{\Gamma_i} (\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda) \cdot (\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i) \\ &\le \|\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda\|_{L^2(\Gamma_i)} \|\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i\|_{L^2(\Gamma_i)} \\ &\lesssim h_{\mathcal{D}}^s |\lambda|_{H^s(\mathcal{F}_\Gamma)} \|\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i\|_{L^2(\Gamma_i)}, \end{aligned} \quad (81) $$

where we have used (80) in the second inequality. Letting $\pi_M^0$ be the $L^2$-orthogonal projector on piecewise constant functions over $M$ and introducing $(\pi_M^0 \mathbf{v}_i)|_{\Gamma_i}$ we have, since $\pi_{\mathcal{F}_\Gamma}^0 (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i} = (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i}$,

$$ \begin{aligned} \|\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i\|_{L^2(\Gamma_i)} &\le \|\mathbf{v}_i - (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i}\|_{L^2(\Gamma_i)} + \|\pi_{\mathcal{F}_\Gamma}^0 (\mathbf{v}_i - (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i})\|_{L^2(\Gamma_i)} \\ &\le 2\|\mathbf{v}_i - (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i}\|_{L^2(\Gamma_i)}, \end{aligned} $$

where the conclusion comes from the boundedness in $L^2(\Gamma_i)$-norm of $\pi_{\mathcal{F}_\Gamma}^0$. The approximation properties of $\pi_M^0$ (see [16, Theorem 1.45]) then yield

$$ \|\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i\|_{L^2(\Gamma_i)} \lesssim h_{\mathcal{D}}^{\frac{1}{2}} \|\mathbf{v}_i\|_{H^1(\Omega_i^+)} \quad (82) $$

Plugging (82) into (81), dividing by $\|\mathbf{v}_i\|_{H^1(\Omega_i^+)}$ and taking the supremum over $\mathbf{v}_i$ concludes the proof. $\square$

The second property is a bound between $H^{-1/2}$-like norms on $\Gamma$.

**Lemma B.2 (Estimate of the discrete $H^{-1/2}$-norm).** It holds

$$ \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}} \lesssim \|\lambda_{\mathcal{D}}\|_{-1/2, \Gamma}, \quad \forall \lambda_{\mathcal{D}} \in \mathbf{M}_{\mathcal{D}}. $$

Proof. Let $\lambda \in \mathbf{M}_{\mathcal{D}}$. Given the definitions (14) and (57) of the two norms, we only need to prove that, for each $i \in I$, there exists $\mathbf{v}_i \in H^1(\Omega_i^+; \Gamma_i)^d$ such that

$$ \|\mathbf{v}_i\|_{H^1(\Omega_i^+)} \lesssim \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}, i} \quad \text{and} \quad \int_{\Gamma_i} \lambda_{\mathcal{D}} \cdot \mathbf{v}_i = \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}, i}^2 \quad (83) $$
and confirm the established first order error estimates on the $L^2(\Omega)$-norm of the displacement gradient and on the dual norm of the Lagrange multiplier.

## A Construction of $U_{\mathcal{K}_s(i)}$ and $\varpi_{\mathcal{K}_s(i)}$ for the averaged interpolator

For each $K \in \mathcal{M}$ and $s \in \mathcal{V}_K$, we construct an explicit $U_{\mathcal{K}_s(i)}$ and $\varpi_{\mathcal{K}_s(i)}$ that satisfies (24). The construction shows that $U_{\mathcal{K}_s(i)}$ can be entirely contained in any single element in $\mathcal{K}_s(i)$. For simplicity of presentation we assume here that the space dimension is $d = 3$, but the same construction can also be done in 2D.

Let $L \in \mathcal{K}_s(i)$. By mesh regularity there is a simplex $S \subset L$ that contains $s$ as one of its vertices, and that is shape-regular (with regularity factor bounded above by the global mesh regularity factor), and has a measure comparable to the measure of $L$. We fix $U_{\mathcal{K}_s(i)} = S$, then build $\varpi_{\mathcal{K}_s(i)}$ on the reference simplex

$$ \hat{S} = \text{co}\{(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)\} $$

and linearly transport it onto $S$; the relation (24b) are preserved by this linear transport, and the upper bound in (24a) is the same as the one on the reference simplex.

Without loss of generality we can therefore assume that $\mathbf{x}_s = (0, 0, 0)$. The simplex $\hat{S}$ has center of mass $\mathbf{x}_{\hat{S}} = (1/3, 1/3, 1/3)$, and the simplex $S_1 = \frac{1}{2}\hat{S}$ has center of mass $\mathbf{x}_{S_1} = (1/6, 1/6, 1/6)$, so that

$$ \mathbf{x}_s = 2\mathbf{x}_{S_1} - \mathbf{x}_{\hat{S}}. \qquad (80) $$

Let us define $\varpi_{\mathcal{K}_s(i)} = 16\mathbf{1}_{S_1} - \mathbf{1}_{\hat{S}}$, where $\mathbf{1}_A$ is the characteristic function of $A$. Since $|S_1| = |\hat{S}|/8$ we have

$$ \int_{\hat{S}} \varpi_{\mathcal{K}_s(i)} = 16|S_1| - |\hat{S}| = |\hat{S}|, $$

which establishes the first relation in (24b). To prove the second, by definition of the centers of mass we write

$$ \int_{\hat{S}} \mathbf{x} \varpi_{\mathcal{K}_s(i)} = \int_{S_1} 16\mathbf{x} - \int_{\hat{S}} \mathbf{x} = 16|S_1|\mathbf{x}_{S_1} - |\hat{S}|\mathbf{x}_{\hat{S}} = |\hat{S}|(2\mathbf{x}_{S_1} - \mathbf{x}_{\hat{S}}) = |\hat{S}|\mathbf{x}_s, $$

where the conclusion follows from (80).

## B Estimate on fracture norms

We establish here two bounds involving the dual norms on $\Gamma$, that are used in the proof of the error estimates. The first one is an approximation property of the $L^2$-projector $\pi_{\mathcal{F}_\Gamma}^0$ on piecewise constant functions.

**Lemma B.1** (Approximation properties of $\pi_{\mathcal{F}_\Gamma}^0$ in $L^2$ and $H^{-1/2}$-like norms). We have

$$ \|\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda\|_{L^2(\Gamma)} \lesssim h_{\mathcal{D}} |\lambda|_{H^1(\mathcal{F}_\Gamma)} \quad \forall \lambda \in H^1(\mathcal{F}_\Gamma), \qquad (81) $$

and, for all $s \in [0, 1]$,

$$ \|\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda\|_{-1/2, \Gamma} \lesssim h_{\mathcal{D}}^{1+s} |\lambda|_{H^s(\mathcal{F}_\Gamma)} \quad \forall \lambda \in H^s(\mathcal{F}_\Gamma). \qquad (82) $$
where

$$ \|\lambda_D\|_{-1/2, D, i} = \left( \sum_{\sigma \subset \Gamma_i} h_\sigma |\sigma| |\lambda_\sigma|^2 \right)^{1/2} $$

is the restriction to $\Gamma_i$ of $\|\lambda_D\|_{-1/2, D}$. To achieve this, we first define the boundary values of $v_i$ and then lift these boundary values to create $v_i$ itself.

### Step 1: design of bubble functions.

For each $\sigma \subset \Gamma_i$ we take a ball $B_\sigma \subset \sigma$ and a bubble function $b_\sigma \in W^{1,\infty}(\sigma)$ such that

$$ \begin{gathered} B_\sigma \text{ has radius } \geq h_\sigma, \quad \text{dist}(B_\sigma, \partial\sigma) \geq h_\sigma, \\ \int_\sigma b_\sigma = |\sigma|, \quad b_\sigma = 0 \text{ outside } B_\sigma, \quad \|b_\sigma\|_{L^\infty} \leq 1, \quad \|\nabla b_\sigma\|_{L^\infty} \leq h_\sigma^{-1}. \end{gathered} \quad (84) $$

The existence of $B_\sigma$ is ensured by the mesh regularity assumption, and $b_\sigma$ can then be constructed on $B_\sigma$ by scaling the function $x \to \text{dist}(x, \partial B_\sigma)$. We then have

$$ |b\sigma|_{H^{1/2}(\sigma)}^2 = \int_{\sigma} \int_{\sigma} \frac{|b_{\sigma}(\mathbf{x}) - b_{\sigma}(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^d} dxdy \le \|\nabla b_{\sigma}\|_{L^{\infty}}^2 \int_{\sigma} \int_{\sigma} \frac{dxdy}{|\mathbf{x} - \mathbf{y}|^{d-2}}. $$

Using polar coordinates around $x$ (setting $y = x + \rho\xi$ with $\rho > 0$ and $\xi$ of unit length) and recalling that $\sigma$ has dimension $d-1$ we easily get $\int_{\sigma} \frac{dy}{|x-y|^{d-2}} \le h_{\sigma}$. Hence, invoking (84), we infer

$$ |b\sigma|_{H^{1/2}(\sigma)}^2 \le h_{\sigma}^{-1}|\sigma|. \quad (85) $$

For any $x \in B_\sigma$, since $x$ is at distance $\ge h_\sigma$ of $\partial\sigma$ a use of polar coordinates around $x$ yields $\int_{\partial\Omega_i^+ \setminus \sigma} \frac{dy}{|x-y|^d} \le h_\sigma^{-1}$. Using (84), we infer that

$$ \int_{\sigma} \int_{\partial \Omega_i^+ \setminus \sigma} \frac{|b_{\sigma}(\mathbf{x})|^2}{|\mathbf{x} - \mathbf{y}|^d} \le h_{\sigma}^{-1} |\sigma|. \quad (86) $$

### Step 2: construction of the boundary value.

Set $w : \partial\Omega_i^+ \to \mathbb{R}^d$ such that $w|_{\sigma} = h_\sigma \lambda_\sigma b_\sigma$ for all $\sigma \subset \Gamma_i$, and $w|_{\partial\Omega_i^+ \setminus \Gamma_i} = 0$. Then,

$$ \int_{\Gamma_i} \lambda_D \cdot w = \sum_{\sigma \subset \Gamma_i} h_\sigma |\lambda_\sigma|^2 \int_{\sigma} b_\sigma = \sum_{\sigma \subset \Gamma_i} h_\sigma |\lambda_\sigma|^2 |\sigma| = \|\lambda_D\|_{-1/2, D, i}^2. \quad (87) $$

By the support condition in (84), $w$ is continuous (and piecewise $H^1$) on $\partial\Omega_i^+$, so it belongs to $H^{1/2}(\partial\Omega_i^+)^d$. Having in mind to create a lifting of $w$, we want to bound its $H^{1/2}$ norm. Specifically, we will prove that

$$ \|w\|_{H^{1/2}(\partial\Omega_i^+)} \le \|\lambda_D\|_{-1/2, D, i}. \quad (88) $$

We only detail the bound of the seminorm, as bounding the $L^2$ norm of $w$ is straightforward. We have

$$ \begin{aligned} |w|_{H^{1/2}(\partial\Omega_i^+)}^2 &= \sum_{\sigma \subset \partial\Omega_i^+} \sum_{\sigma' \subset \partial\Omega_i^+} \int_{\sigma} \int_{\sigma'} \frac{|w(\mathbf{x}) - w(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^d} dxdy \\ &= \sum_{\sigma \subset \Gamma_i} \int_{\sigma} \int_{\sigma} \frac{|w(\mathbf{x}) - w(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^d} dxdy + \sum_{\sigma \subset \partial\Omega_i^+} \sum_{\sigma' \subset \partial\Omega_i^+, \sigma' \ne \sigma} \int_{\sigma} \int_{\sigma'} \frac{|w(\mathbf{x}) - w(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^d} dxdy \end{aligned} $$
Proof. Let $\lambda \in H^s(\mathcal{F}_\Gamma)$. The approximation properties of $L^2$-orthogonal projectors on polynomial spaces give, for all $s \in [0, 1]$ and $\sigma \in \mathcal{F}_\Gamma$, denoting by $\pi_\sigma^0$ the $L^2$-orthogonal projector on $\mathbb{P}^0(\sigma)^d$,

$$ \|\lambda - \pi_{\sigma}^{0} \lambda\|_{L^2(\sigma)} \lesssim h_{\mathcal{D}}^{s} |\lambda|_{H^{s}(\sigma)} $$

(this estimate is established in [16, Theorem 1.45] for the scalar case and $s \in \{0, 1\}$, and can be obtained for vector-valued functions and $s \in (0, 1)$ working component by component and interpolating between $s = 0$ and $s = 1$). Squaring, summing over $\sigma \in \mathcal{F}_\Gamma$ and taking the square root, we infer

$$ \|\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda\|_{L^2(\Gamma)} \lesssim h_{\mathcal{D}}^s |\lambda|_{H^s(\mathcal{F}_\Gamma)}. \quad (83) $$

The case $s = 1$ corresponds to (81).

We now turn to (82). Given the definition (17), we only need to bound the norm on $\Gamma_i$ for each $i \in I$. For all $\mathbf{v}_i \in H_0^1(\Omega_i^+; \Gamma_i) \setminus \{0\}$, by orthogonality property of $\pi_{\mathcal{F}_\Gamma}^0$ and a Cauchy–Schwarz inequality we have

$$ \begin{aligned} \int_{\Gamma_i} (\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda) \cdot \mathbf{v}_i &= \int_{\Gamma_i} (\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda) \cdot (\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i) \\ &\le \|\lambda - \pi_{\mathcal{F}_\Gamma}^0 \lambda\|_{L^2(\Gamma_i)} \|\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i\|_{L^2(\Gamma_i)} \\ &\lesssim h_{\mathcal{D}}^s |\lambda|_{H^s(\mathcal{F}_\Gamma)} \|\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i\|_{L^2(\Gamma_i)}, \end{aligned} \quad (84) $$

where we have used (83) in the second inequality. Letting $\pi_M^0$ be the $L^2$-orthogonal projector on piecewise constant functions over $M$ and introducing $(\pi_M^0 \mathbf{v}_i)|_{\Gamma_i}$ we have, since $\pi_{\mathcal{F}_\Gamma}^0 (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i} = (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i}$,

$$ \begin{aligned} \|\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i\|_{L^2(\Gamma_i)} &\le \|\mathbf{v}_i - (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i}\|_{L^2(\Gamma_i)} + \|\pi_{\mathcal{F}_\Gamma}^0 (\mathbf{v}_i - (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i})\|_{L^2(\Gamma_i)} \\ &\le 2\|\mathbf{v}_i - (\pi_M^0 \mathbf{v}_i)|_{\Gamma_i}\|_{L^2(\Gamma_i)}, \end{aligned} $$

where the conclusion comes from the boundedness in $L^2(\Gamma_i)$-norm of $\pi_{\mathcal{F}_\Gamma}^0$. The approximation properties of $\pi_M^0$ (see [16, Theorem 1.45]) then yield

$$ \|\mathbf{v}_i - \pi_{\mathcal{F}_\Gamma}^0 \mathbf{v}_i\|_{L^2(\Gamma_i)} \lesssim h_{\mathcal{D}}^{\frac{1}{2}} \|\mathbf{v}_i\|_{H^1(\Omega_i^+)} \quad (85) $$

Plugging (85) into (84), dividing by $\|\mathbf{v}_i\|_{H^1(\Omega_i^+)}$ and taking the supremum over $\mathbf{v}_i$ concludes the proof. $\square$

The second property is a bound between the discrete $H^{-1/2}$ and the $H^{-1/2}$-like norms on $\Gamma$.

**Lemma B.2 (Estimate of dual norms on $\Gamma$).** Recalling the definitions (17) of the $H^{-1/2}$-like norm and (60) of the discrete $H^{-1/2}$-norm on $\Gamma$, it holds

$$ \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}} \lesssim \|\lambda_{\mathcal{D}}\|_{-1/2, \Gamma}, \quad \forall \lambda_{\mathcal{D}} \in \mathbf{M}_{\mathcal{D}}. $$

Proof. Let $\lambda \in \mathbf{M}_{\mathcal{D}}$. Given the definitions (17) and (60) of the two norms, we only need to prove that, for each $i \in I$, there exists $\mathbf{v}_i \in H_0^1(\Omega_i^+; \Gamma_i)$ such that

$$ \|\mathbf{v}_i\|_{H^1(\Omega_i^+)} \lesssim \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}, i} \quad \text{and} \quad \int_{\Gamma_i} \lambda_{\mathcal{D}} \cdot \mathbf{v}_i = \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}, i}^2 \quad (86) $$
$$ = \mathfrak{I}_1 + \mathfrak{I}_2. \qquad (89) $$

To bound $\mathfrak{I}_1$, we use $\mathbf{w}|_{\sigma} = h_{\sigma}\lambda_{\sigma}b_{\sigma}$ and (85):

$$ \mathfrak{I}_1 = \sum_{\sigma \subset \Gamma_i} |\mathbf{w}|_{\sigma}|_{H^{1/2}(\sigma)}^2 \lesssim \sum_{\sigma \subset \Gamma_i} h_{\sigma}|\sigma| |\lambda_{\sigma}|^2 = \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}, i}^2. \qquad (90) $$

We now turn to $\mathfrak{I}_2$, for which we write

$$ \begin{aligned} \mathfrak{I}_2 &\le 2 \sum_{\sigma \subset \partial \Omega_i^+} \sum_{\sigma' \subset \partial \Omega_i^+, \sigma' \ne \sigma} \int_{\sigma} \int_{\sigma'} \frac{|\mathbf{w}(\mathbf{x})|^2}{|\mathbf{x}-\mathbf{y}|^d} dxdy + 2 \sum_{\sigma \subset \partial \Omega_i^+} \sum_{\sigma' \subset \partial \Omega_i^+, \sigma' \ne \sigma} \int_{\sigma} \int_{\sigma'} \frac{|\mathbf{w}(\mathbf{y})|^2}{|\mathbf{x}-\mathbf{y}|^d} dxdy \\ &= 4 \sum_{\sigma \subset \partial \Omega_i^+} \int_{\sigma} \int_{\partial \Omega_i^+ \setminus \sigma} \frac{|\mathbf{w}(\mathbf{x})|^2}{|\mathbf{x}-\mathbf{y}|^d} dxdy \\ &\lesssim \sum_{\sigma \subset \partial \Omega_i^+} h_{\sigma} |\sigma| |\lambda_{\sigma}|^2 = \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}, i}^2, \end{aligned} \qquad (91) $$

where the second line follows gathering the two terms in the right-hand side of the first line by symmetric roles of $\sigma$ and $\sigma'$, and the last line from $\mathbf{w}|_{\sigma} = h_{\sigma}\lambda_{\sigma}b_{\sigma}$ and (86).

Plugging (90) and (91) into (89) concludes the proof of (88).

*Step 3: conclusion.*

Since $\mathbf{w} \in H^{1/2}(\partial \Omega_i^+)^d$ we can find a lifting $\mathbf{v}_i \in H^1(\Omega_i^+)^d$ of $\mathbf{w}$ such that $\|\mathbf{v}_i\|_{H^1(\Omega_i^+)} \lesssim \|\mathbf{w}\|_{H^{1/2}(\partial \Omega_i^+)}$. Since $\mathbf{w}$ vanishes outside $\Gamma_i$, we have $\mathbf{v}_i \in H^1(\Omega_i^+; \Gamma_i)$. By (88) and (87), the relations (83) hold and the proof is complete. $\square$

## Acknowledgement

We acknowledge the partial support from the joint laboratory IFPEN-Inria Convergence HPC/AI/HPDA for the energetic transition in realizing this project. We would like to thank Isabelle Faille and Guillaume Enchéry from IFPEN for fruitful discussions and comments during the elaboration of this work. Jérôme Droniou additionally acknowledges the funding of the European Union (ERC Synergy, NEMESIS, project number 101115663). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

## References

[1] Odd Andersen, Halvor M. Nilsen, and Xavier Raynaud. Virtual element method for geomechanical simulations of reservoir models. *Computational Geosciences*, 21(5):877–893, 2017.

[2] Laurence Beaude, Franz Chouly, Mohamed Laaziri, and Roland Masson. Mixed and nitsche's discretizations of coulomb frictional contact-mechanics for mixed dimensional poromechanical models. *Computer Methods in Applied Mechanics and Engineering*, 413:116124, 2023.

[3] L. Beirão Da Veiga, F. Brezzi, and L.D. Marini. Virtual elements for linear elasticity problems. *SIAM Journal on Numerical Analysis*, 51:794–812, 2013.
where

$$ \|\lambda_D\|_{-1/2, \mathcal{D}, i} := \left( \sum_{\sigma \subset \Gamma_i} h_\sigma |\sigma| |\lambda_\sigma|^2 \right)^{1/2} $$

is the restriction to $\Gamma_i$ of $\|\lambda_D\|_{-1/2, \mathcal{D}}$. To achieve this, we first define the boundary values of $v_i$ and then lift these boundary values to create $v_i$ itself.

Step 1: design of bubble functions.

For each $\sigma \subset \Gamma_i$ we take a $(d-1)$-dimensional ball $\mathcal{B}_\sigma \subset \sigma$ and a bubble function $b_\sigma \in W^{1,\infty}(\sigma)$ such that

$$ \mathcal{B}_\sigma \text{ has radius } \ge h_\sigma, \quad \text{dist}(\mathcal{B}_\sigma, \partial\sigma) \ge h_\sigma, \\ \int_\sigma b_\sigma = |\sigma|, \quad b_\sigma = 0 \text{ outside } \mathcal{B}_\sigma, \quad \|b_\sigma\|_{L^\infty} \le 1, \quad \|\nabla b_\sigma\|_{L^\infty} \le h_\sigma^{-1}. \qquad (87) $$

The existence of $\mathcal{B}_\sigma$ is ensured by the mesh regularity assumption, and $b_\sigma$ can then be constructed on $\mathcal{B}_\sigma$ by scaling the function $\mathbf{x} \to \text{dist}(\mathbf{x}, \partial\mathcal{B}_\sigma)$. We then have

$$ |b_\sigma|^2_{H^{1/2}(\sigma)} = \int_\sigma \int_\sigma \frac{|b_\sigma(\mathbf{x}) - b_\sigma(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^d} d\mathbf{x}d\mathbf{y} \le \|\nabla b_\sigma\|_{L^\infty}^2 \int_\sigma \int_\sigma \frac{d\mathbf{x}d\mathbf{y}}{|\mathbf{x} - \mathbf{y}|^{d-2}}. $$

Using polar coordinates around $\mathbf{x}$ (setting $\mathbf{y} = \mathbf{x} + \rho\xi$ with $\rho > 0$ and $\xi$ of unit length) and recalling that $\sigma$ has dimension $d-1$ we easily get $\int_\sigma \frac{dy}{|\mathbf{x}-\mathbf{y}|^{d-2}} \le h_\sigma$. Hence, invoking (87), we infer

$$ |b_\sigma|^2_{H^{1/2}(\sigma)} \le h_\sigma^{-1} |\sigma|. \qquad (88) $$

For any $\mathbf{x} \in \mathcal{B}_\sigma$, since $\mathbf{x}$ is at distance $\ge h_\sigma$ of $\partial\sigma$ a use of polar coordinates around $\mathbf{x}$ yields $\int_{\partial\Omega_i^+ \setminus \sigma} \frac{dy}{|\mathbf{x}-\mathbf{y}|^d} \le h_\sigma^{-1}$. Using (87), we infer that

$$ \int_{\sigma} \int_{\partial \Omega_i^+ \setminus \sigma} \frac{|b_{\sigma}(\mathbf{x})|^2}{|\mathbf{x} - \mathbf{y}|^d} \le h_{\sigma}^{-1} |\sigma|. \qquad (89) $$

Step 2: construction of the boundary value.

Set $\mathbf{w} : \partial\Omega_i^+ \to \mathbb{R}^d$ such that $\mathbf{w}|_\sigma = h_\sigma\lambda_\sigma b_\sigma$ for all $\sigma \subset \Gamma_i$, and $\mathbf{w}|_{\partial\Omega_i^+ \setminus \Gamma_i} = 0$. Then,

$$ \int_{\Gamma_i} \lambda_D \cdot \mathbf{w} = \sum_{\sigma \subset \Gamma_i} h_\sigma |\lambda_\sigma|^2 \int_\sigma b_\sigma = \sum_{\sigma \subset \Gamma_i} h_\sigma |\lambda_\sigma|^2 |\sigma| = \|\lambda_D\|_{-1/2, \mathcal{D}, i}^2 \qquad (90) $$

By the support condition in (87), $\mathbf{w}$ belongs to $H^1 \subset H^{1/2}(\partial\Omega_i^+)$. Having in mind to create a lifting of $\mathbf{w}$, we want to bound its $H^{1/2}$-norm. Specifically, we will prove that

$$ \|\mathbf{w}\|_{H^{1/2}(\partial\Omega_i^+) \to \mathcal{D}} \le \|\lambda_D\|_{-1/2, \mathcal{D}, i}. \qquad (91) $$

We only detail the bound of the seminorm, as bounding the $L^2$-norm of $\mathbf{w}$ is straightforward. We have

$$ |\mathbf{w}|_{H^{1/2}(\partial\Omega_i^+)}^2 = \sum_{\sigma \subset \partial\Omega_i^+} \sum_{\sigma' \subset \partial\Omega_i^+} \int_{\sigma} \int_{\sigma'} \frac{|\mathbf{w}(\mathbf{x}) - \mathbf{w}(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^d} d\mathbf{x}d\mathbf{y} $$
[4] F. Ben Belgacem and Y. Renard. Hybrid finite element methods for the signorini problem. *Mathematics of Computation*, 72:1117–1145, 2003.

[5] R. L. Berge, I. Berre, E. Keilegavlen, J. M. Nordbotten, and B. Wohlmuth. Finite volume discretization for poroelastic media with fractures modeled by contact mechanics. *International Journal for Numerical Methods in Engineering*, 121:644–663, 2019.

[6] F. Bonaldi, K. Brenner, J. Droniou, and R. Masson. Gradient discretization of two-phase flows coupled with mechanical deformation in fractured porous media. *Computers and Mathematics with Applications*, 98:40–68, 2021.

[7] F. Bonaldi, J. Droniou, R. Masson, and A. Pasteau. Energy-stable discretization of two-phase flows in deformable porous media with frictional contact at matrix-fracture interfaces. *Journal of Computational Physics*, 455:Paper No. 110984, 28, 2022.

[8] Bonaldi, Francesco, Brenner, Konstantin, Droniou, Jérôme, Masson, Roland, Pasteau, Antoine, and Trenty, Laurent. Gradient discretization of two-phase poro-mechanical models with discontinuous pressures at matrix fracture interfaces. *ESAIM: M2AN*, 55(5):1741–1777, 2021.

[9] W. M. Boon and J. M. Nordbotten. Mixed-dimensional poromechanical models of fractured porous media. *Acta Mechanica*, 2022.

[10] A. Borio, F. Hamon, N. Castelletto, J.A. White, and R.S. Settgast. Hybrid mimetic finite-difference and virtual element formulation for coupled poromechanics. *Computer Methods in Applied Mechanics and Engineering*, 383:113917, 2021.

[11] H. Chi, L. Beirão da Veiga, and G.H. Paulino. Some basic formulations of the virtual element method (vem) for finite deformations. *Computer Methods in Applied Mechanics and Engineering*, 318:148–192, 2017.

[12] Franz Chouly, Alexandre Ern, and Nicolas Pignet. A hybrid high-order discretization combined with nitsche's method for contact and tresca friction in small strain elasticity. *SIAM Journal on Scientific Computing*, 42(4):A2300–A2324, 2020.

[13] Julien Coulet, Isabelle Faille, Vivette Girault, Nicolas Guy, and Frédéric Nataf. A fully coupled scheme using virtual element method and finite volume for poroelasticity. *Computational Geosciences*, 24:381–403, 2020.

[14] L Beirão Da Veiga, Carlo Lovadina, and David Mora. A virtual element method for elastic and inelastic problems on polytope meshes. *Computer methods in applied mechanics and engineering*, 295:327–346, 2015.

[15] D. Di Pietro and S. Lemaire. An extension of the Crouzeix-Raviart space to general meshes with application to quasi-incompressible linear elasticity and Stokes flow. *Mathematics of Computation*, 84:1–31, 2015.

[16] Daniele Antonio Di Pietro and Jérôme Droniou. *The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications*, volume 19 of *Modeling, Simulation and Applications*. Springer International Publishing, 2020.

[17] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. *The Gradient Discretisation Method*, volume 82 of *Mathematics & Applications*. Springer, 2018.

[18] Jérôme Droniou, Guillaume Enchéry, Isabelle Faille, Ali Haidar, and Roland Masson. A bubble vem–fully discrete polytopal scheme for mixed-dimensional poromechanics with frictional con-
$$ 
\begin{aligned}
&= \sum_{\sigma \subset \Gamma_i} \int_{\sigma} \int_{\sigma} \frac{|\mathbf{w}(\mathbf{x}) - \mathbf{w}(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^d} dxdy + \sum_{\sigma \subset \partial \Omega_i^+} \sum_{\sigma' \subset \partial \Omega_i^+, \sigma' \neq \sigma} \int_{\sigma} \int_{\sigma'} \frac{|\mathbf{w}(\mathbf{x}) - \mathbf{w}(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^d} dxdy \\
&= \mathfrak{T}_1 + \mathfrak{T}_2.
\end{aligned} 
\quad (92) $$

To bound $\mathfrak{T}_1$, we use $\mathbf{w}|_{\sigma} = h_{\sigma} \lambda_{\sigma} b_{\sigma}$ and (88):

$$ \mathfrak{T}_1 = \sum_{\sigma \subset \Gamma_i} |\mathbf{w}|_{\sigma}|_{H^{1/2}(\sigma)}^2 \lesssim \sum_{\sigma \subset \Gamma_i} h_{\sigma} |\sigma| |\lambda_{\sigma}|^2 = \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}, i}^2 \quad (93) $$

We now turn to $\mathfrak{T}_2$, for which we write

$$ 
\begin{aligned}
\mathfrak{T}_2 &\le 2 \sum_{\sigma \subset \partial \Omega_i^+} \sum_{\sigma' \subset \partial \Omega_i^+, \sigma' \ne \sigma} \int_{\sigma} \int_{\sigma'} \frac{|\mathbf{w}(\mathbf{x})|^2}{|\mathbf{x} - \mathbf{y}|^d} dxdy + 2 \sum_{\sigma \subset \partial \Omega_i^+} \sum_{\sigma' \subset \partial \Omega_i^+, \sigma' \ne \sigma} \int_{\sigma} \int_{\sigma'} \frac{|\mathbf{w}(\mathbf{y})|^2}{|\mathbf{x} - \mathbf{y}|^d} dxdy \\
&= 4 \sum_{\sigma \subset \partial \Omega_i^+} \int_{\sigma} \int_{\partial \Omega_i^+ \setminus \sigma} \frac{|\mathbf{w}(\mathbf{x})|^2}{|\mathbf{x} - \mathbf{y}|^d} dxdy \\
&\lesssim \sum_{\sigma \subset \partial \Omega_i^+} h_{\sigma} |\sigma| |\lambda_{\sigma}|^2 = \|\lambda_{\mathcal{D}}\|_{-1/2, \mathcal{D}, i}^2,
\end{aligned} 
\quad (94) $$

where the second line follows gathering the two terms in the right-hand side of the first line by symmetric roles of $\sigma$ and $\sigma'$, and the last line from $\mathbf{w}|_{\sigma} = h_{\sigma} \lambda_{\sigma} b_{\sigma}$ and (89).

Plugging (93) and (94) into (92) concludes the proof of (91).

### Step 3: conclusion.

Since $\mathbf{w} \in H^{1/2}(\partial \Omega_i^+)$ we can find a lifting $\mathbf{v}_i \in H^1(\Omega_i^+)$ of $\mathbf{w}$ such that $\|\mathbf{v}_i\|_{H^1(\Omega_i^+)} \lesssim \|\mathbf{w}\|_{H^{1/2}(\partial \Omega_i^+)}$. Since $\mathbf{w}$ vanishes outside $\Gamma_i$, we have $\mathbf{v}_i \in H_0^1(\Omega_i^+; \Gamma_i)$. By (91) and (90), the relations (86) hold and the proof is complete. $\square$

## Acknowledgement

We acknowledge the partial support from the joint laboratory IFPEN-Inria Convergence HPC/AI/HPDA for the energetic transition which fully funded the postdoctoral fellowship of Ali Haidar. We would like to thank Isabelle Faille and Guillaume Enchéry from IFPEN for fruitful discussions and comments during the elaboration of this work.

This research pas partially funded by the European Union (ERC Synergy, NEMESIS, project number 101115663). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

## References

[1] O. Andersen, H. M. Nilsen, and X. Raynaud. “Virtual element method for geomechanical simulations of reservoir models”. In: *Computational Geosciences* 21.5 (2017), pp. 877–893.

[2] L. Beaude, F. Chouly, M. Laaziri, and R. Masson. “Mixed and Nitsche’s discretizations of Coulomb frictional contact-mechanics for mixed dimensional poromechanical models”. In: *Computer Methods in Applied Mechanics and Engineering* 413 (2023), p. 116124.
tact at matrix-fracture interfaces. *Computer Methods in Applied Mechanics and Engineering*, 422:116838, 2024.

[19] Guillaume Drouet and Patrick Hild. An accurate local average contact method for nonmatching meshes. *Numerische Mathematik*, 136(2):467–502, 2017.

[20] Guillaume Enchéry and Léo Agélas. Coupling linear virtual element and non-linear finite volume methods for poroelasticity. *Comptes Rendus. Mécanique*, 2023. Online first.

[21] A. Franceschini, N. Castelletto, J.A. White, and H.A. Tchelepi. Algebraically stabilized Lagrange multiplier method for frictional contact mechanics with hydraulically active fractures. *Computer Methods in Applied Mechanics and Engineering*, 368:113161, 2020.

[22] T. T. Garipov, M. Karimi-Fard, and H.A. Tchelepi. Discrete fracture model for coupled flow and geomechanics. *Computational Geosciences*, 20(1):149–160, 2016.

[23] T.T. Garipov and M.H. Hui. Discrete fracture modeling approach for simulating coupled thermo-hydro-mechanical effects in fractured reservoirs. *International Journal of Rock Mechanics and Mining Sciences*, 122:104075, 2019.

[24] P. Hansbo and M.G. Larson. Discontinuous Galerkin and the Crouzeix–Raviart element: Application to elasticity. *ESAIM: Mathematical Modelling and Numerical Analysis*, 37:63–72, 2003.

[25] Haslinger J., I. Hlaváček, and J. Necas. Numerical methods for unilateral problems in solid mechanics. In *Finite Element Methods (Part 2)*, Numerical Methods for Solids (Part 2), volume 4 of *Handbook of Numerical Analysis*, pages 313–485. Elsevier, 1996.

[26] E. Keilegavlen and J.M. Nordbotten. Finite volume methods for elasticity with weak symmetry. *International Journal for Numerical Methods in Engineering*, 112:939–962, 2017.

[27] Morteza Nejati, Adriana Paluszny, and Robert W. Zimmerman. A finite element framework for modeling internal frictional contact in three-dimensional fractured media using unstructured tetrahedral meshes. *Computer Methods in Applied Mechanics and Engineering*, 306:123–150, 2016.

[28] A.V. Phan, J.A.L. Napier, L.J. Gray, and T. Kaplan. Symmetric-Galerkin BEM simulation of fracture with frictional contact. *International Journal for Numerical Methods in Engineering*, 57:835–851, 2003.

[29] Ivar Stefansson, Inga Berre, and Eirik Keilegavlen. A fully coupled numerical model of thermo-hydro-mechanical processes and fracture contact mechanics in porous media. *Computer Methods in Applied Mechanics and Engineering*, 386:114122, 2021.

[30] P. Wriggers, W. T. Rust, and B. D. Reddy. A virtual element method for contact. *Computational Mechanics*, 58(6):1039–1050, 2016.

[31] Peter Wriggers, Fadi Aldakheel, and Blaž Hudobivnik. *Virtual Elements for Fracture Processes*, pages 243–315. Springer International Publishing, Cham, 2024.
[3] L. Beirão Da Veiga, F. Brezzi, and L. Marini. "VIRTUAL ELEMENTS FOR LINEAR ELASTICITY PROBLEMS". In: *SIAM Journal on Numerical Analysis* 51 (2013), pp. 794–812.

[4] F. Ben Belgacem and Y. Renard. "Hybrid Finite Element Methods for the Signorini problem". In: *Mathematics of Computation* 72 (2003), pp. 1117–1145.

[5] R. L. Berge, I. Berre, E. Keilegavlen, J. M. Nordbotten, and B. Wohlmuth. "Finite volume discretization for poroelastic media with fractures modeled by contact mechanics". In: *International Journal for Numerical Methods in Engineering* 121 (2019), pp. 644–663. doi: 10.1002/nme.6238.

[6] F. Bonaldi, K. Brenner, J. Droniou, and R. Masson. "Gradient discretization of two-phase flows coupled with mechanical deformation in fractured porous media". In: *Computers and Mathematics with Applications* 98 (2021), pp. 40–68. URL: https://arxiv.org/abs/2004.09860.

[7] F. Bonaldi, J. Droniou, R. Masson, and A. Pasteau. "Energy-stable discretization of two-phase flows in deformable porous media with frictional contact at matrix-fracture interfaces". In: *Journal of Computational Physics* 455 (2022), Paper No. 110984, 28. ISSN: 0021-9991. DOI: 10.1016/j.jcp.2022.110984. URL: https://doi.org/10.1016/j.jcp.2022.110984.

[8] F. Bonaldi, K. Brenner, J. Droniou, R. Masson, A. Pasteau, and L. Trenty. "Gradient discretization of two-phase poro-mechanical models with discontinuous pressures at matrix fracture interfaces". In: *ESAIM: M2AN* 55.5 (2021), pp. 1741–1777. DOI: 10.1051/m2an/2021036. URL: https://doi.org/10.1051/m2an/2021036.

[9] W. M. Boon and J. M. Nordbotten. "Mixed-dimensional poromechanical models of fractured porous media". In: *Acta Mechanica* (2022). DOI: 10.1007/s00707-022-03378-1. URL: https://doi.org/10.1007/s00707-022-03378-1.

[10] A. Borio, F. Hamon, N. Castelletto, J. White, and R. Settgast. "Hybrid mimetic finite-difference and virtual element formulation for coupled poromechanics". In: *Computer Methods in Applied Mechanics and Engineering* 383 (2021), p. 113917. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2021.113917. URL: https://www.sciencedirect.com/science/article/pii/S0045782516309094.

[11] H. Chi, L. B. da Veiga, and G. Paulino. "Some basic formulations of the virtual element method (VEM) for finite deformations". In: *Computer Methods in Applied Mechanics and Engineering* 318 (2017), pp. 148–192. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2016.12.020. URL: https://www.sciencedirect.com/science/article/pii/S0045782516309094.

[12] F. Chouly, A. Ern, and N. Pignet. "A Hybrid High-Order Discretization Combined with Nitsche's Method for Contact and Tresca Friction in Small Strain Elasticity". In: *SIAM Journal on Scientific Computing* 42.4 (2020), A2300–A2324.

[13] J. Coulet, I. Faille, V. Girault, N. Guy, and F. Nataf. "A fully coupled scheme using virtual element method and finite volume for poroelasticity". In: *Computational Geosciences* 24 (2020), pp. 381–403.

[14] L. B. Da Veiga, C. Lovadina, and D. Mora. "A virtual element method for elastic and inelastic problems on polytope meshes". In: *Computer methods in applied mechanics and engineering* 295 (2015), pp. 327–346.

[15] D. Di Pietro and S. Lemaire. "An extension of the Crouzeix-Raviart space to general meshes with application to quasi-incompressible linear elasticity and Stokes flow". In: *Mathematics of Computation* 84 (2015), pp. 1–31.
[16] D. A. Di Pietro and J. Droniou. *The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications*. Vol. 19. Modeling, Simulation and Applications. Springer International Publishing, 2020, xxxi+525p. ISBN: 978-3-030-37202-6 (Hardcover), 978-3-030-37203-3 (eBook). DOI: 10.1007/978-3-030-37203-3. URL: https://hal.archives-ouvertes.fr/hal-021518

[17] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. *The Gradient Discretisation Method*. Vol. 82. Mathematics & Applications. Springer, 2018, 511p. ISBN: 978-3-319-79041-1 (Softcover) 978-3-319-79042-8 (eBook). DOI: 10.1007/978-3-319-79042-8.

[18] J. Droniou, G. Enchéry, I. Faille, A. Haidar, and R. Masson. "A bubble VEM-fully discrete polytopal scheme for mixed-dimensional poromechanics with frictional contact at matrix-fracture interfaces". In: *Computer Methods in Applied Mechanics and Engineering* 422 (2024), p. 116838. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2024.116838. URL: https://www.sciencedirect.com/science/article/pii/S004578252400094X.

[19] G. Drouet and P. Hild. "An accurate local average contact method for nonmatching meshes". In: *Numerische Mathematik* 136.2 (2017), pp. 467-502.

[20] G. Enchéry and L. Agélas. "Coupling linear virtual element and non-linear finite volume methods for poroelasticity". en. In: *Comptes Rendus. Mécanique* (2023). Online first. DOI: 10.5802/crmeca.225.

[21] A. Franceschini, N. Castelletto, J. White, and H. Tchelepi. "Algebraically stabilized Lagrange multiplier method for frictional contact mechanics with hydraulically active fractures". In: *Computer Methods in Applied Mechanics and Engineering* 368 (2020), p. 113161.

[22] T. T. Garipov, M. Karimi-Fard, and H. Tchelepi. "Discrete fracture model for coupled flow and geomechanics". In: *Computational Geosciences* 20.1 (2016), pp. 149-160. DOI: 10.1007/s10596-015-9554-z. URL: https://doi.org/10.1007/s10596-015-9554-z.

[23] T. Garipov and M. Hui. "Discrete Fracture Modeling approach for simulating coupled thermo-hydro-mechanical effects in fractured reservoirs". In: *International Journal of Rock Mechanics and Mining Sciences* 122 (2019), p. 104075. ISSN: 1365-1609. DOI: https://doi.org/10.1016/j.ijrmms.2019.100010. URL: https://www.sciencedirect.com/science/article/pii/S1365160918311018.

[24] P. Hansbo and M. Larson. "DISCONTINUOUS GALERKIN AND THE CROUZEIX-RAVIART ELEMENT: APPLICATION TO ELASTICITY". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 37 (2003), pp. 63-72.

[25] H. J., I. Hlaváček, and J. Necas. "Numerical methods for unilateral problems in solid mechanics". In: *Finite Element Methods (Part 2), Numerical Methods for Solids (Part 2)*. Vol. 4. Handbook of Numerical Analysis. Elsevier, 1996, pp. 313-485. DOI: https://doi.org/10.1016/S1570-8659(96) URL: https://www.sciencedirect.com/science/article/pii/S1570865996800056.

[26] E. Keilegavlen and J. Nordbotten. "Finite volume methods for elasticity with weak symmetry". In: *International Journal for Numerical Methods in Engineering* 112 (2017), pp. 939-962.

[27] M. Nejati, A. Paluszny, and R. W. Zimmerman. "A finite element framework for modeling internal frictional contact in three-dimensional fractured media using unstructured tetrahedral meshes". In: *Computer Methods in Applied Mechanics and Engineering* 306 (2016), pp. 123-150. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2016.03.028. URL: https://www.sciencedirect.com/science/article/pii/S004578251630113X.

[28] A. Phan, J. Napier, L. Gray, and T. Kaplan. "Symmetric-Galerkin BEM simulation of fracture with frictional contact". In: *International Journal for Numerical Methods in Engineering* 57 (2003), pp. 835-851.
[29] I. Stefansson, I. Berre, and E. Keilegavlen. "A fully coupled numerical model of thermo-hydro-mechanical processes and fracture contact mechanics in porous media". In: *Computer Methods in Applied Mechanics and Engineering* 386 (2021), p. 114122. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2021.114122. URL: https://www.sciencedirect.com/science/a

[30] P. Wriggers, W. T. Rust, and B. D. Reddy. "A virtual element method for contact". In: *Computational Mechanics* 58.6 (2016), pp. 1039–1050.

[31] P. Wriggers, F. Aldakheel, and B. Hudobivnik. "Virtual Elements for Fracture Processes". In: *Virtual Element Methods in Engineering Sciences*. Cham: Springer International Publishing, 2024, pp. 243–315. ISBN: 978-3-031-39255-9. DOI: 10.1007/978-3-031-39255-9_10. URL: https://doi.org/10.1007/978-3-031-39255-9_10.