Dual virtual element method for discrete fractures networksA. Fumagalli, E.
Keilegavlen

# Dual virtual element method for discrete fractures networks††thanks:
Submitted to the editors DATE. This work was partially funded by the ANIGMA
project from the Research Council of Norway (project no. 244129/E20) through
the ENERGIX program.

Alessio Fumagalli Department of Mathematics, University of Bergen, Bergen,
Norway<EMAIL_ADDRESS><EMAIL_ADDRESS>Eirik
Keilegavlen22footnotemark: 2

###### Abstract

Discrete fracture networks is a key ingredient in the simulation of physical
processes which involve fluid flow in the underground, when the surrounding
rock matrix is considered impervious. In this paper we present two different
models to compute the pressure field and Darcy velocity in the system. The
first allows a normal flow out of a fracture at the intersections, while the
second grants also a tangential flow along the intersections. For the
numerical discretization, we use the mixed virtual element method as it is
known to handle grid elements of, almost, any arbitrary shape. The flexibility
of the discretization allows us to loosen the requirements on grid
construction, and thus significantly simplify the flow discretization compared
to traditional discrete fracture network models. A coarsening algorithm, from
the algebraic multigrid literature, is also considered to further speed up the
computation. The performance of the method is validated by numerical
experiments.

###### keywords:

Porous media, Discrete fracture network, Interface model, Virtual element
method.

76S05, 65N08, 65N30

##  1 Introduction

In many hard rocks fractures, naturally occurring or engineered, are of a
paramount importance to understand and simulate flow paths. The construction
of efficient simulation models for flow in fracture networks is therefore of
relevance for applications such as energy recovery and storage, waste disposal
(nuclear and CO2), just to name a few [4, 29, 25].

Flow in the fracture planes is usually modeled by Darcy’s law [32, 29], or
Forchheimer’s law if the Reynold number is sufficiently high [21, 30], and it
may also be of crucial interest to include flow along intersections in the
models [40, 20, 22]. Because of hydraulic aperture, which is several order of
magnitude smaller than other characteristic sizes of the problem, from a
modeling perspective a fracture may be represented as a two-dimensional plane
(generally manifold) [1, 32, 2], embedded in a three-dimensional domain. The
set of, possibly intersecting, fractures forms the network where the flow may
take place and, usually, is referred as a discrete fracture network (DFN).
Depending on how the fractures were created, they may intersect in an
arbitrary manner, leading to highly complex simulation geometries. The flow
properties of a fracture are determined by its geo-mechanical and geo-chemical
history [18], and there can be significant heterogeneities as well as
permeability anisotropy in the fracture plane. Infilling processes or
geological movements may alter also the composition and orientation of the
material presents in the intersections, leaving the latter as a privileged
patterns or obstacles for the flow.

One of the main challenge in DFN simulations in complex fracture geometries is
the construction of the computational mesh, combined with the subsequent
discretization of the flow equations. Depending on the properties of the
chosen numerical method, this may put significant constraints on the mesh
algorithm, in particular in meshing of fracture intersections [33, 26, 27].
Compared to two-dimensional gridding problems, which by themselves can be non-
trivial, gridding of DFN problems is complicated by the requirement, commonly
made, that the grids in different fracture planes are mutually conforming.
That is, hanging nodes are not allowed in fracture intersections. The gridding
can be simplified by allowing for a relaxed interpretation of the fracture
planes [28], in effect locally assigning a curvature to the fracture plane.

If this requirement is loosened, and some of the burden is transferred to an
appropriately chosen numerical method, the gridding problem becomes
significantly simpler and standard algorithms may be applied to create the
global mesh. For instance, the extended finite element method handles
discontinuities internal to elements by enriching the approximation space [24,
16, 22]. This removes all requirements on the mesh conforming to fracture
intersections, and thus significantly simplifies the gridding [22, 11, 19,
35].

A reasonable compromise between burdening mesh generation and discretization
is to require that meshes conform to fracture intersections, but also allow
for hanging nodes. This approach allows for the independent construction of a
set of bi-dimensional meshes in the fracture planes, which is considerably
less difficult than gridding the whole DFN model fully coupled. See Fig. 1 and
Fig. 2 as an example.

Figure 1: Example of pressure field in a realistic DFN with several
intersecting fractures, represented as ellipses. The fractures are represented
by one co-dimensional objects. The blue line represents an injection well,
while the red line a production well. Note that the represented mesh is
conforming at the intersections. A detailed representation of the grid is
reported in Fig. 2.

Figure 2: Four different details of the grid for the example represented in
Fig. 1. The method considered is able to handle these intersection among
fractures as conforming.

To further increase the applicability of discretization scheme, it should then
allow for general polygonal cells, including hanging nodes. In particular, we
note that the virtual element method (VEM) has successfully been applied in
the DFN setting [10, 9]. See [8, 6, 14, 5, 3, 7] for a presentation of VEM in
various contexts.

In this work, we consider different models to describe the single-phase flow
in the network of fractures, with and without the tangential contribution of
the intersections included. A reduced model is thus provided, with suitable
normal and tangential effective permeabilities, to describe the flow in a two-
codimensional framework. We consider a dual formulation of the generalized
Darcy problem where the virtual element method is extend to handle the
proposed modes. Local mass conservation is thus guaranteed and the normal
flux, from each edge in the mesh, may be directly used to simulate the heat or
tracer transport in the DFN. A gridding strategy is also presented to generate
a conforming mesh, in the sense that the edges may be split to ensure
conformity along the intersections of the fractures. The proposed scheme is
thus a flexible and robust tool to perform simulations in this context. We
present several numerical tests to validate the approach, focusing the study
on the order of convergence for the pressure and the velocity in different
scenarios.

The paper is organized as follows. The mathematical models and their analyses
to describe the fluid flow in the DFN are presented in Section 2.The Section 3
is devoted to prove the well posedness of the continuous models considered. In
Section 4 we review the implementation of the mixed virtual element method,
and its implementation in a DFN setting. The gridding and coarsening strategy
are presented in Section 5. Experimental results presented in Section 6 mainly
consider the error decay of the numerical method presented. Conclusions follow
in Section 7. Appendix A briefly describes the coarsening strategy while
Appendix B contains tables related to the examples.

##  2 Continuous model

This section is devoted to the presentation of the physical model. In Section
2.1 the setting for single-phase flow is introduced. Section 2.2 presents the
model for a single fracture, while the different couplings among fractures are
described in Section 2.3.

###  2.1 Single-phase flow in DFN

We make use of the symbol $\sharp(A)\in\mathbb{N}$ indicating the counting
measure of $A$. We define $\mathcal{N}$ the set of indexes associated to the
network of fractures, one value identify a fracture. Let us consider
$\sharp(\mathcal{N})$ distinct and planar domains $\Omega_{i}$, for
$i\in\mathcal{N}$, embedded in $\mathbb{R}^{3}$ such that their union is
indicated by $\Omega\vcentcolon=\cup_{i\in\mathcal{N}}\Omega_{i}$. Each
polygon represents a fracture and $\Omega$ is the discrete fracture network.
We indicate by $\mathcal{I}$ the set of indexes associated to the intersection
among fractures, one value identify one intersection. We define also the set
of one co-dimension intersections among fractures as $\gamma$, which can be
viewed also as a disjoint union of $\sharp(\mathcal{I})$ lines $\gamma_{k}$
such that $\gamma=\cup_{k\in\mathcal{I}}\gamma_{k}$. We indicate, for each
fracture, the finite set of indexes
$\mathcal{G}_{i}\vcentcolon=\left\\{k\in\mathbb{N}:\gamma_{k}\cap\Omega_{i}\neq\emptyset\right\\}$.
We suppose that if $k\in\mathcal{G}_{i}$ then $\gamma_{k}$ belongs to the
internal part of $\Omega_{i}$, that is, we assume fractures do not intersect
along fracture boundaries. We indicate by $\mathcal{P}$ the set of indexes
associated with the intersection of one-codimensional objects, one value
identifies one intersection. Finally the set of two co-dimensional
intersections among fractures is indicated by $\xi$, which can be viewed as a
union of $\sharp(\mathcal{P})$ points $\xi_{l}$ such that
$\xi=\cup_{l\in\mathcal{P}}\xi_{l}$. We suppose that if
$\xi_{l}\subset\gamma_{k}$, for some $l$ and $k$, then $\xi_{l}$ belongs to
the internal part of $\gamma_{k}$. Note that
$\gamma\cup\xi=\cap_{i\in\mathcal{N}}\Omega_{i}$. The aforementioned sets of
indexes are ordered in a natural way. See Fig. 3 as an example. Throughout the
paper generally, subscripts $i$ and $j$ denote quantities related with
fracture planes $\Omega_{i}$, subscript $k$ is associated with one co-
dimensional intersections $\gamma_{k}$, while subscript $l$ is used for two
co-dimensional intersections $\xi_{l}$. We indicate data and unknowns defined
on $\gamma$ with $\hat{\cdot}$.

$\Omega_{1}$$\Omega_{2}$$\Omega_{3}$$\gamma_{1}$$\gamma_{3}$$\xi$$\gamma_{1}$

Figure 3: Example of a DFN with three intersecting fractures, represented as
ellipses. The fractures form three one co-dimensional objects, the straight
lines, and one two co-dimensional object, the black dot.

Our objective is to compute the pressure field, indicated by $p$ and
$\hat{p}$, and the Darcy velocity field, indicated by $\bm{u}$ and
$\hat{\bm{u}}$, in $\Omega$ such that the Darcy equations are fulfilled. To
simplify the presentation we start by considering a single fracture
$\Omega_{i}$ without intersections. The normal of the fracture is indicated by
$\bm{n}_{i}$, and we define the normal projection matrix as
$N_{i}\vcentcolon=\bm{n}_{i}\otimes\bm{n}_{i}$ and the tangential projection
matrix as $T_{i}\vcentcolon=I-N_{i}$. In the following, for each fracture, we
make use of tangential divergence and gradient, defined as

$\displaystyle\nabla_{T_{i}}\cdot\vcentcolon=T_{i}:\nabla\quad\text{ and
}\quad\nabla_{T_{i}}\vcentcolon=T_{i}\nabla.$

In the sequel, we will sometimes drop the subscript on $T_{i}$ when there
should be no room for confusion.

###  2.2 Fracture flow

Following the work of [32, 20] we assume that the fracture permeability can be
written as a full elliptic second order tensor for the permeability in the
tangential space. We adopt a reduced order, or hybrid dimensional, model, to
describe $p_{i}$ and $\bm{u}_{i}$ in the fracture. The model explicitly takes
into account the hydraulic aperture $d_{i}$ of the fracture, which can change
in space along the fracture. We require that exist $d_{*}>0$ such that
$d_{i}\geq d_{*}$ for all $i$, i.e., the fracture does not degenerate. In the
forthcoming models we assume that the permeability and fracture aperture are
constant in each cell of the computational grid, this is a practical
assumption when the discretization is applied. The generalized Darcy equation
and the conservation of mass in $\Omega_{i}$ can be written as

(1) $\displaystyle\begin{aligned} &\bm{u}_{i}+\lambda_{i}\nabla
p_{i}=\bm{0}\\\ &\nabla\cdot\bm{u}_{i}=f_{i}\end{aligned}\quad\text{ in
}\Omega_{i}\quad\text{and}\quad p_{i}=g_{i}\text{ on }\partial\Omega_{i},$

where, for each $\Omega_{i}$, $\lambda_{i}\vcentcolon=d_{i}k_{i}$ is the
effective permeability, $k_{i}$ the tangential permeability, $f_{i}$ the
scalar sink or source, and $g_{i}$ the boundary condition. For simplicity of
the proofs, we assume Dirichlet conditions on all fractures, although it would
suffice to require only Dirichlet conditions on a non-zero measure portion of
the boundary. The case of where some fractures with only Neumann conditions
are imposed is studied in the numerical examples. The proofs can also be
extended to this case to the price of increased technicalities. Under these
assumptions, problem Eq. 1 is well posed, see [15].

###  2.3 Coupling Flow Between Fractures

Now, consider two fractures $\Omega_{i}$ and $\Omega_{j}$ such that a one co-
dimensional intersection occurs $\gamma_{k}=\Omega_{i}\cap\Omega_{j}$. We
indicate by $\bm{n}_{k}$ the set of unit vectors orthogonal to $\gamma_{k}$.
In the internal part of $\Omega_{i}$, i.e. in $\Omega_{i}\setminus\gamma_{k}$,
equation Eq. 1 is applied and for each internal point of $\gamma_{k}$, far
from its possible internal ending, we locally identify two sides of
$\Omega_{i}$ indicated $\Omega_{i}^{+}$ and $\Omega_{i}^{-}$. Since this
definition is local, for simplicity we keep uniform notation $\Omega_{i}^{+}$
and $\Omega_{i}^{-}$ for “the same side” of $\Omega_{i}$. We indicate data and
unknowns restricted to one side of a fracture with a $+$ or $-$ superscript.
For convenience $\partial\Omega_{i}$ represents the outer part of the
boundary, i.e., without the intersection. We indicate by $T_{i}\bm{n}_{k}$ the
unit normal vector of $\gamma_{k}$ lying in the plane of $\Omega_{i}$ and
pointing from $\Omega_{i}^{+}$ to $\Omega_{i}^{-}$. $T_{j}\bm{n}_{k}$ is
defined similarly. See Fig. 4 as an example.

$T_{j}\bm{n}_{k}$$T_{i}\bm{n}_{k}$$\Omega_{j}$$\Omega_{i}^{+}$$\Omega_{i}^{-}$$\gamma_{k}$

Figure 4: Example of two fractures with the notation introduced in the
derivation of the models.

Following [1, 2], at the intersection we assume the coupling conditions along
$\gamma_{k}$ as

(2) $\displaystyle\begin{aligned} &\sum_{m=i,j}\left\llbracket{\bm{u}_{m}\cdot
T_{m}\bm{n}_{k}}\right\rrbracket_{\gamma_{k}}=0\\\
&p_{i}^{+}|_{\gamma_{k}}=p_{i}^{-}|_{\gamma_{k}}=p_{j}^{+}|_{\gamma_{k}}=p_{j}^{-}|_{\gamma_{k}}\end{aligned}\quad\text{
on }\gamma_{k},$

where $\cdot|_{\gamma_{k}}$ is the trace operator from a fracture,
$\Omega_{i}$ or $\Omega_{j}$, to $\gamma_{k}$ and
$\left\llbracket{\cdot}\right\rrbracket_{\gamma_{k}}$ denotes the jump
operator across $\gamma_{k}$ defined as $\left\llbracket{\bm{u}_{m}\cdot
T_{m}\bm{n}_{k}}\right\rrbracket_{\gamma_{k}}\vcentcolon=\bm{u}_{m}^{+}\cdot
T_{m}\bm{n}_{k}|_{\gamma_{k}}-\bm{u}_{m}^{-}\cdot
T_{m}\bm{n}_{k}|_{\gamma_{k}}$, for $m=i,j$. Thus, while we allow for flow
from $\Omega_{i}$ to $\Omega_{j}$, the permeability of the intersection is so
high, or the width of the fracture so small, that the pressure drop is
negligible, see [32] for further considerations.

In subsurface flow, intersections between fractures, or faults, may be
significant conduits for fluid flow [36, 34]. It is therefore of interest to
extend the model Eq. 2 to allow for tangential flow along $\gamma_{k}$ and for
pressure jumps over the intersection, as studied in e.g. [20, 37, 12]. Flow in
the intersection is still modeled by a generalized Darcy’s law with an
additional source term from the surrounding fractures

(3a) $\displaystyle\begin{aligned}
&\hat{\bm{u}}_{k}+\hat{\lambda}_{k}\nabla\hat{{p}}_{k}=\bm{0}\\\
&\nabla\cdot\hat{\bm{u}}_{k}=\hat{f}_{k}+\sum_{m=i,j}\left\llbracket{\bm{u}_{m}\cdot
T_{m}\bm{n}_{k}}\right\rrbracket_{\gamma_{k}}\end{aligned}\text{in
}\gamma_{k}\quad\text{and}\quad\begin{aligned} &\hat{p}_{k}=\hat{g}_{k}&\text{
on }\partial\gamma_{k}\\\ &\hat{\bm{u}}_{k}\cdot\bm{n}_{\gamma_{k}}=0&\text{
on }\partial\gamma_{k}^{in}\end{aligned}$ where
$\hat{\lambda}_{k}\vcentcolon=d_{i}d_{j}\hat{k}(\gamma_{k})$ is the effective
permeability along $\gamma_{k}$, $\hat{k}(\gamma_{k})$ is the tangential
permeability of $\gamma_{k}$, and $\hat{g}_{k}$ is the boundary condition on
the, possibly empty, outer part (i.e., included in $\partial\Omega_{i}$ or
$\partial\Omega_{j}$) of the boundary of $\gamma_{k}$, for simplicity,
indicated by $\partial\gamma_{k}$. We define as $\gamma_{k}^{in}$ the,
possible empty, part of the boundary of $\gamma_{k}$ such that it lies
internally to $\Omega_{i}$ or $\Omega_{j}$. We identify by
$\bm{n}_{\gamma_{k}}$ the unit vector tangential to the intersection pointing
outward with respect to $\partial\gamma_{k}^{in}$. Note that $d_{i}$,
respectively $d_{j}$, is the aperture of fracture $\Omega_{i}$, respectively
$\Omega_{j}$. Their product represents the measure of the area normal to the
intersection. The differential operators are now defined along $\gamma_{k}$
and can be written in terms of local coordinates. The coupling condition Eq. 2
is generalized, for $m=i,j$, as (3b) $\displaystyle\begin{aligned}
&\bm{u}_{m}^{+}\cdot
T_{m}\bm{n}_{k}|_{\gamma_{k}}=-\widetilde{\lambda}_{m}\left(\hat{p}_{k}-p_{m}^{+}|_{\gamma_{k}}\right)\\\
&\bm{u}_{m}^{-}\cdot
T_{m}\bm{n}_{k}|_{\gamma_{k}}=-\widetilde{\lambda}_{m}\left(p_{m}^{-}|_{\gamma_{k}}-\hat{p}_{k}\right)\end{aligned}\quad\text{in
}\gamma_{k}$

where $\widetilde{\lambda}_{m}\vcentcolon=\widetilde{k}(\gamma_{k})/d_{m}$ is
the effective permeability of the fracture $\Omega_{m}$ in the direction
normal to the intersection $\gamma_{k}$, and $\widetilde{k}(\gamma_{k})$ is
the normal permeability of $\gamma_{k}$. Note that it is possible to have
different normal permeability for each fracture and even for each side of a
fracture, however to avoid confusion we assume a unique value for
$\widetilde{k}(\gamma_{k})$. Following [16, 20] the coupling conditions Eq. 3b
can be written as

$\displaystyle\begin{aligned}
&\left\llbracket{p_{m}}\right\rrbracket_{\gamma_{k}}=\widetilde{\lambda}_{m}^{-1}\left(\bm{u}_{m}^{+}\cdot
T_{m}\bm{n}_{k}|_{\gamma_{k}}+\bm{u}_{m}^{-}\cdot
T_{m}\bm{n}_{k}|_{\gamma_{k}}\right)\\\ &\left\llbracket{\bm{u}_{m}\cdot
T_{m}\bm{n}_{k}}\right\rrbracket_{\gamma_{k}}=\widetilde{\lambda}_{m}\left(p_{m}^{+}+p_{m}^{-}\right)\end{aligned}\quad\text{in
}\gamma_{k}$

where the jump of pressure and velocity across $\gamma_{k}$ are explicitly
written.

Finally, considering a single two co-dimensional intersection between
fractures, we have $\xi_{l}=\gamma_{k}\cap\gamma_{n}$. Following the idea
proposed in [12] we require the same model presented in Eq. 2 at the
intersection, we have

(4) $\displaystyle\begin{aligned}
&\sum_{m=k,n}\left\llbracket{\hat{\bm{u}}_{m}\cdot
T_{m}\bm{n}_{l}}\right\rrbracket_{\xi_{l}}=0\\\
&\hat{p}_{k}^{+}|_{\xi_{l}}=\hat{p}_{k}^{-}|_{\xi_{l}}=\hat{p}_{n}^{+}|_{\xi_{l}}=\hat{p}_{n}^{-}|_{\xi_{l}}\end{aligned}\quad\text{in
}\xi_{l}$

where $\cdot|_{\xi_{l}}$ is the trace operator from a one co-dimensional
intersection to $\xi_{l}$ and $T_{m}\bm{n}_{l}$ is a unit vector, pointing
outward from $\xi_{l}$, and lying on $\gamma_{m}$. While it is possible to
extend the model to allow for jumps in pressure and velocity over $\xi_{l}$,
see [20, 22], this is not considered herein.

We can now define two problems defined on the network $\Omega$ which will be
used in the sequel. The first problem considers continuous coupling conditions
at the fracture intersections, it is the following.

######  Problem 1 (DFN Darcy flow - continuous coupling).

Considering the aforementioned data and assumptions on $\Omega$, the model
problem with continuous coupling conditions is to find $(\bm{u},p)$ in
$\Omega$ such that Eq. 1 is valid for all $i\in\mathcal{N}$ and Eq. 2 is valid
for all $k\in\mathcal{I}$.

The second problem allows a flow tangential to the fracture intersection and a
possible pressure jump across the one co-dimensional intersections.

######  Problem 2 (DFN Darcy flow - discontinuous coupling).

Considering the aforementioned data and assumptions on $\Omega$, the model
problem with discontinuous coupling conditions is to find $(\bm{u},p)$ in
$\Omega$ such that Eq. 1 is valid for all $i\in\mathcal{N}$, Eq. 3 is valid
for all $k\in\mathcal{I}$, and Eq. 4 is valid for all $l\in\mathcal{P}$.

######  Remark 1.

In the case that $\gamma_{k}$ is the one co-dimensional intersection of more
than two fractures, coupling conditions Eq. 2 can be generalized extending the
summation of the normal fluxes on all the fractures and assuming the pressure
continuity at the intersection. A similar consideration can be done for the
coupling conditions Eq. 3 where the summation is extended to all the fractures
involved and the definition of $\hat{\lambda}_{k}$ should include the measure
of the cross section of $\gamma_{k}$, now simply represented as a rectangle.

######  Remark 2.

In the case that $\xi_{l}$ is the two co-dimensional intersection of more than
two fractures, coupling conditions Eq. 4 can be generalized assuming the
pressure continuity at the intersection for all the one co-dimensional
objects.

##  3 Well posedness

In this part the weak formulation of the problem 1 and 2 and the related
functional spaces are introduced. We provide also well posedness results for
both problems.

###  3.1 DFN Darcy flow - continuous coupling

We define $\|\cdot\|_{E}:L^{2}(E)\rightarrow\mathbb{R}$ as the usual
$L^{2}$-norm and $\left(\cdot,\cdot\right)_{E}:L^{2}(E)\times
L^{2}(E)\rightarrow\mathbb{R}$ which is the usual $L^{2}$-scalar product. To
simplify the presentation we use this notation for both scalar and vector
functions. We indicate by $[\cdot]_{i}$ and $[\cdot]_{ij}$ the $i$-th,
respectively $ij$-th, component of the vector, respectively of the matrix, in
the square brackets. We assume also pressure boundary conditions assigned to
each $\partial\Omega_{i}$ by a function $g_{i}\in
H^{\frac{1}{2}}(\partial\Omega_{i})$, a source term $f_{i}\in
L^{2}(\Omega_{i})$, and $[\lambda_{i}]_{jk}\in L^{\infty}(\Omega_{i})$. Given
a regular domain $E\in\Omega_{i}$, we define the functional spaces
$Q(E)\vcentcolon=L^{2}(E)$ for the scalar fields and
$V(E)\vcentcolon=\left\\{\bm{v}\in[L^{2}(E)]^{3}:\nabla\cdot\bm{v}\in
L^{2}(E)\right\\}$ for the vector fields, which are Hilbert spaces. It is
worth to remember that we are dealing with differential operators on the
tangent planes of the fractures. With an abuse of notation we write
$V_{i}=V(\Omega_{i})$ and $Q_{i}=Q(\Omega_{i})$ for each fracture. We endow
$Q_{i}$ with the usual $L^{2}$-norm and $V_{i}$ with the $H_{\rm div}$-norm on
the tangent space, defined as
$\|\bm{v}\|_{V_{i}}^{2}\vcentcolon=\|\bm{v}\|_{\Omega_{i}}^{2}+\|\nabla\cdot\bm{v}\|_{\Omega_{i}}^{2}$.
The global spaces for the network $\Omega$ are defined as
$Q(\Omega)\vcentcolon=\prod_{i\in\mathcal{N}}Q_{i}$ with norm and
$V(\Omega)\vcentcolon=\prod_{i\in\mathcal{N}}V_{i}$, with an abuse of notation
we indicate $V=V(\Omega)$ and $Q=Q(\Omega)$. Their norms are defined as
$\|q\|_{Q}^{2}\vcentcolon=\sum_{i\in\mathcal{N}}\|q_{i}\|_{Q_{i}}^{2}$ and
$\|\bm{v}\|_{V}^{2}\vcentcolon=\sum_{i\in\mathcal{N}}\|\bm{v}_{i}\|_{V_{i}}^{2}$.
Following the standard procedure we can derive the weak formulation of 1 which
requires the definition of the following bilinear forms

$\displaystyle a\left(\cdot,\cdot\right):V\times V\rightarrow\mathbb{R}:\quad
a\left(\bm{u},\bm{v}\right)\vcentcolon=\sum_{i\in\mathcal{N}}a_{i}\left(\bm{u}_{i},\bm{v}_{i}\right),\,\,a_{i}\left(\bm{u}_{i},\bm{v}_{i}\right)\vcentcolon=\left(\lambda_{i}^{-1}\bm{u}_{i},\bm{v}_{i}\right)_{\Omega_{i}}$
$\displaystyle b\left(\cdot,\cdot\right):V\times Q\rightarrow\mathbb{R}:\quad
b\left(\bm{u},q\right)\vcentcolon=\sum_{i\in\mathcal{N}}b_{i}\left(\bm{u}_{i},q_{i}\right),\,\,b_{i}\left(\bm{u}_{i},q_{i}\right)\vcentcolon=-\left(\nabla\cdot\bm{u}_{i},q_{i}\right)_{\Omega_{i}}$

and the following functionals which include the boundary data of the problem
and the scalar source term

$\displaystyle B\left(\cdot\right):V\rightarrow\mathbb{R}:\quad
B\left(\bm{v}\right)\vcentcolon=-\sum_{i\in\mathcal{N}}\langle\bm{v}_{i}\cdot\bm{n}_{\partial\Omega_{i}},g_{i}\rangle_{\partial\Omega_{i}}$
$\displaystyle F\left(\cdot\right):Q\rightarrow\mathbb{R}:\quad
F\left(q\right)\vcentcolon=-\sum_{i\in\mathcal{N}}\left(q_{i},f_{i}\right)_{\Omega_{i}}$

where the dual pairing is defined as
$\langle\cdot,\cdot\rangle_{\partial\Omega_{i}}:H^{-\frac{1}{2}}(\partial\Omega_{i})\times
H^{\frac{1}{2}}(\partial\Omega_{i})\rightarrow\mathbb{R}$ and
$\bm{n}_{\partial\Omega_{i}}$ is the unit normal pointing outward from
$\partial\Omega_{i}$ and tangent to $\Omega_{i}$. With this definitions we
introduce the following problem.

######  Problem 3 (Weak formulation of 1).

The weak formulation of 1 is to find $(\bm{u},p)\in V\times Q$ such that

(5) $\displaystyle\begin{aligned}
&a\left(\bm{u},\bm{v}\right)+b\left(\bm{v},p\right)=B\left(\bm{v}\right)&\forall\bm{v}\in
V\\\ &b\left(\bm{u},q\right)=F\left(q\right)&\forall q\in Q\end{aligned}$

######  Theorem 3.1.

3 is well posed.

######  Proof 3.2.

Using standard arguments, mainly Cauchy-Schwarz inequality, it is possible to
prove the boundedness, and then the continuity, of the bilinear forms $a$ and
$b$ on their spaces as well as the functionals $F$ and $B$. To show coercivity
of $a$, on the kernel of $b$ we consider a function $\bm{w}\in V$ such that
$b\left(\bm{w},q\right)=0$ for all $q\in Q$. We then have
$\nabla_{T_{i}}\cdot\bm{w}_{i}=0$ almost everywhere in $\Omega_{i}$ and
$\|\bm{w}\|_{V}^{2}=\sum_{i\in\mathcal{N}}\|\bm{w}_{i}\|_{\Omega_{i}}^{2}$,
thus

$\displaystyle
a\left(\bm{w},\bm{w}\right)=\sum_{i\in\mathcal{N}}\left(\lambda_{i}^{-1}\bm{w}_{i},\bm{w}_{i}\right)_{\Omega_{i}}\gtrsim\|\bm{w}\|_{V}^{2}.$

To simplify the proof we now consider only two fractures, being the case of
multiple fractures a straightforward extension where the following techniques
are properly applied to each fracture intersection. To verify the inf-sup
condition on $b$ we can have two possibilities: either the intersection does
not reach the boundary of the fracture or its contrary. In both cases, let $q$
be a function in $Q$ and consider the following auxiliary problem

$\displaystyle\begin{aligned}
&-\nabla_{T_{i}}\cdot\nabla_{T_{i}}\varphi_{i}=q_{i}&\text{in }\Omega_{i}\\\
&\varphi_{i}=0&\text{on
}\partial\Omega_{i}\cup\left\\{\gamma\cap\Omega_{i}\right\\}\end{aligned},$

which are decoupled problems, one for each fracture. If the intersection
completely cuts $\Omega_{i}$ then, following [17] the previous problem admits
a unique solution on each disconnected part of $\Omega_{i}$, called
$\Omega_{i}^{*}$. We have $\varphi_{i}^{*}\in H^{2}(\Omega_{i}^{*})$ $\forall
i$, such that
$\|\varphi_{i}^{*}\|_{H^{2}(\Omega_{i}^{*})}\lesssim\|q_{i}\|_{\Omega_{i}^{*}}$.
With an abuse of notation, indicating the $H^{2}$-broken norm with the same
symbol, we obtain
$\|\varphi_{i}\|_{H^{2}(\Omega_{i})}\lesssim\|q_{i}\|_{\Omega_{i}}$. We
consider now $\bm{v}\in V$ such that
$\bm{v}|_{\Omega_{i}^{*}}=\nabla_{T_{i}}\varphi_{i}^{*}$, we have
$-\nabla_{T_{i}}\cdot\bm{v}|_{\Omega_{i}^{*}}=q_{i}$ in $\Omega_{i}^{*}$ and

$\displaystyle\|\bm{v}\|_{V}^{2}=\sum_{i\in\mathcal{N}}\|\nabla_{T_{i}}\varphi_{i}\|^{2}_{\Omega_{i}}+\|q_{i}\|^{2}_{\Omega_{i}}\lesssim\|q\|_{Q}^{2}.$

With this choice of $\bm{v}$ we obtain the boundedness from below of the
bilinear form $b$

$\displaystyle
b\left(\bm{v},q\right)=\sum_{i\in\mathcal{N}}-\left(\nabla_{T_{i}}\cdot\bm{v},q\right)_{\Omega_{i}}=\|q\|_{Q}^{2}\gtrsim\|q\|_{Q}\|\bm{v}\|_{V}.$

Following [15] we conclude that in the case where the intersection boundaries
coincide with the boundary of the fracture planes, 3 is well posed. We
consider now the problem of partially immersed intersection, i.e. , one ending
part touches a fracture boundary and the other is immersed in the fracture.
The fully immersed intersection is a straightforward extension, where the same
technique is applied to both the ending parts. If the intersection is immersed
in the fracture planes, following [17, 23] the previous auxiliary problem
admits a unique solution $\varphi_{i}\in
H^{\frac{3}{2}-\epsilon}(\Omega_{i})$, $\epsilon>0$ and for all
$i\in\mathcal{N}$. Define a regularized kernel
$\theta\in{C}^{\infty}(T\Omega_{i})$ with compact support in $T\Omega_{i}$,
that is, the tangent space of $\Omega_{i}$, centred in the immersed part of
$\gamma$. We introduce also the corresponding family of mollifiers
$\theta_{\zeta}(x)=\zeta^{-n}\theta(x/\zeta)$ for $\zeta>0$ and
$n\in\mathbb{N}$. We consider now $\bm{v}_{\zeta}\in V\cap
C^{\infty}(T\Omega_{i})$ such that
$[\bm{v}_{\zeta}]_{j}=\phi_{\zeta}\ast[\nabla_{T_{i}}\varphi_{i}]_{j}$ for
$j=1,2$. We have

$\displaystyle\|\bm{v}_{\zeta}\|_{V}^{2}=\|\theta_{\zeta}\ast\nabla_{T_{i}}\varphi_{i}\|^{2}_{\Omega_{i}}+\|\nabla_{T_{i}}\cdot\theta_{\zeta}\ast\nabla_{T_{i}}\varphi_{i}\|^{2}_{\Omega_{i}}$

where, with an abuse of notation, the mollifier is applied to each component
of $\bm{v}_{i}$. Using the properties of the mollifiers, we can bound the
first term as

$\displaystyle\|\theta_{\zeta}\ast\nabla_{T_{i}}\varphi_{i}\|_{\Omega_{i}}\leq\|\nabla_{T_{i}}\varphi_{i}\|_{\Omega_{i}}\leq\|\varphi_{i}\|_{H^{1}(\Omega_{i})}\leq\|q\|_{\Omega_{i}},$

while the second term can be estimated as

$\displaystyle\|\nabla_{T_{i}}\cdot\theta_{\zeta}\ast\nabla_{T_{i}}\varphi_{i}\|_{\Omega_{i}}=\|\theta_{\zeta}\ast\nabla_{T_{i}}\cdot\nabla_{T_{i}}\varphi_{i}\|_{\Omega_{i}}\leq\|\nabla_{T_{i}}\cdot\nabla_{T_{i}}\varphi_{i}\|_{\Omega_{i}}=\|q\|_{\Omega_{i}}.$

We obtain $\|\bm{v}_{\zeta}\|_{V}\lesssim\|q_{i}\|_{Q}$. With this choice we
have for a $\Omega_{i}$

$\displaystyle
b_{i}\left(\bm{v}_{\zeta},q_{i}\right)=-\left(\nabla_{T_{i}}\cdot\bm{v}_{\zeta},q_{i}\right)_{\Omega_{i}}=-\left(\nabla_{T_{i}}\cdot\theta_{\zeta}\ast\nabla_{T_{i}}\varphi_{i},q_{i}\right)_{\Omega_{i}}=$
$\displaystyle=-\left(\theta_{\zeta}\ast\nabla_{T_{i}}\cdot\nabla_{T_{i}}\varphi_{i},q_{i}\right)_{\Omega_{i}}=\left(\theta_{\zeta}\ast
q_{i},q_{i}\right)_{\Omega_{i}}.$

Using again the property of the mollifiers we obtain that for
$n\rightarrow\infty$: $\bm{v}_{\zeta}\xrightarrow{L^{2}}\bm{v}$, with
$\bm{v}\in V$, and $\theta_{\zeta}\ast q_{i}\xrightharpoonup{L^{2}}q_{i}$.
Following [15] we conclude also in this case that 3 is well posed.

###  3.2 DFN Darcy flow - discontinuous coupling

We consider now the functional setting to present the weak formulation of
problem 2. Referring to Eq. 3a and Eq. 3b, we assume pressure boundary
conditions assigned to each ending point $\partial\gamma_{k}$ by a scalar
$\hat{g}_{i}\in\mathbb{R}$, a source term defined in $\hat{f}_{k}\in
L^{2}(\gamma_{k})$, an effective tangential permeability in
$\hat{\lambda}_{k}\in L^{\infty}(\gamma_{k})$, and an effective normal
permeability with regularity $\widetilde{\lambda}_{m}\in
L^{\infty}(\gamma_{k})$. Motivated by e.g. [32], we introduce a new family of
spaces, one for each fracture $\Omega_{i}$ by

$\displaystyle W_{i}\vcentcolon=\left\\{\bm{v}\in V_{i}:\bm{v}^{+}\cdot
T_{i}\bm{n}_{k}|_{\gamma_{k}}\in L^{2}(\gamma_{k})\text{ and }\bm{v}^{-}\cdot
T_{i}\bm{n}_{k}|_{\gamma_{k}}\in L^{2}(\gamma_{k}),\forall
k\in\mathcal{G}_{i}\right\\}$

and their composition for the vector fields
$W\vcentcolon=\prod_{i\in\mathcal{N}}W_{i}$, which are Hilbert spaces endowed
with norms

$\displaystyle\|\bm{v}\|_{W_{i}}^{2}\vcentcolon=\|\bm{v}\|_{V_{i}}^{2}+\sum_{k\in\mathcal{G}_{i}}\|\bm{v}^{+}\cdot
T_{i}\bm{n}_{k}\|_{\gamma_{k}}^{2}+\|\bm{v}^{-}\cdot
T_{i}\bm{n}_{k}\|_{\gamma_{k}}^{2}\quad\text{ and
}\quad\|\bm{v}\|_{W}^{2}\vcentcolon=\sum_{i\in\mathcal{N}}\|\bm{v}_{i}\|_{W_{i}}^{2}.$

It is worth to notice that we require more regularity on the intersections for
$W_{i}$ than $V_{i}$, where implicitly we assume $H^{-1/2}$-regularity, in
order to properly take into account the coupling conditions Eq. 3b. This
assumption is related to Robin-type boundary conditions for problem in mixed
form, see [32] for a more detailed investigation. For the intersections we
introduce the functional spaces for both scalar fields as
$\hat{Q}_{k}\vcentcolon=L^{2}(\gamma_{k})$, with the usual $L^{2}$-norm, and
vector fields as

$\displaystyle\hat{V}_{k}\vcentcolon=\left\\{\hat{\bm{v}}\in[L^{2}(\gamma_{k})]^{3}:\nabla\cdot\hat{\bm{v}}\in
L^{2}(\gamma_{k})\right\\}\quad\text{with}\quad\|\hat{\bm{v}}\|_{\hat{V}_{k}}^{2}\vcentcolon=\|\hat{\bm{v}}\|_{\gamma_{k}}^{2}+\|\nabla\cdot\hat{\bm{v}}\|_{\gamma_{k}}^{2}.$

The global spaces are $\hat{V}\vcentcolon=\prod_{k\in\mathcal{I}}\hat{V}_{k}$,
with norm
$\|\hat{\bm{v}}\|_{\hat{V}}^{2}\vcentcolon=\sum_{k\in\mathcal{I}}\|\hat{\bm{v}}_{k}\|_{\hat{V}_{k}}^{2}$
for the vector fields, and
$\hat{Q}\vcentcolon=\prod_{k\in\mathcal{I}}\hat{Q}_{k}$, with norm
$\|\hat{q}\|_{\hat{Q}}^{2}\vcentcolon=\sum_{k\in\mathcal{I}}\|\hat{q}_{k}\|^{2}_{\hat{Q}_{k}}$
for the scalar fields. Finally the spaces for the coupled problem are
$U\vcentcolon=W\times\hat{V}$, with induced norm from $W$ and $\hat{V}$, and
$O\vcentcolon=Q\times\hat{Q}$, with induced norm from $Q$ and $\hat{Q}$. All
the aforementioned spaces are Hilbert spaces. We introduce the bilinear forms
associated with the intersections as

$\displaystyle\hat{a}\left(\cdot,\cdot\right):\hat{V}\times\hat{V}\rightarrow\mathbb{R}:\quad\hat{a}\left(\hat{\bm{u}},\hat{\bm{v}}\right)\vcentcolon=\sum_{k\in\mathcal{I}}\hat{a}_{k}\left(\hat{\bm{u}}_{k},\hat{\bm{v}}_{k}\right),\,\,\hat{a}_{k}\left(\hat{\bm{u}}_{k},\hat{\bm{v}}_{k}\right)\vcentcolon=\left(\hat{\lambda}_{k}^{-1}\hat{\bm{u}}_{k},\hat{\bm{v}}_{k}\right)_{\gamma_{k}}$
$\displaystyle\hat{b}\left(\cdot,\cdot\right):\hat{V}\times\hat{Q}\rightarrow\mathbb{R}:\quad\hat{b}\left(\hat{\bm{u}},\hat{q}\right)\vcentcolon=\sum_{k\in\mathcal{I}}\hat{b}_{k}\left(\hat{\bm{u}}_{k},\hat{q}_{k}\right),\,\,\hat{b}_{k}\left(\hat{\bm{u}}_{k},\hat{q}_{k}\right)\vcentcolon=-\left(\nabla\cdot\hat{\bm{u}}_{k},\hat{q}_{k}\right)_{\gamma_{k}}.$

The global bilinear forms for the coupled problem are
$\alpha\left(\cdot,\cdot\right):U\times U\rightarrow\mathbb{R}$ and
$\beta\left(\cdot,\cdot\right):U\times O\rightarrow\mathbb{R}$

$\displaystyle\alpha\left(\left(\bm{u},\hat{\bm{u}}\right),\left(\bm{v},\hat{\bm{v}}\right)\right)\vcentcolon=a\left(\bm{u},\bm{v}\right)+\hat{a}\left(\hat{\bm{u}},\hat{\bm{v}}\right)+cc^{+}_{1}(\bm{u},\bm{v})+cc^{-}_{1}(\bm{u},\bm{v})$
$\displaystyle\beta\left(\left(\bm{u},\hat{\bm{u}}\right),\left(q,\hat{q}\right)\right)\vcentcolon=b\left(\bm{u},q\right)+\hat{b}\left(\hat{\bm{u}},\hat{q}\right)+cc_{2}(\bm{u},\hat{q})$

where the bilinear forms associated with the coupling conditions are

(6) $\displaystyle cc^{+}_{1}\left(\cdot,\cdot\right):V\times
V\rightarrow\mathbb{R}:\quad
cc^{+}_{1}(\bm{u},\bm{v})\vcentcolon=\sum_{i\in\mathcal{N}}\sum_{j\in\mathcal{G}_{i}}\left(\widetilde{\lambda}_{i}^{-1}\bm{u}_{i}^{+}\cdot
T_{i}\bm{n}_{j},\bm{v}_{i}^{+}\cdot T_{i}\bm{n}_{j}\right)_{\gamma_{j}}$ (7)
$\displaystyle
cc_{2}\left(\cdot,\cdot\right):V\times\hat{Q}\rightarrow\mathbb{R}:\quad
cc_{2}(\bm{u},\hat{q})\vcentcolon=\sum_{i\in\mathcal{N}}\sum_{j\in\mathcal{G}_{i}}\left(\left\llbracket{\bm{u}_{i}\cdot
T_{i}\bm{n}_{j}}\right\rrbracket_{\gamma_{j}},\hat{q}_{j}\right)_{\gamma_{j}},$

$cc_{1}^{-}(\cdot,\cdot)$ follow immediatly. The functional $B$ is extended
naturally on $U$ and we define

$\displaystyle\Phi:O\rightarrow\mathbb{R}:\quad\Phi\left(\left(q,\hat{q}\right)\right)\vcentcolon=-\sum_{i\in\mathcal{N}}\left(q_{i},f_{i}\right)_{\Omega_{i}}-\sum_{k\in\mathcal{I}}\left(\hat{q}_{k},\hat{f}_{k}\right)_{\gamma_{k}}$

######  Problem 3.3 (Weak formulation of 2).

The weak formulation of 2 is to find
$(\left(\bm{u},\hat{\bm{u}}\right),\left(p,\hat{p}\right))\in U\times O$ such
that

(8) $\displaystyle\begin{aligned}
&\alpha\left(\left(\bm{u},\hat{\bm{u}}\right),\left(\bm{v},\hat{\bm{v}}\right)\right)+\beta\left(\left(\bm{v},\hat{\bm{v}}\right),\left(p,\hat{p}\right)\right)=B\left(\left(\bm{v},\hat{\bm{v}}\right)\right)&\forall\left(\bm{v},\hat{\bm{v}}\right)\in
U\\\
&\beta\left(\left(\bm{u},\hat{\bm{u}}\right),\left(q,\hat{q}\right)\right)=\Phi\left(\left(q,\hat{q}\right)\right)&\forall\left(q,\hat{q}\right)\in
O\end{aligned}$

######  Theorem 3.4.

3.3 is well posed.

######  Proof 3.5.

Using standard arguments, mainly Cauchy-Schwarz inequality, it is possible to
prove the boundedness, and then the continuity, of the bilinear forms $\alpha$
and $\beta$ on their spaces as well as the functionals $\Phi$ and $B$.
Considering a function $\left(\bm{w},\hat{\bm{w}}\right)\in U$ such that
$\beta\left(\left(\bm{w},\hat{\bm{w}}\right),\left(q,\hat{q}\right)\right)=0$
for all $\left(q,\hat{q}\right)\in O$, we have
$\nabla_{T_{i}}\cdot\bm{w}_{i}=0$ almost everywhere in $\Omega_{i}$ and
$\nabla\cdot\hat{\bm{w}}_{k}=0$ almost everywhere in $\gamma_{k}$. The norm of
$\left(\bm{w},\hat{\bm{w}}\right)$ becomes

$\displaystyle\|\left(\bm{w},\hat{\bm{w}}\right)\|_{U}^{2}=\sum_{i\in\mathcal{N}}\|\bm{w}_{i}\|_{\Omega_{i}}^{2}+\sum_{k\in\mathcal{G}_{i}}\|\bm{w}^{+}_{i}\cdot
T_{i}\bm{n}_{k}\|_{\gamma_{k}}^{2}+\|\bm{w}^{-}_{i}\cdot
T_{i}\bm{n}_{k}\|_{\gamma_{k}}^{2}+\sum_{k\in\mathcal{I}}\|\hat{\bm{w}}_{k}\|_{\gamma_{k}}^{2}.$

It is possible to show that
$\alpha\left(\left(\bm{w},\hat{\bm{w}}\right),\left(\bm{w},\hat{\bm{w}}\right)\right)\gtrsim\|\left(\bm{w},\hat{\bm{w}}\right)\|_{U}^{2}$,
thus $\alpha$ is coercive on the kernel of $\beta$.

Given a function $(q,\hat{q})\in O$, taking inspiration from [32, 20], we
introduce the following auxiliary problems for fractures

$\displaystyle\begin{aligned}
&-\nabla_{T_{i}}\cdot\nabla_{T_{i}}\varphi_{i}=q_{i}&\text{in }\Omega_{i}\\\
&\nabla_{T_{i}}\varphi_{i}^{+}\cdot T_{i}\bm{n}_{j}=\hat{q}_{j}&\text{on
}\gamma_{j},j\in\mathcal{G}_{i}\\\ &\nabla_{T_{i}}\varphi_{i}^{-}\cdot
T_{i}\bm{n}_{j}=-\hat{q}_{j}&\text{on }\gamma_{j},j\in\mathcal{G}_{i}\\\
&\left\llbracket{\varphi_{i}}\right\rrbracket_{\gamma_{j}}=0&\text{on
}\gamma_{j},j\in\mathcal{G}_{i}\\\ &\varphi_{i}=0&\text{on
}\partial\Omega_{i}\end{aligned}$

and for the intersections

$\displaystyle\begin{aligned}
&-\nabla\cdot\nabla\hat{\varphi}_{k}=\hat{q}_{k}&\text{in }\gamma_{k}\\\
&\hat{\varphi}_{k}=0&\text{on }\partial\gamma_{k}\end{aligned}.$

We consider only the case where the intersection reaches the fracture
boundary; intersection endpoints inside the fracture plane can be handled
similar to the proof of Theorem 3.1. Following [17] the previous problems
admits a unique solution on each disconnected part of $\Omega_{i}$, indicated
by $\Omega_{i}^{*}$. In the first case we have $\varphi_{i}^{*}\in
H^{2}(\Omega_{i}^{*})$ $\forall i$, such that
$\|\varphi_{i}^{*}\|_{H^{2}(\Omega_{i}^{*})}\lesssim\|q_{i}\|_{\Omega_{i}^{*}}+\sum_{j\in\mathcal{G}_{i}}\|\hat{q}_{j}\|_{\gamma_{j}}$.
We consider now $\bm{v}\in W$ such that
$\bm{v}|_{\Omega_{i}^{*}}=\nabla_{T_{i}}\varphi_{i}^{*}$, we have
$-\nabla_{T_{i}}\cdot\bm{v}|_{\Omega_{i}^{*}}=q_{i}$ in $\Omega_{i}^{*}$,
$\bm{v}|_{\Omega_{i}^{*}}\cdot T_{i}\bm{n}_{j}=\hat{q}_{j}$ for all
$j\in\mathcal{G}_{i}$, and

$\displaystyle\|\bm{v}_{i}\|_{W_{i}}^{2}=\|\nabla_{T_{i}}\varphi_{i}\|^{2}_{\Omega_{i}}+\|q_{i}\|^{2}_{\Omega_{i}}+2\sum_{j\in\mathcal{G}_{i}}\|\hat{q}_{j}\|_{\gamma_{j}}^{2}\lesssim\|q_{i}\|_{Q_{i}}^{2}+\sum_{j\in\mathcal{G}_{i}}\|\hat{q}_{j}\|_{Q_{j}}^{2}\leq\|q\|_{Q}^{2}+\|\hat{q}\|_{\hat{Q}}^{2}.$

For the second family of problems we obtain the existence of
$\hat{\varphi}_{k}\in H^{2}(\gamma_{k})$ $\forall k$, such that
$\|\hat{\varphi}_{k}\|_{H^{2}(\gamma_{k})}\lesssim\|\hat{q}_{k}\|_{\gamma_{k}}$.
Considering $\hat{\bm{v}}\in\hat{V}$ such that
$\hat{\bm{v}}_{k}=\nabla\hat{\varphi}_{k}$, we have
$-\nabla\cdot\hat{\bm{v}}_{k}=\hat{q}_{k}$ and

$\displaystyle\|\hat{\bm{v}}\|_{\hat{V}}^{2}=\sum_{k\in\mathcal{I}}\|\nabla_{T_{i}}\hat{\varphi}_{k}\|^{2}_{\gamma_{k}}+\|\hat{q}_{k}\|^{2}_{\gamma_{k}}\lesssim\|\hat{q}\|_{\hat{Q}}^{2}.$

With this choice of $\bm{v}$ and $\hat{\bm{v}}$ we obtain the boundedness from
below of the bilinear form $\beta$

$\displaystyle\beta\left(\left(\bm{v},\hat{\bm{v}}\right),\left(q,\hat{q}\right)\right)=\sum_{i\in\mathcal{N}}\|q_{i}\|_{Q_{i}}^{2}+2\sum_{j\in\mathcal{G}_{i}}\|\hat{q}_{j}\|_{\gamma_{j}}^{2}+\sum_{k\in\mathcal{I}}\|\hat{q}_{k}\|_{\gamma_{k}}^{2}\gtrsim\|(q,\hat{q})\|_{O}\|(\bm{v},\hat{\bm{v}})\|_{U}.$

Thus the inf-sup condition is fulfilled, and following [15] we conclude that
3.3 is well posed.

##  4 Discrete approximation

In this part we present the numerical discretization of 3 and 3.3. We extend
the virtual element method for mixed problem presented in [14, 5, 7] to the
DFN setting. For simplicity we start the discretization of 3 considering only
the network composed by a single fracture $\Omega=\Omega_{i}$, the extension
of more than one fracture is trivial but requires a heavier notation, while
the approximation of Eq. 2 will be discussed below. For realistic fracture
networks, the solutions will commonly have low regularity due to
heterogeneities, and we therefore limit ourselves to lowest order methods.

###  4.1 Discrete setting

Let us consider $\mathcal{T}(\Omega)$ a tessellation of a single fracture
$\Omega$ into non-overlapping polygons $E$ and the union of the edges $e$ as
$\mathcal{E}(\Omega)=\left\\{e\in\partial E\right\\}$. We denote by $h_{E}$
the diameter of $E$, by $\bm{x}_{E}$ the centre of $E$, by $\bm{n}_{e}$ with
$e\in\mathcal{E}(E)$ the unit normal of $e$ pointing outward with respect to
the internal part of $E$, and by $h$ the maximum diameter in
$\mathcal{T}(\Omega)$. When there is no room for confusion we write
$\mathcal{T}$ instead of $\mathcal{T}(\Omega)$ and $\mathcal{E}$ instead of
$\mathcal{E}(\Omega)$. We assume that there exist $\rho_{E}\in\mathbb{R}^{+}$
such that $E$ is star-shaped with respect to every point of a disc of radius
$\rho_{E}h_{E}$ and $\left|e\right|\geq\rho_{E}h_{E}$, for all
$e\in\mathcal{E}(E)$. An example of elements is presented in Fig. 5.

$E_{1}$$E_{2}$$E_{3}$

Figure 5: Example of three elements admitted in the discretization. The blue
points represent the centre $\bm{x}_{E}$ for each element $E$, the red points
define the edges $e$, and the green arrows are the normal vector to the edges.
$E_{1}$ is a rectangular element but composed of 14 edges and $E_{2}$ is a
non-convex element.

We introduce the following local spaces for the element $E$ of $\mathcal{T}$:
for the scalar fields $Q_{h}(E)\vcentcolon=\left\\{q\in
Q(E):q\in\mathbb{P}_{0}(E)\right\\}$ and for the vector fields

$\displaystyle V_{h}(E)\vcentcolon=\left\\{\bm{v}\in
V(E):\bm{v}\cdot\bm{n}_{e_{i}}|_{e_{i}}\in\mathbb{P}_{0}(e_{i})\forall
e_{i}\in\mathcal{E}(E),\nabla\cdot\bm{v}\in\mathbb{P}_{0}(E),\nabla\times\bm{v}=\bm{0}\right\\}.$

It is worth to notice that we do not reconstruct the velocity inside an
element, since in most of the applications the normal flux from a face is the
physical quantity of interest. The global spaces for the fracture $\Omega$ are
defined by

$\displaystyle Q_{h}(\Omega)\vcentcolon=\left\\{q\in Q(\Omega):q|_{E}\in
Q_{h}(E)\,\forall E\in\mathcal{T}\right\\}$ $\displaystyle
V_{h}(\Omega)\vcentcolon=\left\\{\bm{v}\in V(\Omega):\bm{v}|_{E}\in
V_{h}(E)\,\forall E\in\mathcal{T}\right\\}.$

The degrees of freedom, simply d.o.f.’s , for $Q_{h}(\Omega)$ are constant
value in each element of the mesh, while for $\bm{v}\in V_{h}(\Omega)$ we
consider

$\displaystyle\left(\bm{v}\cdot\bm{n}_{e},q\right)_{e}\quad\forall
e\in\mathcal{E}\text{ and }\forall q\in\mathbb{P}_{0}(e).$

We indicate by $N_{\textit{dof}}=\sharp\mathcal{E}(E)$ the number of d.o.f.’s
for a generic element $E$. To simplify the presentation, we construct the
approximation locally to each element $E\in\mathcal{T}$ and then use the fact
that the velocity d.o.f.’s are single valued for all $e\in\mathcal{E}$, the
global approximation can be build up. When it is clear from the context, we
consider the bilinear form introduced in the previous section restricted to an
element of the mesh. For the element $E\in\mathcal{T}$ we define the local
canonical base for both spaces: for $Q_{h}(E)$ is trivial, while for
$V_{h}(E)$ we indicate with
$\operatorname*{base}(V_{h}(E))=\left\\{\bm{\varphi}_{\omega}\right\\}_{\omega=1}^{N_{\textit{dof}}}$
such that $\bm{\varphi}_{\omega}\in V_{h}(E)$ and given $\bm{v}\in V_{h}(E)$
we have
$\bm{v}=\sum_{\omega=1}^{N_{\textit{dof}}}[\bm{v}]_{\omega}\bm{\varphi}_{\omega}$,
with $[\bm{v}]_{\omega}\in\mathbb{R}$. As a consequence of these facts we have
$\left(\bm{\varphi}_{\omega}\cdot\bm{n}_{e_{j}},1\right)_{e_{j}}=\delta_{\omega
j}$ which means
$\left|e_{\omega}\right|\bm{\varphi}_{\omega}\cdot\bm{n}_{e_{\omega}}=1$, as
well as $\left(\nabla\cdot\bm{\varphi}_{\omega},1\right)_{E}=1$ which means
$\left|E\right|\nabla\cdot\bm{\varphi}_{\omega}=1$. With these choices, it is
possible to compute directly the bilinear form $b(\cdot,\cdot)$ as well as the
functionals $F(\cdot)$ and $B(\cdot)$. However the computation of
$a(\cdot,\cdot)$ is more tricky and requires the introduction of a projection
operator $\Pi_{0}$, which gives an approximation of $a(\cdot,\cdot)\approx
a_{h}(\cdot,\cdot)$. We will see that the explicit evaluation of the basis
functions $\bm{\varphi}_{i}$ is not needed in the internal part of $E$, but
the request that the order of convergence is preserved when using
$a_{h}(\cdot,\cdot)$ instead of $a(\cdot,\cdot)$. We introduce the following
local space

$\displaystyle{\mathcal{V}}(E)\vcentcolon=\left\\{\bm{v}\in
V_{h}(E):\bm{v}=\lambda\nabla v,\text{ for }v\in\mathbb{P}_{1}(E)\right\\},$

where $\nabla$ stands for the tangential gradient on the current fracture. We
define the projection operator $\Pi_{0}:V(E)\rightarrow\mathcal{V}(E)$ such
that $a(\bm{v}-\Pi_{0}\bm{v},\bm{w})=0$ for all $\bm{w}\in\mathcal{V}(E)$. To
introduce the discrete bilinear forms, we start by having for
$\bm{u},\bm{v}\in V_{h}(E)$

(9) $\displaystyle
a(\bm{u},\bm{v})=a(\Pi_{0}\bm{u},\Pi_{0}\bm{v})+a((I-\Pi_{0})\bm{u},(I-\Pi_{0})\bm{v}),$

since the cross terms are zero due to $\Pi_{0}$.

###  4.2 Computation of the local $H_{\rm div}$-mass matrix

We focus now on the approximation of the first term in (9). Since
$\Pi_{0}\bm{u}\in V_{h}(E)$, it can be expanded in the canonical base
$\Pi_{0}\bm{u}=\sum_{\omega=1}^{N_{\textit{dof}}}[\Pi_{0}\bm{u}]_{\omega}\bm{\varphi}_{\omega}$
with $\bm{\varphi}_{\omega}=\lambda\nabla\varphi_{\omega}$,
$\varphi_{\omega}\in\mathbb{P}_{1}(E)$. The trial function
$\Pi_{0}\bm{v}\in\mathcal{V}$ then $\exists v\in\mathbb{P}_{1}(E)$ such that
$\Pi_{0}\bm{v}=\lambda\nabla v$. Moreover we can write $\varphi_{\omega}$, as
well as $v$, using a monomial expansion as

$\displaystyle\varphi_{\omega}=\sum_{m=1}^{2}s_{\omega}^{m}m_{m}\quad\text{with}\quad
m_{m}(\bm{x})\vcentcolon=\dfrac{[\bm{x}]_{m}-[\bm{x}_{E}]_{m}}{h_{E}}\quad\text{and}\quad
s_{\omega}^{m}\in\mathbb{R}.$

It is worth to notice that the coordinates in the definition of the monomials
are defined in the tangential space of $\Omega$. The reason why
$m_{0}(\bm{x})=1$ is not included will become clear below. We obtain for
$\omega$

$\displaystyle\left(\lambda^{-1}\bm{\varphi}_{\omega},\bm{v}\right)_{E}=-\left(\nabla\cdot\bm{\varphi}_{\omega},v\right)_{E}+\sum_{e\in\mathcal{E}(E)}\left(\bm{\varphi}\cdot\bm{n}_{e},v\right)_{e}=-\dfrac{1}{\left|E\right|}\left(1,v\right)_{E}+\dfrac{1}{\left|e_{\omega}\right|}\left(1,v\right)_{e_{\omega}},$

which is computable using the monomial expansion of $v$. Indeed, by using the
monomial expansion for both $\varphi_{\omega}$ and $v$, we obtain

$\displaystyle\left(\lambda^{-1}\bm{\varphi}_{\omega},\lambda\nabla
m_{j}\right)_{E}=\left(\nabla\varphi_{\omega},\lambda\nabla
m_{j}\right)_{E}=\sum_{i=1}^{2}s_{\omega}^{i}\left(\nabla m_{i},\lambda\nabla
m_{j}\right)_{E}\quad\text{for}\quad j=1,2$

the scalar product in the previous expression is computable. Defining the
matrices and vectors $G\in\mathbb{R}^{2\times 2}$,
$\bm{f}_{\omega},\bm{s}_{\omega}\in\mathbb{R}^{2}$, $F\in\mathbb{R}^{2\times
N_{\textit{dof}}}$, and $\Pi^{*}\in\mathbb{R}^{2\times N_{\textit{dof}}}$ as

$\displaystyle[G]_{ij}\vcentcolon=\left(\lambda\nabla m_{i},\nabla
m_{j}\right)_{E},\quad[\bm{f}_{\omega}]_{i}\vcentcolon=-\dfrac{1}{\left|E\right|}\left(1,m_{i}\right)_{E}+\dfrac{1}{\left|e_{\omega}\right|}\left(1,m_{i}\right)_{e_{\omega}},$
$\displaystyle\bm{s}_{\omega}\vcentcolon=[s^{1}_{\omega},s^{2}_{\omega}]^{\top}=G^{-1}\bm{f}_{\omega},\quad
F\vcentcolon=[\bm{f}_{1}|\ldots|\bm{f}_{N_{\textit{dof}}}],\quad\Pi^{*}\vcentcolon=G^{-1}F,$

we have $s_{\omega}^{k}=[\Pi^{*}]_{\omega k}$. The final expression thus reads

$\displaystyle
a(\Pi_{0}\bm{\varphi}_{\omega},\Pi_{0}\bm{\varphi}_{\theta})=\sum_{i=1}^{N_{\textit{dof}}}\sum_{j=1}^{N_{\textit{dof}}}s_{\omega}^{i}s_{\theta}^{j}\left(\lambda\nabla
m_{i},\nabla
m_{j}\right)_{E}=\sum_{i=1}^{N_{\textit{dof}}}\sum_{j=1}^{N_{\textit{dof}}}[\Pi^{*}]_{\omega
i}[G]_{ij}[\Pi^{*}]_{\theta j}=$
$\displaystyle=[(\Pi^{*})^{\top}G\Pi^{*}]_{\omega\theta}.$

For the second, stability, term in (9), we follow [14, 5, 7] and approximate
the term as

$\displaystyle
a((I-\Pi_{0})\bm{\varphi}_{\omega},(I-\Pi_{0})\bm{\varphi}_{\theta})\approx
s\left((I-\Pi_{0})\bm{\varphi}_{\omega},(I-\Pi_{0})\bm{\varphi}_{\theta}\right)\vcentcolon=$
$\displaystyle\vcentcolon=\varsigma\sum_{i=1}^{N_{\textit{dof}}}\textit{dof}_{i}\left((I-\Pi_{0})\bm{\varphi}_{\omega}\right)\textit{dof}_{i}\left((I-\Pi_{0})\bm{\varphi}_{\theta}\right),$

where $s(\cdot,\cdot):V_{h}(E)\times V_{h}(E)\rightarrow\mathbb{R}$ is the
bilinear form associated to the stabilization and $\varsigma\in\mathbb{R}$ is
a suitable parameter that will be explained later on. We give an expression of
$\textit{dof}_{i}\left((I-\Pi_{0})\bm{\varphi}_{\omega}\right)$. Since
$\mathcal{V}\subset V_{h}(E)$ we expand each projected element of the
canonical base on the canonical base itself, i.e.
$\Pi_{0}\bm{\varphi}_{\omega}=\sum_{i=1}^{N_{\textit{dof}}}\pi_{\omega}^{i}\bm{\varphi}_{i}$
for $\omega=1,\ldots,N_{\textit{dof}}$ with
$\pi_{\omega}^{j}=\textit{dof}_{j}(\Pi_{0}\bm{\varphi}_{\omega})$, which is
also

$\displaystyle\Pi_{0}\bm{\varphi}_{\omega}=\sum_{j=1}^{2}s_{\omega}^{j}\lambda\nabla
m_{j}=\sum_{j=1}^{2}s_{\omega}^{j}\sum_{i=1}^{N_{\textit{dof}}}\textit{dof}_{i}(\lambda\nabla
m_{j})\bm{\varphi}_{i}=\sum_{i=1}^{N_{\textit{dof}}}\left(\sum_{j=1}^{2}s_{\omega}^{j}\textit{dof}_{i}(\lambda\nabla
m_{j})\right)\bm{\varphi}_{i},$

obtaining $\pi_{\omega}^{i}$ equal to the term in the brackets. This is
computable introducing the matrix $D\in\mathbb{R}^{2\times 2}$, with
$[D]_{ij}\vcentcolon=\textit{dof}_{i}(\lambda\nabla m_{j})$, as
$\pi_{\omega}^{i}=[D\Pi^{*}]_{i\omega}$. For the stabilization we finally
obtain
$s\left((I-\Pi_{0})\bm{\varphi}_{\omega},(I-\Pi_{0})\bm{\varphi}_{\theta}\right)=\varsigma[(I-D\Pi^{*})^{\top}(I-D\Pi^{*})]_{\omega\theta}$.
Introducing the bilinear form $a_{h}(\cdot,\cdot):V_{h}(E)\times
V_{h}(E)\rightarrow\mathbb{R}$ as

$\displaystyle
a_{h}(\bm{\varphi}_{\omega},\bm{\varphi}_{\theta})\vcentcolon=a\left(\Pi_{0}\bm{\varphi}_{\omega},\Pi_{0}\bm{\varphi}_{\theta}\right)+s\left((I-\Pi_{0})\bm{\varphi}_{\omega},(I-\Pi_{0})\bm{\varphi}_{\theta}\right)=$
$\displaystyle=\left[(\Pi^{*})^{\top}G\Pi^{*}+\varsigma(I-D\Pi^{*})^{\top}(I-D\Pi^{*})\right]_{\omega\theta},$

the local approximation of 3 for a single fracture is

######  Problem 4.1 (Local discrete formulation of 3).

the discrete approximation of the weak problem in $E$ is find
$\left(\bm{u},p\right)\in V_{h}(E)\times Q_{h}(E)$ such that

$\displaystyle\begin{aligned}
&a_{h}\left(\bm{u},\bm{\varphi}_{\theta}\right)+b\left(p,\bm{\varphi}_{\theta}\right)=-\left(\bm{\varphi}_{\theta}\cdot\bm{n}_{\partial\Omega},g\right)_{\partial\Omega\cap
e_{\theta}}\quad\forall\bm{\varphi}_{\theta}\in\operatorname*{base}(V_{h}(E))\\\
&b\left(1,\bm{u}\right)=-(f_{i},1)_{E}\end{aligned}.$

###  4.3 Fracture intersection

We consider two intersecting fractures $\Omega_{i}$ and $\Omega_{j}$ such that
$\gamma_{k}=\Omega_{i}\cap\Omega_{j}$. The general case is extendible with an
analogous procedure. We discretize $\gamma_{k}$ as the union of consecutive
edges of $\mathcal{E}(\Omega_{i})\cap\mathcal{E}(\Omega_{j})$, indicated by
$\mathcal{T}(\gamma_{k})$. For each fracture we double the velocity d.o.f.’s
in $\mathcal{T}(\gamma_{k})$. We enforce Eq. 2 by using Lagrange multipliers
along $\mathcal{T}(\gamma_{k})$, one for each edge involved.

If two intersecting fractures are present, then the coupled model 2 can be
considered, and we need to introduce a proper discretization also for the
intersection as well as for the coupling condition Eq. 3b. The procedure for
the former is similar to the derivation of the discrete system for the
fractures but in a mono-dimensional framework, i.e. in $\mathcal{T}(\gamma)$.
For simplicity we consider a single two-codimensional object indicated by
$\gamma$, the extension to multiple intersections is trivial but requires a
heavier notation. The approximation of Eq. 4 will be discussed below. We
consider the discrete spaces
$\hat{Q}_{h}(E)\vcentcolon=\left\\{\hat{q}\in\hat{Q}(E):\,\hat{q}\in\mathbb{P}_{0}(E)\right\\}$
with $E\in\mathcal{T}(\gamma)$ for the pressure, for the velocity we have

$\displaystyle\hat{V}_{h}(E)\vcentcolon=\left\\{\hat{\bm{v}}\in\hat{V}(E):\hat{\bm{v}}\cdot\bm{n}_{e_{i}}|_{e_{i}}\in\mathbb{R}\,\forall
e_{i}\in\mathcal{E}(E),\nabla\cdot\hat{\bm{v}}\in\mathbb{P}_{0}(E),\nabla\times\hat{\bm{v}}=\bm{0}\right\\},$

where $\bm{n}_{e_{i}}$ is the unit tangential vector of $\gamma$ which point
outward from $E$. The global spaces $\hat{Q}_{h}(\gamma)$ and
$\hat{V}_{h}(\gamma)$ follow naturally. The d.o.f.’s for $\hat{Q}_{h}(\Omega)$
are piece-wise constant in each element of $\mathcal{T}(\gamma)$, while for
$\hat{\bm{v}}\in\hat{V}_{h}(\gamma)$, given an ending point $e$ of an element
$E$, we consider $\hat{\bm{v}}\cdot\bm{n}_{e}|_{e}$. Also in this case, we
introduce a projection operator
$\hat{\Pi}_{0}:\hat{V}(E)\rightarrow\hat{\mathcal{V}}(E)$ such that
$\hat{a}(\hat{\bm{v}}-\hat{\Pi}_{0}\hat{\bm{v}},\hat{\bm{w}})=0$ for all
$\hat{\bm{w}}\in\hat{\mathcal{V}}(E)$ and where

$\displaystyle\hat{\mathcal{V}}(E)\vcentcolon=\left\\{\hat{\bm{v}}\in\hat{V}_{h}(E):\hat{\bm{v}}=\hat{\lambda}\nabla\hat{v},\text{
for }\hat{v}\in\mathbb{P}_{1}(E)\right\\}.$

Following the same process described for an element belonging to a fracture
mesh, we obtain the local matrix formulation for the $H_{\rm div}$-mass matrix

$\displaystyle\hat{a}(\hat{\bm{\varphi}}_{\omega},\hat{\bm{\varphi}}_{\theta})\approx\hat{a}_{h}(\hat{\bm{\varphi}}_{\omega},\hat{\bm{\varphi}}_{\theta})=\left[(\hat{\Pi}^{*})^{\top}\hat{G}\hat{\Pi}^{*}+\hat{\varsigma}(I-\hat{D}\hat{\Pi}^{*})^{\top}(I-\hat{D}\hat{\Pi}^{*})\right]_{\omega\theta}$

where $\hat{\bm{\varphi}}_{\omega},\hat{\bm{\varphi}}_{\theta}$ are elements
of the base for $\hat{V}_{h}(E)$, $\hat{\varsigma}\in\mathbb{R}$ is a proper
stabilization parameter, and $\hat{\Pi}^{*},\hat{G}$, and $\hat{D}$ are
suitable matrices. The approximation of the $\hat{b}(\cdot,\cdot)$ bilinear
form follows from the definition of the d.o.f.’s for the velocity and the
pressure. The discrete formulation of the coupling condition presented in Eq.
6 can be derived using the d.o.f.’s introduced previously. In particular, for
each intersection $\gamma_{j}$ of a fracture $\Omega_{i}$, we double the
d.o.f. associated with the velocity obtaining the discrete form of
$\bm{u}^{+}\cdot T_{i}\bm{n}_{j}|_{\gamma_{j}}$ and $\bm{u}^{-}\cdot
T_{i}\bm{n}_{j}|_{\gamma_{j}}$. With this choice the implementation of Eq. 6
is immediate.

In the case where two intersections meet in a point $\xi$, coupling conditions
Eq. 4 should be adopted. For each intersection we double the reduced velocity
d.o.f.’s in $\xi$ and we enforce the condition by using Lagrange multipliers
in $\xi$.

###  4.4 Stabilization term

To conclude, we discuss now the stabilization parameters $\varsigma$ and
$\hat{\varsigma}$ introduced previously. Following [14, 5, 7], to obtain a
proper error decay we require that exist
$\iota_{*},\iota^{*},\hat{\iota}_{*},\hat{\iota}^{*}\in\mathbb{R}^{+}$,
independent from the discretization size, such that

(12) $\displaystyle\begin{gathered}\iota_{*}a(\Pi_{0}\bm{v},\Pi_{0}\bm{v})\leq
s((I-\Pi_{0})\bm{v},(I-\Pi_{0})\bm{v})\leq\iota^{*}a(\Pi_{0}\bm{v},\Pi_{0}\bm{v})\qquad\forall\bm{v}\in
V(\Omega)\\\
\hat{\iota}_{*}\hat{a}(\Pi_{0}\hat{\bm{v}},\Pi_{0}\hat{\bm{v}})\leq\hat{s}((I-\hat{\Pi}_{0})\hat{\bm{v}},(I-\hat{\Pi}_{0})\hat{\bm{v}})\leq\hat{\iota}^{*}\hat{a}(\hat{\Pi}_{0}\hat{\bm{v}},\hat{\Pi}_{0}\hat{\bm{v}})\qquad\forall\hat{\bm{v}}\in\hat{V}(\Omega)\end{gathered}$

The stability term introduced previously for the fractures fulfils
automatically the request, however for highly heterogeneous fractures in the
permeability a proper scaling is recommended, e.g.
$\varsigma_{i}=\|\lambda^{-1}_{i}\|_{L^{\infty}(\Omega_{i})}$ for each
fracture $\Omega_{i}$. To fulfil the second request in Eq. 12 we can compute
explicitly the local matrices involved. Given a segment $E$ of length $h_{E}$
and supposing that the effective permeability $\hat{\lambda}$ is constant in
$E$, we obtain

$\displaystyle(\hat{\Pi}^{*})^{\top}\hat{G}\hat{\Pi}^{*}=\dfrac{h_{E}}{4\hat{\lambda}}\begin{bmatrix}1&-1\\\
-1&1\end{bmatrix}\quad\text{and}\quad(I-\hat{D}\hat{\Pi}^{*})^{\top}(I-\hat{D}\hat{\Pi}^{*})=\dfrac{1}{2}\begin{bmatrix}1&1\\\
1&1\end{bmatrix}.$

We see that choosing $\hat{\varsigma}=h_{E}/\hat{\lambda}$ is enough to scale
properly the stability term for the intersection.

######  Remark 4.2.

It is important to note that if $\bm{u}=\Pi_{0}\bm{u}$ or
$\bm{u}=\Pi_{0}\bm{v}$ we have the identity:
$a_{h}(\bm{u},\bm{v})=a(\bm{u},\bm{v})$. The same is valid for the
intersection flow.

##  5 The computational grid

Realistic fracture networks can have a highly complex geometry, and
correspondingly the construction of the computational grid is challenging. In
particular, intersections between fractures add complexity, since for the
types of discretization considered herein, intersections are treated as
constraints in the gridding algorithm. This often leads to a high number of
cells, and depending on intersection geometry, also low quality elements. The
requirement that the grid conforms to intersections can be avoided by using
specialized numerical methods, however their implementation tends to be
tedious.

Since the virtual element method can handle almost any polygon type, we make
two modifications of what can be considered a broadly used DFN meshing
algorithm, both aimed at alleviating the computational cost. Other approaches
can be found in [33, 26, 28], to name a few. Consider the DFN composed by two
intersecting fractures, $\Omega_{i}$ and $\Omega_{j}$, presented in Fig. 6,
the extension to several fractures is straightforward.

$\Omega_{i}$$\Omega_{j}$$\gamma_{k}$

Figure 6: The geometry of the two fractures and their intersection marked in
light green, which are used in this section to present the meshing and
coarsening strategy.

Our first modification is to build each mesh fracture separately and then link
them together through the intersection $\gamma_{k}$. This reduces the meshing
problem to a set of decoupled 2d domain with internal constraints from the
intersection lines, and allows us to apply established 2d meshing software
without adaptation for the coupling between domains. We have applied the
library Triangle, see [38], which is a fast and robust triangular grid
generator that allows for internal constraints. See Fig. 7 on the left as an
example. The independent meshing also means the fractures can be meshed in
parallel.

Figure 7: Left: a triangulated fracture where the intersection is preserved.
In general the conformity of the triangles across the intersection is not
preserved. Right: Grouping of triangular cells after coarsening. Clusters of
triangles with the same colour will form a cell in the coarse mesh.

When a mesh of a fracture is created, our second modification is to coarsen
the grid by cell agglomeration so as to lighten the total computational cost.
Our approach is motivated by the coarsening in algebraic multigrid methods,
e.g. [39], and the coarsening is carried out independently for each fracture.
As a measure of connectedness between, we consider a two-point flux
approximation (TPFA) discretization of Eq. 1, and denote the discretization
matrix $A$. To preserve fracture intersections in the coarse grid, we treat
intersections as a boundary for the TPFA discretization, and we explicitly
prohibit a coarse cell to cover both sides of the end of an intersection. The
coarsening ratio is determined by the parameter $c_{depth}$. Details on the
coarsening algorithm can be found in Appendix A, or in the literature on
algebraic multigrid methods, see e.g. [39].

Fig. 7 on the right shows the clustering map after the application of the
algorithm. We clearly see that the resulting elements may be concave or even
non star-shaped (but a finite union of star-shaped), e.g. the big orange
element on the right part of the mesh. It is worth to notice that the process
of creating and coarsening the meshes is embarrassingly parallel. Once the
mesh of each fracture is created, we simply compute a co-refinement of the
edges lying on $\gamma_{k}$ from both fractures. Splitting the edges of the
fractures at the intersection and considering a coherent numeration of the
elements and edges the global mesh is build up. Refer to Fig. 8 as an example
of final mesh with and without the coarsening strategy. The effect of the
parameter $c_{depth}$ is illustrated in Fig. 9. We also note that in
anisotropic media the coarse grid will to some extent adapt to the
preferential flow directions, see Fig. 10.

Figure 8: The meshes of the two fractures linked together before (left) and
after (right) coarsening. The points at the intersection are highlighted in
white and preserved in the coarsening. The original grid has 93 elements and
153 edges, while the coarse grid has 23 elements and 82 edges. The elements at
the intersection are non-matching but conforming.

An example of the grid overlapped to the pressure field is represented in Fig.
1 in the introduction.

##  6 Examples

In this section we validate the models presented in the previous sections
through several tests and examples. In particular in Section 6.1 we highlight
the potentiality of the coarsening algorithm presented in Section 5 applied in
our context. Section 6.2 contains several tests with a single as well as
multiple fractures to put in evidence the error decay and others properties of
the numerical solution. In Section 6.3 we present a numerical example showing
the importance to use model presented in 2. Finally in Section 6.4 we present
the solution on a realistic geometry. Our implementation is carried out within
the framework of the Matlab Reservoir Simulation Toolbox (MRST) [31].

We make use of the following

$\displaystyle\text{sparsity}=\dfrac{\text{nnz}(A)}{\text{size}(A)^{2}},$

with $\text{nnz}(A)$ the number of non-zero values of the matrix $A$ and
$\text{size}(A)$ the number of rows (or columns) of $A$, as a measure for the
sparsity of a matrix.

###  6.1 Coarsening

In this subsection we investigate the ability of the coarsening algorithm,
introduced in Section 5 to generate coarse meshes where the direction of the
cells reflects the underlying physic. To that end, we consider two problems
with different effective permeability: isotropic and anisotropic, the latter
in a heterogeneous setting. In both cases we start from a single fracture
discretized by a triangular grid, then using the algorithm we coarsen the mesh
with different level of refinement: $c_{depth}=1$, $c_{depth}=3$, and
$c_{depth}=5$. The fracture is the unit square rotate by $\pi/4$ along the $z$
axis.

In the first case the effective permeability is the identity tensor. Fig. 9
shows the meshes obtained with different level of coarsening.

Figure 9: Example of meshes used for the example in Section 6.1 with isotropic
effective permeability tensor. From the left: triangular mesh, the coarse
meshes with $c_{depth}=1$, $c_{depth}=3$, and $c_{depth}=5$, respectively.

The resulting cells of the meshes do not have a preferred alignment. The
relative sparsity of the VEM matrix for each coarse level is: $1.4\cdot
10^{-3}$ for the triangular grid, $3.9\cdot 10^{-3}$ for $c_{depth}=1$,
$2.3\cdot 10^{-2}$ for $c_{depth}=3$, and $2\cdot 10^{-1}$ for $c_{depth}=5$,
which is a natural effect of this coarsening algorithm. Finally the number of
edges (minimum, average, and maximum) and cells for each level is: $(3,5,6)$
and 628 for $c_{depth}=1$, $(8,12,15)$ and 122 for $c_{depth}=3$, and
$(29,35,41)$ and 11 for $c_{depth}=5$. The resulting meshes respect the
theoretical requests for the VEM.

The second case deals with an anisotropic and heterogeneous effective
permeability tensor. We divide the fracture into four sectors: top-left part
($x\leq 0$ and $y>0$), top-right part ($x>0$ and $y>0$), bottom-left part
($x\leq 0$ and $y\leq 0$), and bottom-right part ($x>0$ and $y\leq 0$). In the
first and last sector we set $\lambda={\rm diag}[1,100]$, while in the others
we have $\lambda={\rm diag}[100,1]$. Fig. 10 shows the meshes obtained with
different level of coarsening.

Figure 10: Example of meshes used for the example in Section 6.1 for the
anisotropic permeability tensor. From the left: triangular mesh, the coarse
meshes with $c_{depth}=1$, $c_{depth}=3$, and $c_{depth}=5$, respectively.

The resulting cells of the meshes have a preferred alignment related to the
underlying permeability which respect both the anisotropy and the
heterogeneity. The sparsity of the VEM matrix for each coarse level is:
$1.4\cdot 10^{-3}$ for the triangular grid, $4.5\cdot 10^{-3}$ for
$c_{depth}=1$, $1.2\cdot 10^{-2}$ for $c_{depth}=3$, and $3.1\cdot 10^{-2}$
for $c_{depth}=5$, which is a natural effect of this coarsening algorithm.
Finally the number of edges (minimum, average, and maximum) and cells for each
level is: $(3,5,8)$ and 557 for $c_{depth}=1$, $(4,8,12)$ and 236 for
$c_{depth}=3$, and $(5,14,24)$ and 95 for $c_{depth}=5$. Compared with the
previous case we notice that the number of cells and the variability in the
number of edges is greater, which is needed to preserve the anisotropy of the
final mesh. The smaller sparsity of the matrices is a direct consequence of
this fact. Also in this case the resulting meshes respect the theoretical
requests for the VEM.

We conclude that the coarsening presented in Section 5 is a valid tool to
create a coarse mesh which is adapted with the underlying physic. In the next
examples we consider this coarsening strategy in the computation of error
decay.

###  6.2 Convergence evidence

In this part we present tests that show numerical evidence of convergence of
the approximate to the exact solution as well as the order of convergence. In
the first part, Section 6.2.1, a single fracture is considered with different
types of meshes. In the second example two intersecting fractures, with
coupling conditions described in 1 are applied. Finally two fractures are
considered where the flow is described also in their intersection. In the
following experiments we make use of the relative errors for both the pressure
and the projected velocity, we use the following expressions:

$\displaystyle
err(p)=\dfrac{\|p-p_{ex}\|_{\Omega}}{\|p_{ex}\|_{\Omega}}\quad\text{and}\quad
err(\Pi_{0}\bm{u})=\dfrac{\|\Pi_{0}\bm{u}-\bm{u}_{ex}\|_{\Omega}}{\|\bm{u}_{ex}\|_{\Omega}}$

where $p_{ex}$ and $\bm{u}_{ex}$ are the exact solutions computed in the
centre of each cell. The same type of errors are computed also for the
intersection when the model described in 2 is considered. In addition,
furthers details are included in Appendix B.

####  6.2.1 Single fracture

We consider a single fracture $\Omega$ constructed applying a rotation matrix
to the points of the unit square $[0,1]^{2}$. The rotation is $\pi/4$ along
the $x$ axis. We assume unit effective permeability $\lambda$ and source term
equal to

$\displaystyle
f(x,y,z)=7\,z-4\,\sin\\!\left(\pi\,y\right)+2\,{\pi}^{2}\,y^{2}\,\sin\\!\left(\pi\,y\right)-8\,\pi\,y\,\cos\\!\left(\pi\,y\right).$

We set pressure boundary conditions such that the exact solution of the
problem is the following

$\displaystyle
p_{ex}(x,y,z)=x^{2}\,z+4\,y^{2}\,\sin\\!\left(\pi\,y\right)-3\,z^{3}.$

The exact Darcy velocity, in Cartesian coordinates, is

$\displaystyle\bm{u}_{ex}(x,y,z)=\left[\begin{array}[]{c}-2\,x\,z\\\
0.5(9\,z^{2}-x^{2})-4\,y\,\sin\\!\left(\pi\,y\right)-2\,\pi\,y^{2}\,\cos\\!\left(\pi\,y\right)\\\
0.5(9\,z^{2}-x^{2})-4\,y\,\sin\\!\left(\pi\,y\right)-2\,\pi\,y^{2}\,\cos\\!\left(\pi\,y\right)\end{array}\right].$

It is possible to verify that $0\leq p_{ex}\leq 1.337$ in $\Omega$. We
consider four different families of grids to analyse the order of convergence
for both $p$ and $\Pi_{0}\bm{u}$, first build in the unit square and then
mapped in $\Omega$. We have: Cartesian grids, “coarse grids” where the
algorithm from Section 5 is applied with $c_{depth}=2$ on Cartesian grids,
triangular grids, and “random grids” where, starting from Cartesian grids, the
internal nodes are randomly moved. In the latter a sanity check of the mesh is
performed to avoid degenerate cases. See Fig. 11 as an example. It is worth to
notice that in the case of “coarse grids” and “random grids” the cells may not
be convex.

Figure 11: Example of meshes used for the example in Section 6.2.1. From the
left: coarse mesh, triangular mesh, and random mesh. The pattern of the coarse
mesh remains similar for subsequent refinements. The meshes are presented in
the $(x,y)$ plane. The Cartesian mesh is omitted.

In Fig. 12 we report the order of convergence for both $p$ and $\Pi_{0}\bm{u}$
for the four mesh families and we can observe that, in all the cases, the
order of convergence respect the decay derived from the theory in [7]. Table 1
in Appendix B reports the detailed values for this example, as well as the
sparsity of the matrix, and the minimum and maximum values of the solution. It
is interesting to notice that in some cases the discrete maximum or minimum
principle is violated.

Figure 12: In the figure the error decay for both $p$, in the left plot, and
$\Pi_{0}\bm{u}$, in the right plot, for the example in Section 6.2.1. We add
the references $\mathcal{O}(h)$ and $\mathcal{O}(h^{2})$ to aid the
comparison.

####  6.2.2 Two fractures with intersection

In the second example we analyse the error decay for two intersecting
fractures, where 1 is considered. For this test we adapt the problem in
Subsubsection 5.3.1 of [10] to our setting. Each fracture is build from an
ellipses discretized by 8 segments. We consider two orthogonal fractures,
depicted in Fig. 13, defined as

$\displaystyle\Omega_{1}=\left\\{(x,y,z)\in\mathbb{R}^{3}:z^{2}+4y^{2}\leq
1,x=0\right\\}$
$\displaystyle\Omega_{2}=\left\\{(x,y,z)\in\mathbb{R}^{3}:x^{2}+4y^{2}\leq
1,z=0\right\\}.$

Figure 13: Example of meshes and numerical solutions for the example in
Section 6.2.2. From the left: triangular mesh, coarse mesh with $c_{depth}=2$,
coarse mesh with $c_{depth}=4$, and coarse mesh with $c_{depth}=5$. The
solution has range in $[0,3.6]$ with a “Blue to Red Rainbow” colour map.

We assume pressure boundary condition on both $\partial\Omega_{1}$ and
$\partial\Omega_{2}$, unit effective permeability $\lambda$ and source as

(15)
$\displaystyle\begin{gathered}f(x,y,z)=\begin{dcases*}8y(1-y)-8(z-1)^{2}&for
$z\geq 0$\\\ 8y(1-y)-8(z+1)^{2}&for $z<0$\end{dcases*}\text{in }\Omega_{1}\\\
f(x,y,z)=\begin{dcases*}8y(1-y)-8(x+\zeta)^{2}&for $x\geq 0$\\\
8y(1-y)-8(x-\zeta)^{2}&for $x<0$\end{dcases*}\text{in
}\Omega_{2}\end{gathered}$

With this data the exact solutions $p_{ex}$ and $\bm{u}_{ex}$ are

(16c)
$\displaystyle\begin{gathered}p_{ex}(x,y,z)=\begin{dcases*}4y(1-y)(z-1)^{2}&for
$z\geq 0$\\\ 4y(1-y)(z+1)^{2}&for $z<0$\end{dcases*}\text{in }\Omega_{1}\\\
p_{ex}(x,y,z)=\begin{dcases*}4y(1-y)(x+\zeta)^{2}&for $x\geq 0$\\\
4y(1-y)(x-\zeta)^{2}&for $x<0$\end{dcases*}\text{in }\Omega_{2}\end{gathered}$
(16f)
$\displaystyle\begin{gathered}\bm{u}_{ex}(x,y,z)=\begin{dcases*}\left[0;\,\,\,8y(z-1)^{2};\,\,\,-8y(1-y)(z-1)\right]^{\top}&for
$z\geq 0$\\\ \left[0;\,\,\,8y(z+1)^{2};\,\,\,-8y(1-y)(z+1)\right]^{\top}&for
$z<0$\end{dcases*}\text{in }\Omega_{1}\\\
\bm{u}_{ex}(x,y,z)=\begin{dcases*}\left[-8y(1-y)(x+\zeta);\,\,\,8y(x+\zeta)^{2};\,\,\,0\right]^{\top}&for
$x\geq 0$\\\
\left[-8y(1-y)(x-\zeta);\,\,\,8y(x-\zeta)^{2};\,\,\,0\right]^{\top}&for
$x<0$\end{dcases*}\text{in }\Omega_{2},\end{gathered}$

we assume $\zeta=1$. We consider four different families of discretization to
analyse the order of convergence for both $p$ and $\Pi_{0}\bm{u}$ in both
fractures. We have: triangular grids, “coarse grids” with $c_{depth}=2$,
$c_{depth}=4$, and $c_{depth}=5$. In all the coarse cases we apply the
algorithm from Section 5 on triangular grids. See Fig. 13 as an example. It is
worth to notice that in the case of “coarse grids” the cells may not be convex
and even they may be not star-shaped.

Figure 14: In the figure the error decay for both $p$, in the left plot, and
$\Pi_{0}\bm{u}$ for the first example in Section 6.2. We add the references
$\mathcal{O}(h)$ and $\mathcal{O}(h^{2})$ to aid the comparison.

In Fig. 14 we report the $L^{2}$-relative errors for both the projected
velocity $\Pi_{0}\bm{u}$ and the pressure $p$. We observe that the error decay
for the latter is quadratic and the former is linear with respect to the mesh
size. Moreover, for the same $h$, the error is higher for the triangular grid
and smaller when $c_{depth}$ increases. One possible explanation is that in
the first case each pressure d.o.f. is linked with 3 velocity d.o.f.’s , while
for $c_{depth}=2$ the average number is 6, for $c_{depth}=4$ we have 20 edges,
and for $c_{depth}=5$ the average number of edges is 37. This fact can help
the accuracy for the meshes with more velocity d.o.f.’s . Note that the grid
size $h$ reported is the average mesh size, which in some cases exhibits an
erratic reduction compared to the steady decrease of the grid size in the
triangular grid. Table 2 in Appendix B reports the detailed values for this
example, as well as the sparsity of the matrix, the minimum and maximum values
of the solution, number of edges for each cell.

######  Remark 6.1.

In this test the meshes at the one co-dimensional intersection are matching
just because the procedure to construct the two meshes is the same, except the
mapping from the reference domain.

####  6.2.3 Two fractures with intersection flow

In this example we analyse the error decay for the coupled system 2, where two
fractures intersect. The fractures $\Omega_{1}$ and $\Omega_{2}$ are the same
as the test shown in 6.2.2, while the intersection $\gamma$ is defined as
$\gamma=\left\\{(x,y,z)\in\mathbb{R}^{3}:x=0,\,z=0,\,y\in[0,1]\right\\}$. The
analytical solution of the problem $p_{ex}$ and $\bm{u}_{ex}$ is Eq. 16 with
$\zeta=-1$ for the fractures and $\hat{p}_{ex}$ and $\hat{\bm{u}}_{ex}$ for
the intersection

$\displaystyle\hat{p}_{ex}(x,y,z)=5y(1-y)\quad\text{in
}\gamma\quad\text{and}\quad\hat{\bm{u}}_{ex}(x,y,z)=\left[0;\,\,\,5+10y;\,\,\,0\right]^{\top}\quad\text{in
}\gamma.$

We consider pressure boundary conditions on both fractures as well as on the
intersections. The effective permeability for the fractures and for the
intersection is unitary, while the normal effective permeability for the
intersection is equal to $\widetilde{\lambda}=8$. The source terms are chosen
as Eq. 15 with $\zeta=-1$ for the fractures and $\hat{q}$ for the intersection

$\displaystyle\hat{q}(x,y,z)=10+32y(1-y)\quad\text{in }\gamma.$

We consider a family of triangular meshes and two of coarse meshes. The latter
are constructed from a triangular mesh with different level of coarsening:
$c_{depth}=2$ and $c_{depth}=4$. Fig. 15 shows a graphical representation of
the solution obtained, where also the pressure in the intersection is
reported.

Figure 15: Example of meshes and numerical solutions for the example in
Section 6.2.3. From the left: triangular mesh, coarse mesh with $c_{depth}=2$,
coarse mesh with $c_{depth}=4$, and pressure in $\gamma$ from the
$c_{depth}=2$ case. The solution has range in $[0,1]$ in the fractures and in
$[0,1.25]$ in the intersection, a “Blue to Red Rainbow” colour map is used.

In Fig. 16 we report the $L^{2}$-relative errors for both the projected
velocity $\Pi_{0}\bm{u}$ and the pressure $p$. We observe that the error decay
for the latter is quadratic and the former is linear with respect to the mesh
size, for the three families of mashes.

Figure 16: In the figure the error decay for both $p$, in the left plot, and
$\Pi_{0}\bm{u}$ for the third example in Section 6.2. We add the references
$\mathcal{O}(h)$ and $\mathcal{O}(h^{2})$ to easy the comparison.

In Fig. 17 we report the $L^{2}$-relative errors for both the projected
velocity $\hat{\Pi}_{0}\hat{\bm{u}}$ and the pressure $\hat{p}$ for the
solution in the intersection $\gamma$. The mesh size is now referred to the
intersection mesh. Also in this case the slope of the pressure error is nearly
quadratic for the three families, but also the velocity (especially for the
triangular and coarse-$c_{depth}=4$ families) has an error that behave
quadratically. This is a possible super-convergence property of the numerical
scheme in one space dimension.

Figure 17: In the figure the error decay for both $\hat{p}$, in the left plot,
and $\hat{\Pi}_{0}\hat{\bm{u}}$ for the third example in Section 6.2. We add
the references $\mathcal{O}(h)$ and $\mathcal{O}(h^{2})$ to easy the
comparison.

Table 3 and Table 4 in Appendix B reports the detailed values for this
example, as well as the sparsity of the matrix, the minimum and maximum values
of the solution, number of edges for each cell.

We conclude the section commenting that in all the tests performed the error
decay, for both the pressure and velocity, is coherent with the analyses
presented in [7].

###  6.3 Modeling of intersection flow

Having established the convergence properties of the virtual element method
with and without flow in the intersection, we next focus on the modeling
aspects of 1 and 2. To that end, we consider a network of four fractures,
defined by

$\displaystyle\Omega_{1}=$
$\displaystyle\left\\{(x,y,z)\in\mathbb{R}^{3}:0\leq x\leq 1,y=0,0\leq z\leq
1\right\\}$ $\displaystyle\Omega_{2}=$ $\displaystyle rot\big{(}$
$\displaystyle\left\\{(x,y,z)\in\mathbb{R}^{3}:1/5\leq x-z\leq 6/5,-1/5\leq
x+z\leq 4/5,y=0\right\\},$ $\displaystyle 2\pi/3,(1/2,0,0),(0,0,1)\big{)}$
$\displaystyle\Omega_{3}=$ $\displaystyle rot\big{(}$
$\displaystyle\left\\{(x,y,z)\in\mathbb{R}^{3}:\nu\leq x\leq
1/2-\nu,y=0,1/2+\nu\leq z\leq 1+\nu,\right\\},$
$\displaystyle\pi/6,(1/2,0,1/2),(1,0,-1)\big{)}$ $\displaystyle\Omega_{4}=$
$\displaystyle rot\big{(}$
$\displaystyle\left\\{(x,y,z)\in\mathbb{R}^{3}:1/2+\nu\leq x\leq
1+\nu,y=0,1/2+\nu\leq z\leq 1+\nu,\right\\},$
$\displaystyle\pi/6,(1/2,0,1/2),(1,0,-1)\big{)}.$

Here we have introduced the short-form $\nu=1/(5\sqrt{2})$ and
$rot(\omega,\theta,v,w)$ denotes the rotation of domain $\omega$ with an angle
$\theta$ around the line parametrized as $l(t)=v+wt$. The rotation angles are
immaterial to the solution provided they are non-zero; the values chosen are
motivated by clarity of visualization.

A source of unit strength is assigned in $(1/2,0,1/2)$, that is, in
$\Omega_{1}$. For the other fractures, we assign homogeneous Dirichlet
conditions on the faces in the corners furthest away from the injection
points. On the remaining boundary edges we impose no-flux conditions. In terms
of sources and boundary conditions, the problem is symmetric with respect to
the fractures $\Omega_{2}$, $\Omega_{3}$ and $\Omega_{4}$. We note that, due
to the Neumann condition, this problem is outside the class considered in
Section 3.

In the intersection between $\Omega_{1}$ and $\Omega_{2}$, the effective
normal permeability $\tilde{\lambda}_{m}$ is set to $10^{-7}$. The
intersection between $\Omega_{1}$ and $\Omega_{3}$ has tangential effective
permeability $\hat{\lambda}_{k}=10^{-10}$, effectively employing 1 here, while
for the intersection of $\Omega_{1}$ and $\Omega_{4}$,
$\hat{\lambda}_{k}=10^{10}$. All other permeabilities have unit value.

$\Omega_{2}$$\Omega_{1}$$\Omega_{3}$$\Omega_{4}$

$\Omega_{2}$$\Omega_{1}$$\Omega_{3}$$\Omega_{4}$

Figure 18: Left: Pressure profile obtained for the family of four intersecting
fractures. Right: Cell-wise velocity magnitudes, and also velocity vectors.
Arrows close to the source and the outlets are not shown.

Pressure and velocity profiles are shown in Fig. 18. Due to the low
$\tilde{\lambda}_{m}$, $\Omega_{2}$ is effectively sealed off from the other
fractures, with a pressure that is virtually equal to the outlet value in the
entire fracture. The pressure difference between $\Omega_{3}$ and $\Omega_{4}$
is negligible. However, as seen from the velocity field, flow into
$\Omega_{4}$ is channelized into the intersection. For $\Omega_{3}$ there is
no flow in the intersection, but a wider sweep of the fracture plane itself
compared to $\Omega_{4}$.

###  6.4 Realistic DFN

In this final example, we illustrate the robustness of the approach for
general fracture geometries. To that end, we consider a stochastically
generated network of 60 fractures using the software described in [41]. The
network is shown in Fig. 19. The fracture permeability is set to unity, and
for simplicity we only consider the model presented in 1. Each fracture is
discretized by 8 boundary edges, we impose no-flux boundary condition on all
the edges of the fractures except for two of them. A pressure $p_{1}=1$ and
$p_{2}=0$ is thus imposed. Again we note that the no-flow conditions are not
covered by the analysis in Section 3.

Figure 19: Geometry used for the example in Section 6.4. Each fracture is
coloured by its identification number and the pressure boundary conditions are
highlight using magenta segments.

The stochastic generation of the fracture network leads to intersection
configurations that may pose difficulties for meshing, in particular if the
fracture planes are coupled during grid creation. Our approach of independent
gridding of each fracture avoids the linkage of nodes in different fracture
planes, but it still create small cells close to short constraints and almost
parallel constraints. The coarsening algorithm, here applied with
$c_{depth}=2$ increases the cell size. The polygonal mesh has $24444$ cells
and $78574$ edges, compared with $106809$ cells for the triangular grid. The
impact of the coarsening is illustrated in Fig. 20.

Figure 20: Details of the mesh used in Section 6.4. Left column: Original
triangular grid. Right: Polygonal grid obtained after coarsening.

Figure 21: Solution obtained with a realistic discrete fracture network. The
figure depicts both pressure and velocity, the latter represented by magenta
arrows. In the right the full system is represented while in the left a detail
close to the outflow boundary is shown. The pressure has range in $[0,1]$ and
a “Blue to Red Rainbow” colour map is used. The velocity is restricted to ease
the interpretation and the arrows are scaled with the velocity magnitude.

Fig. 21 shows the pressure solution and the projected velocity field, together
a zoom-in in the vicinity of the outflow. As indicated by the figures, the
numerical method captures high flow velocities in the in- and outflow
fractures, as well as interaction between the fracture planes. It is
worthwhile to consider how other conservative methods could have discretized
the network in Fig. 19. A key property of the virtual element method, and the
closely related mimetic finite difference method, is the ability to handle
hanging nodes and general polygon grids. A discretization by e.g. mixed finite
elements would have required not only the triangular grid with about four
times as many cells, but also techniques to eliminate hanging nodes at the
intersections, which may further increase the cell number. Another possibility
is to consider the class of finite volume methods which are, generally, more
robust with respect to the shape of the cells. See for instance, [13].

##  7 Conclusions

In this paper we presented a novel approximation for a complex network of
fractures using reduced models to describe fracture flows. The reduced models
are a reasonable approximation when the thickness of each fracture is some
orders of magnitude smaller then its other characteristic sizes and smaller
than the typical size of the surrounding rock matrix. In particular two models
are presented to describe the flow at the intersection of fractures, with the
possibility to allow a tangential flow along the intersection. The virtual
element method’s ability to handle general polygons with hanging nodes allows
us to relax constraints on the gridding of fracture intersections, and also
apply grid coarsening. In the examples we saw the error decay of the solution,
both pressure and Darcy velocity, behave as expected in all the cases. Finally
we noticed that the solution behaves as expected and no evidence of serious
contraindications is present, thus the algorithm presented seems a promising
tool for the numerical approximation for this type of problems.

## Acknowledgment

We acknowledge financial support for the ANIGMA project from the Research
Council of Norway (project no. 244129/E20) through the ENERGIX program. The
authors warmly thank the others components of the ANIGMA project team: Eivind
Bastesen, Inga Berre, Simon John Buckley, Casey Nixon, David Peacock, Atle
Rotevatn, Pål Næverlid Sævik, Luisa F. Zuluaga. The authors wish to thank
also: Lourenço Beirão da Veiga, Wietse Boon, Franco Dassi, Luca Formaggia,
Eren Uçar.

##  Appendix A Details of coarsening algorithm

Here we describe the coarsening algorithm in some detail. For further
information, confer [39]. Define two parameters $\epsilon_{str}\in(0,1)$ and
$c_{depth}\in\mathbb{N}$. An element $i$, called $C$-element, in the original
mesh will be connected, in the coarse grid, to some elements in the set
$N_{i}$, called $F$-elements, those affect the error at $i$ most. Since the
algorithm is matrix base the previous request aim to find the $F$-elements $j$
such that $\left|[A]_{ij}\right|$ is the largest in some sense. We say that a
degree of freedom $i$ is “strongly negatively coupled” (SNC) to $j$ if

$\displaystyle-[A]_{ij}\geq\epsilon_{str}\max_{[A]_{ik}<0}\left|[A]_{ik}\right|\quad\text{and
we have}\quad S_{i}\vcentcolon=\left\\{j\in N_{i}:i\text{ is SNC to
}j\right\\}.$

Since the relation of being SNC is not symmetric we introduce
$S_{i}^{\top}\vcentcolon=\left\\{j:i\in S_{j}\right\\}$, which is the set of
elements which are strongly coupled to $i$. Starting from a $C$-element $i$
the algorithm connect all the $F$-elements $j$ that are strongly coupled to
$i$ to create a coarse element. The process is repeated considering another
$C$-element of the mesh. As highlight in [39], the selection of a new
$C$-element $i$ is based on the measure of importance $\lambda_{i}$ in the set
of undecided elements $U$, i.e. nor $C$ or $F$-elements. We define

$\displaystyle\lambda_{i}\vcentcolon=\sharp\left(S_{i}^{\top}\cap
U\right)+2\sharp\left(S_{i}^{\top}\cap F\right)\quad\text{with}\quad i\in U,$

where $F$ is the set of $F$-elements. The next $C$-element will be the element
with greater $\lambda_{i}$. The presented algorithm is repeated for
$c_{depth}$-times where the initial mesh is the coarse mesh obtained at the
previous step. The greater $c_{depth}$ the bigger will be elements in the
final coarse grid.

######  Remark A.1.

As point out in [39] a reasonable value for $\epsilon_{str}$ is $0.25$. We
note also that in the case of permeability with strong anisotropy the final
mesh is able to better represent the underling physic of the problem. See
Section 6.1 for a discussion.

######  Remark A.2.

To avoid pathological shapes of elements during the coarsening process at the
tip of the intersection we impose a priori the triangles, which share the tip
of the intersection and have an edge on the intersection, as $C$-elements.

##  Appendix B Tables

In this section we present the detailed values for the examples in Section 6.
In Table 1 are reported the values for the test of Section 6.2.1. In Table 2
are reported the values for the test of Section 6.2.2.

Coarse mesh family
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  3.162e-1 |  3.340e-2 |  - |  2.713e-1 |  - |  3.434e-3 |  1.228 |  152 |  1.158e-1
2 |  1.581e-1 |  9.774e-3 |  1.773 |  9.639e-2 |  1.493 |  4.958e-4 |  1.2642 |  584 |  3.321e-2
3 |  7.906e-2 |  2.537e-3 |  1.946 |  3.067e-2 |  1.652 |  7.474e-5 |  1.2947 |  2288 |  8.885e-3
4 |  3.953e-2 |  6.348e-4 |  1.998 |  9.360e-3 |  1.712 |  1.049e-5 |  1.3193 |  9056 |  2.292e-3
5 |  1.976e-2 |  1.577e-4 |  2.009 |  2.839e-3 |  1.721 |  1.393e-6 |  1.3308 |  36032 |  5.806e-4
6 |  9.882e-3 |  3.921e-5 |  2.008 |  8.742e-4 |  1.699 |  1.795e-7 |  1.3335 |  143744 |  1.467e-4

Cartesian mesh family
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  1.414e-1 |  4.099e-2 |  - |  1.936e-1 |  - |  -7.991e-4 |  1.283 |  320 |  1.793e-2
2 |  7.071e-2 |  1.061e-2 |  1.949 |  5.923e-2 |  1.709 |  -2.346e-5 |  1.3055 |  1240 |  4.646e-3
3 |  3.535e-2 |  2.682e-3 |  1.984 |  1.715e-2 |  1.788 |  5.599e-7 |  1.3211 |  4880 |  1.234e-3
4 |  1.768e-2 |  6.728e-4 |  1.995 |  4.807e-3 |  1.835 |  2.797e-7 |  1.3291 |  19360 |  3.162e-4
5 |  8.839e-3 |  1.684e-4 |  1.999 |  1.319e-3 |  1.866 |  4.767e-8 |  1.3329 |  77120 |  7.96e-5
6 |  4.419e-4 |  4.211e-5 |  1.999 |  3.567e-4 |  1.887 |  6.738e-9 |  1.335 |  307840 |  2.02e-5

Random mesh family
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  1.171e-1 |  3.4794e-2 |  - |  2.434e-1 |  - |  -8.777e-4 |  1.2988 |  320 |  2.168e-2
2 |  4.675e-2 |  9.8991e-3 |  1.367 |  1.032e-1 |  0.933 |  1.684e-6 |  1.3159 |  1240 |  5.749e-3
3 |  1.975e-2 |  2.5395e-3 |  1.581 |  4.993e-2 |  0.844 |  8.565e-6 |  1.328 |  4880 |  1481e-3
4 |  8.188e-3 |  6.5716e-4 |  1.535 |  2.481e-2 |  0.794 |  1.920e-6 |  1.3335 |  19360 |  3761e-4
5 |  4.04e-3 |  1.6608e-4 |  1.947 |  1.247e-2 |  0.974 |  2.948e-7 |  1.3353 |  77120 |  9.475e-5
6 |  1.843e-3 |  4.2e-5 |  1.752 |  6.218e-3 |  0.887 |  2.4239e-08 |  1.3362 |  307840 |  2.378e-5

Triangular mesh family
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  8.86e-2 |  2.928e-2 |  - |  3.435e-1 |  - |  -1.449e-3 |  1.3089 |  428 |  1.224e-2
2 |  4.221e-2 |  7.238e-3 |  1.885 |  1.572e-1 |  1.024 |  -1.971e-4 |  1.3157 |  1627 |  3.266e-3
3 |  2.029e-2 |  1.775e-3 |  1.919 |  7.426e-2 |  1.024 |  -1.849e-5 |  1.3296 |  6616 |  8.098e-4
4 |  1.014e-2 |  4.265e-4 |  2.055 |  3.597e-2 |  1.045 |  -1.793e-6 |  1.3319 |  26405 |  2.037e-4
5 |  4.594e-3 |  1.078e-4 |  1.738 |  1.807e-2 |  0.869 |  -2.125e-7 |  1.3348 |  105000 |  5.133e-5
6 |  2.204e-3 |  2.771e-5 |  1.849 |  9.167e-3 |  0.924 |  -5.644e-8 |  1.3358 |  408068 |  1.322e-05

Table 1: For each table we report the values for the discretization ($h$),
errors ($err(p)$ and $err({\Pi_{0}\bm{u}})$) and order of convergence
($\mathcal{O}(p)$ and $\mathcal{O}(\Pi_{0}\bm{u})$) for the example in Section
6.2. The last columns are devoted to the minimum and maximum principle, number
of rows (or columns) of the matrix and sparsity. For the name used in each
table, consider the terminology reported in the aforementioned subsection.

Triangular mesh family
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\sharp$ faces |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  2.261e-1 |  3.87e-2 |  - |  3.022e-1 |  - |  3,3,3 |  5.949e-3 |  3.6581 |  646 |  7.96e-3
2 |  1.058e-1 |  8.160e-3 |  2.048 |  1.452e-1 |  0.964 |  3,3,3 |  5.063e-4 |  3.8417 |  3014 |  1.75e-3
3 |  7.461e-2 |  4.209e-3 |  1.898 |  1.082e-1 |  0.843 |  3,3,3 |  5.492e-4 |  3.8612 |  6070 |  8.76e-4
4 |  6.163e-2 |  2.818e-3 |  2.099 |  8.737e-2 |  1.120 |  3,3,3 |  3.625e-4 |  3.8900 |  8902 |  5.99e-4
5 |  3.382e-2 |  8.337e-4 |  2.029 |  4.881e-2 |  0.970 |  3,3,3 |  9.553e-5 |  3.9453 |  29626 |  1.81e-4
6 |  2.404e-2 |  4.235e-4 |  1.984 |  3.447e-2 |  1.019 |  3,3,3 |  3.394e-5 |  3.9620 |  58233 |  9.23e-5

Coarse mesh family with $c_{depth}=2$
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\sharp$ faces |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  4.586e-1 |  9.771e-2 |  - |  3.484e-1 |  - |  5,6,9 |  5.135e-3 |  3.5731 |  268 |  3.659e-2
2 |  2.070e-1 |  2.080e-2 |  1.945 |  1.46e-1 |  1.094 |  5,6,9 |  3.256e-3 |  3.6826 |  1192 |  9.215e-3
3 |  1.474e-1 |  1.11e-2 |  1.85 |  1.016e-1 |  1.066 |  5,6,9 |  7.392e-4 |  3.8070 |  2366 |  4.725e-3
4 |  1.229e-1 |  7.212e-3 |  2.368 |  8.558e-2 |  0.943 |  5,6,9 |  1.016e-3 |  3.8260 |  3422 |  3.304e-3
5 |  6.705e-2 |  2.164e-3 |  1.987 |  4.469e-2 |  1.073 |  5,6,9 |  4.978e-4 |  3.8962 |  11238 |  1.029e-3
6 |  4.775e-2 |  1.093e-3 |  2.011 |  3.168e-2 |  1.013 |  5,6,9 |  1.828e-4 |  3.9359 |  21943 |  5.328e-4

Coarse mesh family with $c_{depth}=4$
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\sharp$ faces |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  4.047e-1 |  2.075e-2 |  - |  8.634e-2 |  - |  17,20,25 |  3.824e-2 |  3.3361 |  758 |  3.960e-2
2 |  2.686e-1 |  9.821e-3 |  1.825 |  5.670e-2 |  1.026 |  13,19,24 |  9.755e-3 |  3.6202 |  1528 |  2.002e-2
3 |  1.973e-1 |  5.200e-3 |  2.06 |  3.255e-2 |  1.798 |  14,20,26 |  6.993e-3 |  3.6896 |  2792 |  1.168e-2
4 |  1.392e-1 |  2.525e-3 |  2.072 |  2.540e-2 |  0.711 |  13,20,26 |  2.891e-3 |  3.7927 |  5354 |  6.346e-3
5 |  9.758e-2 |  1.198e-3 |  2.099 |  1.569e-2 |  1.357 |  14,20,27 |  1.282e-3 |  3.8650 |  10537 |  3.283e-3
6 |  6.917e-2 |  6.154e-4 |  1.936 |  1.121e-2 |  0.976 |  13,20,26 |  6.152e-4 |  3.9082 |  20554 |  1.734e-3

Coarse mesh family with $c_{depth}=5$
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\sharp$ faces |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  7.449e-1 |  4.609e-2 |  - |  2.371e-1 |  - |  29,38,48 |  9.882e-2 |  2.4769 |  520 |  9.531e-2
2 |  4.761e-1 |  2.827e-2 |  1.091 |  1.538e-1 |  0.967 |  28,36,44 |  6.912e-2 |  3.1075 |  890 |  5.795e-2
3 |  3.429e-1 |  1.313e-2 |  2.336 |  8.447e-2 |  1.825 |  18,35,49 |  2.227e-2 |  3.4482 |  1644 |  3.221e-2
4 |  2.515e-1 |  7.052e-3 |  2.005 |  5.368e-2 |  1.462 |  25,37,47 |  2.055e-2 |  3.5211 |  2942 |  2.003e-2
5 |  1.773e-1 |  3.138e-3 |  2.316 |  3.355e-2 |  1.345 |  21,37,49 |  9.645e-3 |  3.6361 |  5725 |  1.056e-2
6 |  1.27e-1 |  1.630e-3 |  1.962 |  2.577e-2 |  0.79 |  24,37,44 |  2.642e-3 |  3.8138 |  10908 |  5.841e-3

Table 2: For each table we report the values for the discretization ($h$),
errors ($err(p)$ and $err({\Pi_{0}\bm{u}})$) and order of convergence
($\mathcal{O}(p)$ and $\mathcal{O}(\Pi_{0}\bm{u})$) for the example in Section
6.2.2. The last columns are devoted to the number of faces for each cell
(minimum, average, and maximum, respectively), the minimum and maximum
principle, number of rows of the matrix and sparsity. For the name used in
each table, consider the terminology reported in the aforementioned
subsection.

Triangular mesh family
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\sharp$ faces |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  2.261e-1 |  5.494e-2 |  - |  3.112e-1 |  - |  3,3,3 |  6.051e-3 |  9.139e-1 |  655 |  7.871e-3
2 |  1.058e-1 |  1.249e-2 |  1.95 |  1.497e-1 |  0.963 |  3,3,3 |  5.141e-4 |  9.606e-1 |  3031 |  1.746e-3
3 |  6.163e-2 |  3.807e-3 |  2.2 |  8.795e-2 |  0.985 |  3,3,3 |  3.632e-4 |  9.723e-1 |  8935 |  5.974e-4
4 |  3.382e-2 |  1.082e-3 |  2.1 |  4.784e-2 |  1.01 |  3,3,3 |  9.558e-5 |  9.809e-1 |  29691 |  1.807e-4
5 |  1.695e-2 |  2.722e-4 |  2 |  2.434e-2 |  0.978 |  3,3,3 |  3.867e-5 |  9.913e-1 |  117311 |  4.588e-5

Coarse mesh family with $c_{depth}=2$
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\sharp$ faces |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  4.586e-1 |  9.705e-2 |  - |  2.802e-1 |  - |  5,6,8 |  4.628e-3 |  8.183e-1 |  277 |  3.499e-2
2 |  2.070e-1 |  2.211e-2 |  1.86 |  1.393e-1 |  0.88 |  5,6,9 |  3.253e-3 |  9.152e-1 |  1209 |  9.036e-3
3 |  1.058e-1 |  5.478e-3 |  2.08 |  6.718e-2 |  1.09 |  5,6,9 |  3.466e-4 |  9.638e-1 |  4573 |  2.466e-3
4 |  5.461e-2 |  1.246e-3 |  2.24 |  3.519e-2 |  0.977 |  5,6,9 |  5.242e-4 |  9.777e-1 |  16761 |  6.928e-4
5 |  3.033e-2 |  4.337e-4 |  1.79 |  1.928e-2 |  1.02 |  5,6,10 |  5.833e-5 |  9.889e-1 |  54617 |  2.147e-4

Coarse mesh family with $c_{depth}=4$
---
$\sharp$ |  $h$ |  $err(p)$ |  $\mathcal{O}(p)$ |  $err({\Pi_{0}\bm{u}})$ |  $\mathcal{O}(\Pi_{0}\bm{u})$ |  $\sharp$ faces |  $\min p$ |  $\max p$ |  size |  sparsity
1 |  2.686e-1 |  1.249e-2 |  - |  6.404e-2 |  - |  13,19,24 |  9.662e-3 |  8.496e-1 |  1561 |  1.927e-2
2 |  1.253e-1 |  2.592e-3 |  2.06 |  2.509e-2 |  1.23 |  14,20,27 |  3.294e-3 |  9.399e-1 |  6649 |  5.112e-3
3 |  8.399e-2 |  1.107e-3 |  2.13 |  1.489e-2 |  1.31 |  14,20,30 |  7.368e-4 |  9.630e-1 |  14385 |  2.438e-3
4 |  6.012e-2 |  5.624e-4 |  2.02 |  9.809e-3 |  1.25 |  12,20,17 |  7.374e-4 |  9.797e-1 |  27537 |  1.294e-3
5 |  4.013e-2 |  2.352e-4 |  2.16 |  6.26e-3 |  1.11 |  13,20,27 |  1.825e-4 |  9.831e-1 |  60851 |  5.944e-4

Table 3: For each table we report the values for the discretization ($h$),
errors ($err(p)$ and $err({\Pi_{0}\bm{u}})$) and order of convergence
($\mathcal{O}(p)$ and $\mathcal{O}(\Pi_{0}\bm{u})$) for the example in Section
6.2.3. The last columns are devoted to the number of faces for each cell
(minimum, average, and maximum, respectively), the minimum and maximum
principle, number of rows of the matrix and sparsity. For the name used in
each table, consider the terminology reported in the aforementioned
subsection.

Triangular mesh family
---
$\sharp$ |  $h$ |  $err(\hat{p})$ |  $\mathcal{O}(\hat{p})$ |  $err({\hat{\Pi}_{0}\hat{\bm{u}}})$ |  $\mathcal{O}(\hat{\Pi}_{0}\hat{\bm{u}})$ |  $\min\hat{p}$ |  $\max\hat{p}$
1 |  1.250e-1 |  4.665e-2 |  - |  1.793e-2 |  - |  3.423e-1 |  1.269
2 |  6.250e-2 |  1.158e-2 |  2.01 |  5.010e-3 |  1.84 |  1.647e-1 |  1.254
3 |  3.125e-2 |  3.229e-3 |  1.84 |  9.500e-4 |  2.4 |  8.042e-2 |  1.252
4 |  1.562e-2 |  8.601e-4 |  1.91 |  1.812e-4 |  2.39 |  3.966e-2 |  1.250
5 |  7.812e-3 |  2.187e-4 |  1.97 |  4.028e-5 |  2.17 |  1.968e-2 |  1.250

Coarse mesh family with $c_{depth}=2$
---
$\sharp$ |  $h$ |  $err(\hat{p})$ |  $\mathcal{O}(\hat{p})$ |  $err({\hat{\Pi}_{0}\hat{\bm{u}}})$ |  $\mathcal{O}(\hat{\Pi}_{0}\hat{\bm{u}})$ |  $\min\hat{p}$ |  $\max\hat{p}$
1 |  1.250e-1 |  7.198e-2 |  - |  8.557e-3 |  - |  3.524e-1 |  1.303
2 |  6.250e-2 |  1.769e-2 |  2.02 |  2.143e-3 |  2 |  1.658e-1 |  1.262
3 |  3.125e-2 |  4.511e-3 |  1.97 |  7.088e-4 |  1.6 |  8.061e-2 |  1.253
4 |  1.538e-2 |  1.246e-3 |  1.82 |  3.032e-4 |  1.2 |  3.967e-2 |  1.251
5 |  7.812e-3 |  3.541e-4 |  1.86 |  1.124e-4 |  1.47 |  1.968e-2 |  1.250

Coarse mesh family with $c_{depth}=4$
---
$\sharp$ |  $h$ |  $err(\hat{p})$ |  $\mathcal{O}(\hat{p})$ |  $err({\hat{\Pi}_{0}\hat{\bm{u}}})$ |  $\mathcal{O}(\hat{\Pi}_{0}\hat{\bm{u}})$ |  $\min\hat{p}$ |  $\max\hat{p}$
1 |  3.125e-2 |  2.717e-3 |  - |  5.785e-3 |  - |  8.066e-2 |  1.246
2 |  1.408e-2 |  9.627e-4 |  1.3 |  1.703e-3 |  1.53 |  3.965e-2 |  1.248
3 |  7.812e-3 |  3.718e-4 |  1.61 |  5.831e-4 |  1.82 |  1.968e-2 |  1.249
4 |  6.135e-3 |  2.059e-4 |  2.45 |  3.448e-4 |  2.17 |  1.968e-2 |  1.25
5 |  3.906e-3 |  8.863e-5 |  1.87 |  1.411e-4 |  1.98 |  9.803e-3 |  1.25

Table 4: For each table we report the values for the discretization ($h$),
errors ($err(\hat{p})$ and $err({\hat{\Pi}_{0}\hat{\bm{u}}})$) and order of
convergence ($\mathcal{O}(\hat{p})$ and
$\mathcal{O}(\hat{\Pi}_{0}\hat{\bm{u}})$) for the example in Section 6.2.3.
The mesh size is referred to the intersection mesh. For the name used in each
table, consider the terminology reported in the aforementioned subsection.

## References

  * [1] C. Alboin, J. Jaffré, J. E. Roberts, X. Wang, and C. Serres, Domain Decomposition for some Transmission Problems in Flow in Porous Media, in Numerical treatment of multiphase flows in porous media (Beijing, 1999), vol. 552 of Lecture Notes in Phys., Springer, Berlin, 2000, pp. 22–34, doi:10.1007/3-540-45467-5_2.
  * [2] L. Amir, M. Kern, V. Martin, and J. E. Roberts, Décomposition de domaine et préconditionnement pour un modèle 3D en milieu poreux fracturé, in Proceeding of JANO 8, 8th conference on Numerical Analysis and Optimization, Dec. 2005.  2005\.
  * [3] P. F. Antonietti, L. B. da Veiga, D. Mora, and M. Verani, A stream virtual element formulation of the Stokes problem on polygonal meshes, SIAM Journal on Numerical Analysis, 52 (2014), pp. 386–404, doi:10.1137/13091141X.
  * [4] J. Bear, C.-F. Tsang, and G. de Marsily, Flow and contaminant transport in fractured rock, Academic Press, San Diego, 1993.
  * [5] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo, H(div) and H(curl)-conforming VEM, ArXiv e-prints, (2014), arXiv:1407.6822.
  * [6] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo, The hitchhiker’s guide to the virtual element method, Mathematical Models and Methods in Applied Sciences, 24 (2014), pp. 1541–1573, doi:10.1142/S021820251440003X.
  * [7] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo, Mixed virtual element methods for general second order elliptic problems on polygonal meshes, ESAIM: M2AN, 50 (2016), pp. 727–747, doi:10.1051/m2an/2015067.
  * [8] L. Beirão da Veiga, F. Brezzi, and D. L. Marini, Virtual elements for linear elasticity problems, SIAM Journal on Numerical Analysis, 51 (2013), pp. 794–812, doi:10.1137/120874746, http://dx.doi.org/10.1137/120874746, arXiv:http://dx.doi.org/10.1137/120874746.
  * [9] M. F. Benedetto, S. Berrone, A. Borio, S. Pieraccini, and S. Scialò, A hybrid mortar virtual element method for discrete fracture network simulations, Journal of Computational Physics, 306 (2016), pp. 148 – 166, doi:http://dx.doi.org/10.1016/j.jcp.2015.11.034.
  * [10] M. F. Benedetto, S. Berrone, S. Pieraccini, and S. Scialò, The virtual element method for discrete fracture network simulations, Computer Methods in Applied Mechanics and Engineering, 280 (2014), pp. 135–156, doi:http://dx.doi.org/10.1016/j.cma.2014.07.016.
  * [11] S. Berrone, S. Pieraccini, and S. Scialò, On simulations of discrete fracture network flows with an optimization-based extended finite element method, SIAM Journal on Scientific Computing, 35 (2013), pp. 908–935, doi:10.1137/120882883.
  * [12] W. M. Boon and J. M. Nordbotten, Robust discretization of flow in fractured porous media, ArXiv e-prints, (2016), arXiv:1601.06977.
  * [13] K. Brenner, J. Hennicker, R. Masson, and P. Samier, Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces, IMA Journal of Numerical Analysis, (2016), doi:10.1093/imanum/drw044.
  * [14] F. Brezzi, R. S. Falk, and D. L. Marini, Basic principles of mixed virtual element methods, ESAIM: M2AN, 48 (2014), pp. 1227–1240, doi:10.1051/m2an/2013138.
  * [15] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, vol. 15 of Computational Mathematics, Springer Verlag, Berlin, 1991\.
  * [16] C. D’Angelo and A. Scotti, A mixed finite element method for Darcy flow in fractured porous media with non-matching grids, Mathematical Modelling and Numerical Analysis, 46 (2012), pp. 465–489, doi:10.1051/m2an/2011148.
  * [17] G. Dziuk, Finite Elements for the Beltrami operator on arbitrary surfaces, in Partial Differential Equations and Calculus of Variations, S. Hildebrandt and R. Leis, eds., vol. 1357 of Lecture Notes in Mathematics, Springer Berlin Heidelberg, 1988, pp. 142–155, doi:10.1007/BFb0082865.
  * [18] B. Faybishenko, P. A. Witherspoon, and S. M. Benson, Dynamics of fluids in fractured rock, vol. 122 of Geophysical monographs, American geophysical union, 2000.
  * [19] B. Flemisch, A. Fumagalli, and A. Scotti, A Review of the XFEM-Based Approximation of Flow in Fractured Porous Media, vol. 12 of SEMA SIMAI Springer Series, Springer International Publishing, Cham, 2016, ch. Advances in Discretization Methods, pp. 47–76, doi:10.1007/978-3-319-41246-7_3.
  * [20] L. Formaggia, A. Fumagalli, A. Scotti, and P. Ruffo, A reduced model for Darcy’s problem in networks of fractures, ESAIM: Mathematical Modelling and Numerical Analysis, 48 (2014), pp. 1089–1116, doi:10.1051/m2an/2013132.
  * [21] N. Frih, J. E. Roberts, and A. Saada, Modeling fractures as interfaces: a model for Forchheimer fractures, Computers and Geosciences, 12 (2008), pp. 91–104, doi:10.1007/s10596-007-9062-x.
  * [22] A. Fumagalli and A. Scotti, An Efficient XFEM Approximation of Darcy Flows in Arbitrarily Fractured Porous Media, Oil and Gas Sciences and Technologies - Revue d’IFP Energies Nouvelles, 69 (2014), pp. 555–564, doi:10.2516/ogst/2013192.
  * [23] P. Grisvard, Singularities in Boundary Value Problems, vol. 22 of Recherches en mathématiques appliquées, Springer-Verlag, 1992, https://books.google.no/books?id=Kjqj-bAx4PYC.
  * [24] A. Hansbo and P. Hansbo, An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 5537–5552, doi:10.1016/S0045-7825(02)00524-8.
  * [25] M.-H. Hui, B. T. Mallison, and K.-T. Lim, An innovative workflow to model fractures in a giant carbonate reservoir, in International Petroleum Technology Conference, 3-5 December, Kuala Lumpur, Malaysia, 2008, doi:doi:10.2523/12572-MS.
  * [26] J. D. Hyman, C. W. Gable, S. L. Painter, and N. Makedonska, Conforming Delaunay Triangulation of Stochastically Generated Three Dimensional Discrete Fracture Networks: A Feature Rejection Algorithm for Meshing Strategy, SIAM Journal on Scientific Computing, 36 (2014), pp. A1871–A1894, doi:10.1137/130942541.
  * [27] J. D. Hyman, S. Karra, N. Makedonska, C. W. Gable, S. L. Painter, and H. S. Viswanathan, dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport, Computers & Geosciences, 84 (2015), pp. 10–19, doi:10.1016/j.cageo.2015.08.001.
  * [28] M. Karimi-Fard and L. Durlofsky, A general gridding, discretization, and coarsening methodology for modeling flow in porous formations with discrete geological features, Advances in Water Resources, 96 (2016), pp. 354–372, doi:10.1016/j.advwatres.2016.07.019.
  * [29] M. Karimi-Fard, B. Gong, and L. J. Durlofsky, Generation of coarse-scale continuum flow models from detailed fracture characterizations, Water Resources Research, 42 (2006), doi:10.1029/2006WR005015, http://dx.doi.org/10.1029/2006WR005015.
  * [30] P. Knabner and J. E. Roberts, Mathematical analysis of a discrete fracture model coupling Darcy flow in the matrix with Darcy-Forchheimer flow in the fracture, ESAIM: Mathematical Modelling and Numerical Analysis, 48 (2014), pp. 1451–1472, doi:10.1051/m2an/2014003.
  * [31] K.-A. Lie, S. Krogstad, I. S. Ligaarden, J. R. Natvig, H. M. Nilsen, and B. Skaflestad, Open-source MATLAB implementation of consistent discretisations on complex grids, Computational Geosciences, 16 (2012), pp. 297–322, doi:10.1007/s10596-011-9244-4.
  * [32] V. Martin, J. Jaffré, and J. E. Roberts, Modeling Fractures and Barriers as Interfaces for Flow in Porous Media, SIAM J. Sci. Comput., 26 (2005), pp. 1667–1691, doi:10.1137/S1064827503429363.
  * [33] H. Mustapha and K. Mustapha, A New Approach to Simulating Flow in Discrete Fracture Networks with an Optimized Mesh, SIAM Journal on Scientific Computing, 29 (2007), pp. 1439–1459, doi:10.1137/060653482.
  * [34] D. C. Peacock, C. W. Nixon, A. Rotevatn, D. J. Sanderson, and L. F. Zuluaga, Interacting faults, Journal of Structural Geology, 97 (2017), pp. 1 – 22, doi:http://dx.doi.org/10.1016/j.jsg.2017.02.008.
  * [35] S. Pieraccini and S. Scialò, Advances in Discretization Methods, vol. 12 of SEMA-SIMAI Springer Series on Applied and Industrial Mathematics, Springer, 2016, ch. On a PDE-Constrained Optimization Approach for Flow Simulations in Fractured Media, pp. 27–45, doi:10.1007/978-3-319-41246-7.
  * [36] A. Rotevatn, S. J. Buckley, J. A. Howell, and H. Fossen, Overlapping faults and their effect on fluid flow in different reservoir types: A LIDAR-based outcrop modeling and flow simulation study, AAPG Bulletin, (2009), doi:10.1306/09300807092.
  * [37] N. Schwenck, B. Flemisch, R. Helmig, and B. Wohlmuth, Dimensionally reduced flow models in fractured porous media: crossings and boundaries, Computational Geosciences, 19 (2015), pp. 1219–1230, doi:10.1007/s10596-015-9536-1.
  * [38] J. R. Shewchuk, Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator, in Applied Computational Geometry: Towards Geometric Engineering, M. C. Lin and D. Manocha, eds., vol. 1148 of Lecture Notes in Computer Science, Springer-Verlag, May 1996, pp. 203–222.  From the First ACM Workshop on Applied Computational Geometry.
  * [39] U. Trottenberg, C. W. Oosterlee, and A. Schüller, Multigrid, Elsevier Academic Press, 2001.
  * [40] C. Xiang, X. Pang, W. Yang, J. Wang, Q. Li, L. Liu, and Y. Li, Hydrocarbon migration and accumulation along the fault intersection zone—a case study on the reef-flat systems of the no.1 slope break zone in the tazhong area, tarim basin, Petroleum Science, 7 (2010), pp. 211–225, doi:10.1007/s12182-010-0021-0.
  * [41] C. Xu and P. Dowd, A new computer code for discrete fracture network modelling, Computers & Geosciences, 36 (2010), pp. 292 – 301, doi:https://doi.org/10.1016/j.cageo.2009.05.012.