# A novel equi-dimensional finite element method for
flow and transport in fractured porous media satisfying
discrete maximum principle and conservation properties

Maria Giuseppina Chiara Nestola<EMAIL_ADDRESS>Marco Favino
<EMAIL_ADDRESS>Institute of Computational Science, Center for
Computational Medicine in Cardiology (CCMC), Università della Svizzera
italiana, Via Giuseppe Buffi 13, 9600 Lugano, Switzerland.  Institute of
Geochemistry and Petrology, ETH Zurich Clausiusstrasse 25, NW F72 8092 Zurich,
Switzerland.  Institute of Earth Science, University of Lausanne, Géopolis -
CH-1015 Lausanne, Switzerland.

###### Abstract

Numerical simulations of flow and transport in porous media usually rely on
hybrid-dimensional models, i.e., the fracture is considered as objects of a
lower dimension compared to the embedding matrix. Such models are usually
combined with non-conforming discretizations as they avoid the inherent
difficulties associated with the generation of meshes that explicitly resolve
fractures-matrix interfaces. However, non-conforming discretizations demand a
more complicated coupling of different sub-models and may require special care
to ensure conservative fluxes. We propose a novel approach for the simulation
of flow and transport problems in fractured porous media based on an equi-
dimensional representation of the fractures. The major challenge for these
types of representation is the creation of meshes which resolve the several
complex interfaces between the fractures and the embedding matrix. To overcome
this difficulty, we employ a strategy based on adaptive mesh refinement (AMR).
The idea at the base of the proposed AMR is to start from an initially uniform
coarse mesh and refine the elements which have non-empty overlaps with at
least one of the fractures. Iterating this process allows to create non-
uniform non-conforming meshes, which do not resolve the interfaces but can
approximate them with arbitrary accuracy. We demonstrate that low-order finite
element (FE) discretizations on adapted meshes are globally and locally
conservative and we suitably adapt an algebraic flux correction technique to
ensure the discrete maximum principle. In particular, we show that the
notorious conditions on M-matrices have to be adapted to the basis functions
defined on non-conforming meshes. Although the proposed applications come from
geophysical applications, the obtained results could be applied to any
diffusion and transport problems, on both conforming and non-conforming
meshes.

###### keywords:

Fractured porous media; Coupled flow and transport problems; Equi-dimensional
model; Adaptive mesh refinement; Discrete maximum principle; Conservative
fluxes.

††journal: Journal of LaTeX Templates

##  1 Introduction

Numerical simulations of flow and transport problems in fractured porous media
are fundamental for numerous geophysical applications, such as geothermal
energy production, hydrocarbon exploration, nuclear waste disposal, and CO2
storage [1, 2, 3, 4, 5, 6]. Fractures can be regarded as heterogeneities
inside an embedding matrix and can create preferential paths and/or barriers
for flow and transport. Fractures can be arranged in complicated networks and
their aperture is several orders of magnitude smaller than the characteristic
size of the embedding matrix. For these reasons, models that represent the
fractures as inclusions with the same geometrical dimension as the embedding
matrix (referred to as equi-dimensional models) are usually considered
unfeasible [7]. In particular, the creation of meshes which explicitly resolve
the complex interfaces between fractures and the embedding matrix is one of
the main difficulties and represents one of the major bottlenecks in the use
of equi-dimensional models.

For overcoming such difficulties, alternative approaches based on implicit
and/or lower-dimensional representations of the fractures have been derived.
Implicit representations consider the effect of the fractures as material
parameters in an effective continuum. This class includes single-continuum
approaches [8, 9] and multi-continuum approaches where governing equations are
defined for each block and are coupled by interaction terms [10]. Lower-
dimensional representations give rise to hybrid-dimensional approaches which
model fractures as objects of a dimension lower than the one of the embedding
matrix and the fracture aperture is treated as a parameter. The most common
approach currently employed is based on the so-called discrete fracture-matrix
(DFM) models, which implicitly treats the smaller fractures and the larger
fractures are usually considered in a hybrid-dimensional approach.

The numerical simulation of both equi- and hybrid-dimensional approaches may
make use of either conforming or non-conforming meshes. The generation of
conforming meshes is particularly complicated and difficult to make automatic
for realistic fracture networks [11, 12, 13, 14, 15, 16].

On the other hand, non-conforming meshes have to be employed in combination
with specific discretization strategies in order to couple the equations
defined on the matrix domain and on the fracture domain such as the unfitted
finite element methods [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29].

For what concerns flow problems in fractured porous media, they have been
approximated by finite volume (FV) [11, 12, 13, 14], finite differences (FD)
[30, 31], virtual finite element (VFE) [32], and finite element (FE) methods
[33, 24, 34]. FV method is well-known to ensure local conservation properties
by construction, i.e., ensure the balance of fluxes between two arbitrary,
complementary subdomains. Conservation properties for FE method hold naturally
for mixed formulations based on Raviart-Thomas discretizations [33, 35, 36,
37, 38, 19, 20], or discontinuous Galerkin (DG) methods [39], but they can
also be proven for enriched continuous Galerkin methods [40, 34], the extended
finite elements (XFEM) [24, 25, 26], and the embedded finite element methods
(EFEM) [23]. Moreover, the virtual elements method (VEM) and mimetic finite
difference (MFD) have been adapted [30, 32] to preserve conservation
properties even in the presence of a highly anisotropic and distorted grid.
For FE methods based on continuous Galerkin (CG) formulations, the computation
of the approximated fluxes relies, in general, on post-processing techniques
of the numerical solution [41, 42, 43, 44] and the resulting fluxes are as
accurate as the ones computed employing discretizations of the same order
based on mixed formulations [43]. In [45], the authors demonstrated that low-
order FE discretizations are conservative by employing an equivalent
formulation of the problem where an auxiliary flux is associated with the
Dirichlet boundary conditions has been introduced. This alternative
formulation allowed to prove that continuous FE formulations are conservative
for any patch of elements in case of conforming meshes.

The transport of a concentration through a fractured porous medium is
described by a pure advection problem, whose velocity field is computed
through the solution of a flow problem. The solution of a transport problem
can be proven to respect the maximum principle and, hence, to be positive for
positive boundary and initial conditions. On the other hand, a violation of
the discrete maximum principle (DMP) for the solutions of discretizations of
transport equations manifests itself in spurious oscillations and negative
values. Satisfying the DMP is, hence, a key requirement when dealing with
transport problems and, more in general, with differential problems involving
positive physical quantities, such as concentrations. Moreover,
discretizations whose solution respect the DMP can be proven to converge
uniformly [46].

Discretizations that provide positive solutions are usually characterized by
stiffness matrices which are monotone or M-matrices. In particular low-order
FE discretizations are characterized by matrices that have positive diagonal
entries and non-negative extra-diagonal entries. A typical example of matrices
that respect these properties come from the discretization of Laplace
operators on grids whose elements do not present obtuse angles. The FE
discretization based on CG formulations of transport problems typically gives
rise to stiffness matrices which violate the DMP. This is due to the presence
of positive extra-diagonal entries in the mass matrix and in the discrete
advection operator and of negative diagonal entries in the discrete advection
operator.

Several approaches have been introduced to try to ensure the DMP and preserve
the positivity of the solution such as 1) artificial diffusion [47], 2)
stream-line diffusion [48], 3) residual-free bubbles [49], 4) local projection
stabilizations [50]. The idea at the base of all those methods is to introduce
an additional diffusion term to the formulation of the transport problem. On
the other hand, algebraic flux correction [51] introduces an algebraic
diffusion operator to compensate for the entries which does not allow to
ensure the DMP. Instead, FV discretizations naturally provide monotone
matrices when employed with upwind schemes.

Most of the discretization methods for flow and transport problems are based
on hybrid-dimensional representations, while few are based on equi-dimensional
representations. Equi-dimensional representations of fractures are less
popular mainly due to the complexity of meshing complicated fractures networks
[30]. The creation of meshes that resolve the interfaces between fractures and
the embedding matrix is a time-consuming and difficult-to-make-automatic
process. Moreover, the resulting meshes are characterized by elongated and/or
distorted elements, due to the small aperture of the fractures. On the other
hand, hybrid dimensional approaches require the development of 1) specific
strategies to couple the equations defined on the fracture domain and matrix
domain to impose the continuity of the unknowns at the interface; 2) specific
algorithms to detect the intersection between the elements of the fracture and
the matrix domain; 3) solutions strategies designed to treat the coupling
terms [24, 26, 23, 52]. Moreover, hybrid-dimensional models neglect cross
fracture phenomena, do not always guarantee local conservation at the matrix-
fracture interfaces [53], may not ensure the DMP.

The focus of this work is to propose a discretization method for the
simulation of the coupled flow and transport problems in fractured media based
on an equi-dimensional representation of the fractures. For creating suitable
meshes, we employ the adaptive mesh refinement (AMR) strategy proposed in
[54]. The idea of this strategy is to start from a uniform coarse mesh that is
not related to the fracture distribution and to refine the elements which have
a non-empty overlap with at least one fracture. Iterating this step allows
generating adapted non-conforming meshes which are refined at the interfaces
between the fractures and the embedding matrix. The resulting meshes are
composed of squared and cubic elements which reduce the ill-conditioning due
to the presence of skewed elements, at the price of having non-constant
material properties over some elements. This approach allows to explicitly
take fractures into account by assigning different material properties during
the assembly of the stiffness matrices.

For transport problems with shock waves, linear and non-linear algebraic
stabilizations have been compared on adapted meshes to preserve monotonicity
on adapted meshes [55]. Differently, from our geometric approach, the authors
employed a specific error indicator to select the elements to refine.

We show that the proposed discretization is conservative on non-conforming
adapted meshes and for any subdomain. In combination with an algebraic flux
correction scheme, such discretization ensures the DMP. Moreover, the CG
formulation does not present the matrix-fracture coupling terms and, hence,
produces linear systems with a positive definite matrix.

The paper is organized as follows. In Section 2, we discuss the DMP on
conforming meshes and introduce the governing equations for the coupled flow
and transport problems. In Section 3, we introduce the FE formulation and
discuss the DMP on adapted (non-conforming) meshes while in Section 4 we study
the conservation properties of the proposed approach. In Section 5, we present
some numerical results, including a real problem with a complex stochastic
fracture network. Finally, we make some concluding remarks in Section 6.

##  2 Discrete maximum principle on conforming meshes

We introduce the FE discretizations of the coupled models which describe the
single-phase flow and the advective transport of a solute in a fractured
porous medium. and their relevant properties in terms of continuum and
discrete maximum principle. We employ models based on an _equi-dimensional_
representation of the fractures [56, 57], i.e., fractures are described as
$d$-dimensional objects embedded in a $d$-dimensional porous matrix. Equi-
dimensional representations do not involve couplings between the
discretizations of different dimensional objects. Hence, they allow for a
natural extension of classical results related to FE methods based on CG
formulations.

###  2.1 Preliminaries definitions

We denote an axis-aligned box domain in $\mathbb{R}^{d}$, i.e., a rectangle
for $d=2$ or a cuboid for $d=3$, by $\Omega$, where $d$ denotes the dimension
of the problem. The length along the $i$-axis is denoted by $L_{i}$. We denote
by $\partial\Omega$ the boundary of $\Omega$. The outward unit vector to
$\Omega$ is denoted by $\underline{n}$. Points of $\Omega$ are denoted by
capital non-bold letters $X$ and their coordinates by
$\mathbf{x}=(x_{1},x_{2},\ldots,x_{d})$.

Over $\Omega$, we introduce a set $\mathcal{F}$ of $d$-dimensional inclusions
$f_{i}\subset\Omega$ with $i=1,2\,\ldots n^{\mathcal{F}}$, where
$n^{\mathcal{F}}=|\mathcal{F}|$ is the number of inclusions. We assume
inclusions are non-axis-aligned boxes. We are particularly interested in the
case where inclusions have one dimension that is much smaller compared to the
other(s). In this case, we will refer to such inclusions as fractures and the
smaller dimension will be referred to as aperture or thickness.

We define $\Omega_{f}:=\bigcup_{i=1}^{n^{\mathcal{F}}}f_{i}$, i.e., the sets
of the inclusions and $\Omega_{m}:=\Omega\backslash\Omega_{f}$, i.e., the set
representing the matrix. Both $\Omega_{m}$ and $\Omega_{f}$ may be non-
connected sets. The set
$\Gamma:=\overline{\Omega}_{m}\cap\overline{\Omega}_{f}$ is the interface
between the two subdomains.

We denote the fracture aperture by $\delta$, the permeability by $k$, and the
porosity by $\phi$. Permeability and porosity are assumed to be real functions
attaining constant positive values over $\Omega_{m}$ and $\Omega_{f}$, such
as, for example,

$k=\begin{cases}k_{m},\quad\text{ in }\Omega_{m},\\\ k_{f},\,\,\quad\text{ in
}\Omega_{f}.\\\ \end{cases}$ (1)

It is worth to point out that the assumptions adopted to introduce the reader
to the topic could be relaxed: $\Omega$ has not to be a box domain; fractures
may have variable apertures and different shapes; the material properties can
be heterogeneous. In particular, heterogeneous properties would require to
replace scalar values by tensor quantities.

We suppose that $\partial\Omega$ can be decomposed in two subsets $\Gamma_{D}$
and $\Gamma_{N}$ with
$\partial\Omega=\overline{\Gamma}_{D}\cup\overline{\Gamma}_{N}$ and
$\mathring{\Gamma}_{D}\cap\mathring{\Gamma}_{N}=\emptyset$. The set
$\Gamma_{D}$ denotes the subset of $\partial\Omega$ where Dirichlet boundary
conditions are imposed. Notice that $\Gamma_{D}$ may differ for flow and
transport problems.

For a set $\Xi$, we denote by the set of square-integrable functions by
$L^{2}(\Xi)$. The scalar product and the norm in $L^{2}(\Xi)$ are
$\|\,\cdot\,\|_{\Xi}$, and $(\,\cdot\,,\,\cdot\,)_{\Xi}$, respectively. Given
an unknown $q$ and Dirichlet boundary conditions $q_{D}$ imposed on
$\Gamma_{D}$, the set $S$ denotes a suitable function space where the solution
of each problem is sought. The space $S$ will be specified in Section 3 for
each of the problems. Moreover, we introduce the following function spaces:

$U=\\{v\in S\,:\,v=q_{D}\text{ on }\Gamma_{D}\\},$ (2)

and

$V=\\{v\in S\,:\,v=0\text{ on }\Gamma_{D}\\}.$ (3)

Over $\Omega$, we introduce a mesh $\mathcal{T}$ with $N_{E}$ elements and
$N_{N}$ nodes. Elements can be triangles or quadrilaterals for $d=2$ and
tetrahedra or hexahedra for $d=3$. We assume that all elements of a mesh are
of the same type. Elements are denoted by $E$ and have vertices, edges, and,
for $d=3$, faces. The term side refers to edges for $d=2$ and to faces for
$d=3$. We assume, for now, that the mesh is conforming, i.e., the intersection
of each pair of elements is empty, a vertex, an edge, or a face.

We let ${S}^{h}$ denote the nodal interpolation space over $\mathcal{T}$,
i.e.,

${S}^{h}=\\{v\in C^{0}(\Omega)\,:\,v|_{E}\in\mathbb{L}_{1}\\},$

where $\mathbb{L}_{1}$ is the space of linear functions on $E$ for triangular
and tetrahedral meshes (usually denoted by $\mathbb{P}_{1}$) or the space of
multi-linear functions on $E$ for quadrilateral or hexahedral meshes (usually
denoted by $\mathbb{Q}_{1}$).

We let $J$ denote the set of all nodal indices $i=1,2,\ldots,N_{N}$ and we let
$D$ be the subset corresponding to nodes located in $\Gamma_{D}$, i.e.,

$D=\\{i\,:\,\mathbf{x}_{i}\in\Gamma_{D}\\},$

and $I=J\setminus D$. We call $N_{i}$ the Lagrangian basis function associated
with the node $\mathbf{x}_{i}$. Observe that the space
${S}^{h}=\mathrm{span}\\{N_{i}\\}_{i\in J}$.

We denote the set of functions in ${S}^{h}$ which satisfy the Dirichlet
boundary condition by ${U}^{h}$ and the set of functions ${S}^{h}$ which
attain zero value on the Dirichlet boundary by ${V}^{h}$. Finally, we
introduce ${G}^{h}=\mathrm{span}\\{N_{i}\\}_{i\in D}$, i.e., the span of the
basis functions associated with the set of nodes in $D$. We observe that the
space ${V}^{h}$ is the span of the basis functions associated with each node
except the ones where Dirichlet conditions are imposed. Moreover, the
following relations hold:

$\begin{array}[]{lcll l}{U}^{h}&=&U\cap{S}^{h},\\\
{V}^{h}&=&V\cap{S}^{h}&=&\mathrm{span}\\{N_{A}\\}_{A\in I},\\\
{S}^{h}&=&{V}^{h}\bigoplus{G}^{h}&=&\mathrm{span}\\{N_{A}\\}_{A\in J},\\\
{U}^{h}&\subset&{S}^{h}.\end{array}$ (4)

###  2.2 The flow problem

For an incompressible single-phase fluid, the flow through a porous medium is
described by

$\begin{array}[]{c c c l}\nabla\cdot\underline{u}&=&0&\quad\text{in
}\Omega,\\\ \underline{u}&=&-k\nabla p&\quad\text{in }\Omega,\\\
p&=&g&\quad\text{on}\,\,\Gamma_{D},\\\ k\nabla
p\cdot\underline{n}&=&h&\quad\text{on}\,\,\Gamma_{N},\end{array}$ (5)

where $\underline{u}$ is the velocity and $p$ the pressure. For the flow
problem, the space $S$ coincides with the Sobolev space $H^{1}(\Omega)$. The
weak formulation of problem (5) reads:

$\begin{array}[]{l}\text{Find }p\in U\,\text{ such that }\\\\[5.69054pt]
{\displaystyle d(p,v)=f(v)\quad\forall v\in V,}\end{array}$ (6)

where

$d(p,v)=\int_{\Omega}k\,\nabla p\cdot\nabla q\,\text{d}V\quad\text{and}\quad
f(v)=\int_{\Gamma_{N}}h\,v\,\text{d}A.$

###### Theorem 1

The following a-priori estimates hold for problem (6):

$\displaystyle\min_{\partial\Omega}p<p<\max_{\partial\Omega}p$
$\displaystyle\qquad\quad(\text{Maximum Principle}),$ $\displaystyle
p|_{\partial\Omega}\geq 0\implies p\geq 0$
$\displaystyle\qquad\quad\text{(Positivity preservation)}.$

With the definitions introduced in Subsection 2.1, the FE approximation to
problem (5) reads

$\begin{array}[]{l}\text{Find }{p}^{h}\in{U}^{h}\,\text{ such that
}\\\\[5.69054pt] {\displaystyle
d({p}^{h},{v}^{h})=f({v}^{h})\quad\forall{v}^{h}\in{V}^{h},}\end{array}$ (7)

where

$d({p}^{h},{v}^{h})=\int_{\Omega}k\,\nabla p_{h}\cdot\nabla
v_{h}\,\text{d}V\quad\text{and}\quad
f(v)=\int_{\Gamma_{N}}h\,v_{h}\,\text{d}A.$

Hence, the flow problem admits the following algebraic representation

$\mathbf{D}\\\ \mathbf{p}\\\ =\mathbf{f},$ (8)

or, explicitly,

$\begin{pmatrix}\mathbf{D}_{II}&\mathbf{D}_{ID}\\\
\mathbf{0}&\mathbf{I}\end{pmatrix}\begin{pmatrix}\mathbf{p}_{A}\\\
\mathbf{p}_{D}\end{pmatrix}=\begin{pmatrix}\mathbf{f}_{I}\\\
\mathbf{g}_{D}\end{pmatrix},$ (9)

where

$\begin{array}[]{l c l c l}[\mathbf{D}_{II}]_{ij}&=&d(N_{j},N_{i})&&i,j\in
I,\\\ [\mathbf{D}_{ID}]_{ij}&=&d(N_{j},N_{i})&&i\in I,\,j\in D,\\\
[\mathbf{f}_{I}]_{i}&=&f(N_{i})&&i\in I.\\\ \end{array}$

###  2.3 The transport problem

The transport problem is coupled to the flow problem through the velocity
$\underline{u}$ as it enters both in the definition of the boundary conditions
and in the transport equation. We let
$\Gamma_{in}:=\\{\mathbf{x}\in\partial\Omega\,:\,\underline{u}\cdot\underline{n}<0\\}$
denote the inflow boundary and
$\Gamma_{out}:=\\{\mathbf{x}\in\partial\Omega\,:\,\underline{u}\cdot\underline{n}\geq
0\\}$ denote the outflow boundary. The unknown is the concentration
$c=c(\mathbf{x},t)$ of a given solute. As we consider only Dirichlet problems,
we have $\Gamma_{in}=\Gamma_{D}$. The equi-dimensional model for the advective
transport problem reads:

$\left\\{\begin{array}[]{r c l l}{\displaystyle\phi\frac{\partial c}{\partial
t}+\nabla\cdot(\underline{u}\,c)}&=&0&\quad\text{in}\,\,\Omega\times\mathscr{I}\\\
c(\,\cdot\,,0)&=&c_{0}&\quad\text{in}\,\,\Omega\\\
c&=&g&\quad\text{on}\,\,\Gamma_{in}\times\mathscr{I}\end{array}\right.$ (10)

Here, $\mathscr{I}=(0,T_{\text{fin}}]$ is the time interval, $c_{0}$ is the
initial condition, and $g$ is the inflow Dirichlet boundary condition.

Defining by $\Sigma$ the set of points where initial and boundary conditions
are prescribed, i.e.,

$\Sigma:=\\{(\mathbf{x},t)\,:\,\mathbf{x}\in\Gamma_{in}\lor t=0\\},$

the following theorem holds.

###### Theorem 2

For problem (10) the following a-priori estimates hold:

$\begin{array}[]{r c l
l}\nabla\cdot\underline{u}=0&\implies&\min_{\Sigma}c\leq
c\leq\max_{\Sigma}c&\quad(\text{Maximum Principle}),\\\\[2.84526pt]
c|_{\Sigma}\geq 0&\implies&c\geq 0&\quad\text{(Positivity
preservation)}.\end{array}$

Theorem 2 ensures that if $c_{0}$ and $g$ are positive, then the solution is
positive all over the space-time domain $\Omega\times\mathscr{I}$. Moreover,
the concentration in the domain cannot be larger than the maximum value
injected at the inflow boundary.

#### Finite element discretization of transport equation

The function $c_{h}=c_{h}(\,\cdot\,,t)\in{U}^{h}$ is the approximation of the
concentration $c$ at time $t$ and $\mathbf{c}(t)$ is the time dependent array
having as components the unknown coefficients $c_{j}$ with respect to the
basis $\\{N_{j}\\}$. The spatial discretization of (10) reads:

$\begin{array}[]{l}\text{For all }t\in\mathscr{I},\text{ find
}c_{h}(\,\cdot\,,t)\in U_{h}\,\text{ such that
}c(\,\cdot\,,0)=c_{0},\,\,\text{and}\\\\[5.69054pt] {\displaystyle
m\left(\dfrac{\partial c_{h}}{\partial
t},q_{h}\right)+a(c_{h},q_{h})=0\quad\forall q_{h}\in V_{h}.}\end{array}$ (11)

where we have used the following discrete bilinear forms:

$\begin{array}[]{ l c l}m\left(\dfrac{\partial c_{h}}{\partial
t},q_{h}\right)&=&{\displaystyle\hphantom{-}\int_{\Omega}\phi\,\dfrac{\partial
c_{h}}{\partial t}\,q_{h}\,\text{d}V,}\\\\[8.53581pt]
a(c_{h},q_{h})&=&{\displaystyle-\int_{\Omega}c_{h}\,\underline{u}\cdot\nabla
q_{h}\,\text{d}V+\int_{\Gamma_{out}}c_{h}\,q_{h}\,\underline{u}\cdot\underline{n}\,\text{d}A\,.}\end{array}$

Problem (11) admits the following algebraic representation

$\mathbf{M}\frac{\text{d}\mathbf{c}}{\text{d}t}+\mathbf{A}\mathbf{c}=\mathbf{0},$
(12)

where the matrices $\mathbf{M}$ and $\mathbf{A}$ are the scaled mass matrix
and the discrete advection operator, respectively, and will be specified
later.

A classic theorem that provides sufficient conditions to ensure the maximum
principle and positivity preservation for the semi-discrete transport problem
(12) is the following [58].

###### Theorem 3

Suppose that:

  * 1.

$[\mathbf{M}]_{ii}>0$,

  * 2.

$[\mathbf{M}]_{ij}=0$ for $i\neq j$,

  * 3.

$[\mathbf{A}]_{ij}\leq 0$ for $i\neq j$.

Then the following a priori estimates hold for the coefficient $c_{i}$:

  * 1a

The semi-Discrete Maximum Principle (DMP) is satisfied:

$\quad\quad{\displaystyle\sum_{j}[\mathbf{A}]_{ij}=0,\quad c_{j}\geq
c_{i},\quad\forall j\neq i\quad\implies\quad\dfrac{dc_{i}}{dt}\leq 0.}$

  * 1b

The positivity preservation is satisfied:

$c_{j}(0)\geq 0,\,\forall j\,\,\,\,\implies c_{i}(t)\geq 0,\forall t>0.$

Employing an implicit Euler scheme, one needs to solve a linear system at each
time step. The fully discrete counterpart of problem (12) is a sparse linear
system of the form

$\begin{pmatrix}\mathbf{B}_{II}&\mathbf{B}_{ID}\\\
\mathbf{0}&\mathbf{I}\end{pmatrix}\begin{pmatrix}\mathbf{c}_{I}^{n+1}\\\
\mathbf{c}_{D}^{n+1}\end{pmatrix}=\begin{pmatrix}\mathbf{M}_{II}&\mathbf{M}_{ID}\\\
\mathbf{0}&\mathbf{0}\end{pmatrix}\begin{pmatrix}\mathbf{c}_{I}^{n}\\\
\mathbf{c}_{D}^{n}\end{pmatrix}+\begin{pmatrix}\mathbf{0}_{I}\\\
\mathbf{g}_{D}\end{pmatrix},$ (13)

where

$[\mathbf{B}]_{ij}=[\mathbf{M}]_{ij}+\Delta t[\mathbf{A}]_{ij}.$

with $\Delta t$ being the time-step size. We point out that the backward Euler
scheme is unconditionally stable with no restriction on the time-step size
[59].

###  2.4 Discrete maximum principle

For FE approximations, a maximum principle for the discrete solutions, in
general, does not descend from the properties of the continuous counterparts
(5) and (10). In this section, we focus on the global _discrete maximum
principle_. For a detailed discussion on local definition, we refer to [51,
58].

###### Definition 1

The solutions to (9) or (13) satisfy the discrete maximum principle if

$\min_{j}g_{j}\leq u_{i}\leq\max_{j}g_{j}.$ (14)

###### Definition 2

The solutions to (9) or (13) are said to be globally positivity-preserving if

$\mathbf{g}\geq 0\implies\mathbf{u}\geq 0.$ (15)

Typical proofs of (14) and (15) are based on the theory of monotone matrices
and in particular on M-matrices.

###### Definition 3

A regular matrix $\mathbf{A}$ is said to be monotone if
$\mathbf{A}^{-1}\geq\mathbf{0}$ or, equivalently, if

$\mathbf{A}\mathbf{x}\geq\mathbf{0}\implies\mathbf{x}\geq\mathbf{0}.$ (16)

###### Definition 4

A monotone matrix $\mathbf{A}$ with $[\mathbf{A}]_{ij}\leq 0$ for $i\not=j$ is
called an M-matrix.

In general, the conditions expressed by Definitions 3 and 4 are quite
complicated to verify. The following theorem provides some mild conditions for
a matrix to be an M-matrix [60].

###### Theorem 4

A matrix $\mathbf{A}$ that

  1. 1.

is diagonally dominant by rows and strictly diagonally dominant for at least
one row,

  2. 2.

satisfies $[\mathbf{A}]_{ij}\leq 0$ for $i\not=j$,

  3. 3.

satisfies $[\mathbf{A}]_{ii}>0$,

is an M-matrix.

####  2.4.1 Discrete Diffusion Operators

The violation of the DMP conditions in Theorem 4 may be caused by the presence
of positive extra-diagonal entries or negative diagonal entries. Negative
diagonal entries can be found only in the discrete advection operator, while
positive extra-diagonal entries can be found both in the the discrete
advection operator and in the mass-matrix. A strategy employed to remove such
entries is based on adding a suitable artificial diffusion operator
$\mathbf{S}$ [58].

###### Definition 5

A symmetric matrix $\mathbf{S}$ is called a discrete diffusion operator if

  * 1.

$[\mathbf{S}]_{ij}=[\mathbf{S}]_{ji}\leq 0$ for $i\not=j$,

  * 2.

$\sum_{i}[\mathbf{S}]_{ij}=0$.

Observe that for a discrete diffusion operator the diagonal entries are
positive and the column sum is negative.

Given a matrix $\mathbf{Q}$, we call $\mathbf{S}^{Q}$ the discrete diffusion
operator for which:

  * 1.

$[\mathbf{S}^{Q}]_{ij}=-\,\max(0,[\mathbf{Q}]_{ij},[\mathbf{Q}]_{ji})\quad\text{for
}i\neq j$,

  * 2.

$[\mathbf{S}^{Q}]_{ii}=-\sum_{i\neq j}[\mathbf{S}^{Q}]_{ij}$.

The matrix $\tilde{\mathbf{Q}}=\mathbf{Q}+\mathbf{S}^{Q}$ is by construction
an M-matrix.

For a matrix $\mathbf{Q}$ arising from the FE discretization of a bilinear
form, the construction of $\mathbf{S}^{Q}$ can be performed during the
assembly of the matrix $\mathbf{Q}$.

###  2.5 Discrete maximum principle for flow problem

For FE discretizations, trying to satisfy the conditions of Theorem 4 leads,
in general, to geometric restrictions to the elements’ shapes. Actually, the
imposition of such conditions for any element-wise contribution to the
stiffness matrix can be translated into specific constraints on the angles
comprised between two different sides of the same element and on the ratio
between the shortest and the largest lengths (i.e., the aspect ratio). In
particular, for two-dimensional problems, only angles smaller then $\pi/2$ are
allowed in a mesh. For quadrilateral meshes with rectangular elements, the
elements have to be of a non-narrow type, i.e., with aspect ratio smaller than
$\sqrt{2}$.

If a mesh does not satisfy those geometrical constraints, one may replace the
operator ${\mathbf{D}}$ with a corresponding stabilized counterpart
$\tilde{\mathbf{D}}=\mathbf{D}+\mathbf{S}^{D}$. The addition of an algebraic
diffusion operator provides a stiffness matrix that satisfies the hypothesis
of Theorem 4.

###  2.6 Discrete maximum principle for the transport problem

Continuous FE discretizations applied to transport problems are characterized
by spurious undershoots and overshoots in the numerical solution as the
matrices arising do not satisfy the hypothesis of Theorem 4. Several
_stabilization_ strategies have been proposed to make $\mathbf{M}$ and
$\mathbf{A}$ consistent with the hypotheses of the theorem. In the next, we
refer to the algebraic flux correction (AFC) scheme proposed by [58]. The
basic idea of this method is to employ an algebraic diffusion operator that
allows ensuring the DMP and preserve of positivity of the solution.

The AFC is a predictor-corrector strategy which consists of two steps:

  1. 1.

advance the solution in time by a low-order scheme that incorporates enough
numerical diffusion to suppress undershoots and overshoots;

  2. 2.

add a correction to achieve higher accuracy without violating the DMP.

In order to satisfy the hypotheses of Theorem 3, one needs to $1)$ approximate
the consistent mass matrix $\mathbf{M}$ by its lumped counterpart
$\mathbf{M}_{L}$ and $2)$ eliminate the off-diagonal positive entries of the
transport operator $\mathbf{A}$ by adding a discrete diffusion operator, i.e.,
replacing $\mathbf{A}$ with $\tilde{\mathbf{A}}=\mathbf{A}+\mathbf{S}^{A}$.

In the process of AFC, the discrete problem (13) can be split into a good
diffusive part of the form

$\mathbf{M}_{L}\mathbf{c}^{n+1}+\Delta
t\,\tilde{\mathbf{A}}\,\mathbf{c}^{n+1}=\mathbf{M}_{L}\mathbf{c}^{n},$

and a bad anti-diffusive part given by

$\mathbf{f}=(\mathbf{M}_{L}-\mathbf{M})\mathbf{c}^{n+1}+\Delta
t\,\mathbf{S}^{A}\mathbf{c}^{n+1}-(\mathbf{M}_{L}-\mathbf{M})\mathbf{c}^{n}.$

The good diffusive counterpart is positivity-preserving but overly diffusive
due to its linearity. Hence, an anti-diffusive correction is required to
recover the accuracy of the original high-order discretization.

The bad anti-diffusive term admits the decomposition
$[\,\mathbf{f}\,]_{i}=\sum_{j\neq i}[\mathbf{F}]_{ij}$ where the components of
the sum represent a numerical flux which attain a local conservation property,
i.e. $[\mathbf{F}]_{ij}+[\mathbf{F}]_{ji}=0$ and are defined as follows:

$[\mathbf{F}]_{ij}=[\mathbf{M}]_{ij}(\mathbf{c}^{n+1}_{i}-\mathbf{c}_{j}^{n+1})-[\mathbf{M}]_{ij}(\mathbf{c}^{n}_{i}-\mathbf{c}_{j}^{n})-\Delta
t\,[\mathbf{S}]_{ij}^{D}(\mathbf{c}^{n+1}_{i}-\mathbf{c}_{j}^{n+1}).$

The anti-diffusive flux contributions $[\mathbf{F}]_{ij}$ are limited in such
a way that overall equations becomes less diffusive but spurious oscillations
are still suppressed. To this aim, the terms $[\mathbf{F}]_{ij}$ are
multiplied by a solution-dependent correction factor $0\leq\alpha\leq 1$ and
then added to the good diffusive part of the problem:

$\mathbf{M}_{L}\mathbf{c}^{n+1}+\Delta
t\,\tilde{\mathbf{A}}\,\mathbf{c}^{n+1}=\mathbf{M}_{L}\mathbf{c}^{n}+\tilde{\mathbf{f}},$

where $[\,\tilde{\mathbf{f}}\,]_{i}=\sum_{j\neq
i}[\mathbf{F}]_{ij}[\mathbf{\alpha}]_{ij}$ is the modified correction vector.
If all the factors $[\mathbf{\alpha}]_{ij}$ are equal to 1 the original high-
order Galerkin discretization is recovered, while the lower order
discretization is attainied for $[\mathbf{\alpha}]_{ij}=0$. A detailed
description of the technique used to compute the coefficients
$[\mathbf{\alpha}]_{ij}$ can be found in [58].

##  3 Discrete maximum principle
on non-conforming meshes

In this section, we present the derivation of FE discretizations which satisfy
the DMP on non-conforming meshes. As the DMP properties of the matrices and
the discrete diffusion operators depend on the basis employed, we show how
basis functions on non-conforming meshes are constructed from the shape
functions usually employed in the assembly of FE matrices. Finally, we present
a detailed study of the assembly of two matrices. For a diffusion problem,
such as the flow problem, we show that on non-conforming meshes with square
elements, the conditions for the DMP may not be satisfied and its
discretization may require stabilization with a discrete diffusion operator to
ensure the DMP. For a transport problem, we show that a discrete diffusion
operator constructed employing the elemental shape functions, which may be
intuitively chosen in the assembly process, and the one constructed employing
the FE basis functions provide different stabilization matrices and that only
the second one ensures the DMP.

###  3.1 Conforming spaces on non-conforming meshes

We now let the possibility of $\mathcal{T}$ be a non-conforming mesh, i.e., a
mesh where a vertex of an element may belong to the interior of an edge or a
face of another element. Such a vertex is said to be hanging. Non-hanging
vertices are instead said to be regular. We restrict ourselves to cases of
$1$-irregular meshes, i.e., the interior of an edge or a face can have at most
one hanging vertex. In particular, the hanging vertex corresponds to the
center of an edge or of a face of another element. Mesh of this kind are
usually obtained when AMR is employed.

As we work with low-order FE method and Lagrangian bases, nodes will be placed
at the vertices and, hence, they inherit the same terminology, i.e. regular
nodes and hanging nodes. We denote the set of all nodes by $J^{H}$ and its
cardinality by $n^{H}$. while we denote the set of regular nodes by $J^{r}$
and its cardinality by $n^{r}$. We define the set of non-conforming basis
functions as

$B^{H}=\\{N^{H}_{i}\,:\,\Omega\to\mathbb{R},\text{ such that
}N^{H}_{i}|_{E}\in\mathbb{L}_{1}\text{ and
}N^{H}_{i}|_{E}(\mathbf{x}_{j})=\delta_{ij}\\},$

and the set of conforming basis functions as

$B^{r}=\\{N^{r}_{i}\,\in C^{0}(\Omega),\text{ such that
}N^{r}_{i}|_{E}\in\mathbb{L}_{1}\text{ and
}N^{r}_{i}|_{E}(\mathbf{x}_{j})=\delta_{ij}\\},$

with $\delta_{ij}$ begin the Kronecker delta. A function $N^{H}_{i}$ is
created by glueing together the elemental shape functions which attains value
one at the node with coordinates $\mathbf{x}_{i}$. We respectively define the
non-conforming and conforming interpolation spaces as

$W^{H}=\mathrm{span}\,B^{H}\quad\text{and}\quad\quad
S^{h}=\mathrm{span}\,B^{r}.$

The space $W^{H}$ contains also some discontinuous functions as degrees of
freedom are associated with all the nodes of the mesh. Its dimension is
$n^{H}$. The space $S^{h}$ contains only continuous functions as degrees of
freedom are associated only with regular nodes of the mesh. Its dimension is
$n^{r}$. The function space $S^{h}$ is spanned by the Lagrangian basis
functions $B^{r}$, but its definition coincides with the one introduced in
(4).

The difference from the case of conforming meshes is that the support of
Lagrangian basis functions $B^{r}$ is not composed only of the elements that
are adjacent to the node. We explain this concept referring to the simple
example of the mesh depicted in Figure 1.

Figure 1: Exemplary of a non conforming mesh with
$J^{thisshouldnotbeused}=\\{13,14\\}$ and
$J^{r}=\\{1,2,3,4,5,6,7,8,9,10,11,12\\}$.

Regular nodes are numbered from $1$ to $12$, while hanging nodes have indices
$13$ and $14$.

(a) $N^{r}_{2}$ (b) $N^{r}_{5}$

Figure 2: Two basis functions of the conforming space $V^{r}$ defined over the
non-conforming mesh in Figure 1.

For the basis function $N^{r}_{2}$ the support is formed by elements $E_{6}$
and $E_{2}$, which are adjacent to $\mathbf{x}_{2}$ but also by element
$E_{8}$. For the basis function $N^{r}_{5}$, the support is formed by the
adjacent elements $E_{2}$, $E_{4}$, $E_{3}$, and $E_{8}$ and by the elements
$E_{6}$ and $E_{7}$.

The sets $B^{r}(E)$ and $B^{H}(E)$ denote the conforming and non-conforming
basis function on the element $E$. Observe that $B^{H}(E)$ is composed of the
elemental shape functions. Observe that $B^{H}(E)$ is composed of the shape
functions defined on $E$. On a single element, both these sets are a basis for
$\mathbb{L}_{1}$. We can introduce a change of basis from $B^{H}(E)$ to
$B^{r}(E)$. For example, for the element $E_{6}$, we have

$\begin{array}[]{l c l}N^{r}_{10}&=&N^{H}_{10},\\\
N^{r}_{2\phantom{0}}&=&N^{H}_{2}+\dfrac{1}{2}N^{H}_{13},\\\
N^{r}_{5\phantom{0}}&=&\dfrac{1}{2}N^{H}_{13},\\\
N^{r}_{12}&=&N^{H}_{12}.\end{array}$ (17)

The basis functions of $S^{h}$ are constructed by removing the basis functions
of $W^{H}$ at the hanging nodes and enriching the basis functions $N^{r}_{a}$
with contributions from the removed ones. Observe that $N^{r}_{2}$ and
$N^{r}_{5}$ attain value $1/2$ at hanging nodes.

#### Restriction and prolongation operators

We consider the evaluation of a linear functional $w$ on an element $E_{2}$,
with respect to $B^{r}(E_{2})$ and $B^{H}(E_{2})$, i.e.,

$\begin{array}[]{l c l l}w_{i}^{r}&=&w(N^{r}_{i}),&i\in
J^{r}_{E_{2}}\\\\[5.69054pt] w_{i}^{H}&=&w(N^{H}_{i}),&i\in
J^{H}_{E_{2}}.\end{array}$

Given the change of basis (17), we can define the elemental restriction
operator $\mathcal{R}_{E_{2}}$ as

$\displaystyle w^{r}_{10}$ $\displaystyle=\phantom{\dfrac{1}{2}}w^{H}_{10}$
$\displaystyle w^{r}_{2\phantom{0}}$
$\displaystyle=\phantom{\dfrac{1}{2}}w^{H}_{2}+\dfrac{1}{2}w^{H}_{13}$
$\displaystyle w^{r}_{5\phantom{0}}$ $\displaystyle=\dfrac{1}{2}w^{H}_{13}$
$\displaystyle w^{r}_{12}$ $\displaystyle=\phantom{\dfrac{1}{2}}w^{H}_{12}.$
(18)

With a little abuse of notation, i.e., denoting the algebraic representation
of the restriction operator as the restriction operator itself, equation (18)
admits the following algebraic representation
$\mathbf{w}^{r}_{E_{2}}=\mathcal{R}_{E_{2}}\mathbf{w}^{H}_{E_{2}}$:

$\begin{pmatrix}w^{r}_{10}\\\\[2.84526pt] w^{r}_{2\phantom{0}}\\\\[2.84526pt]
w^{r}_{5\phantom{0}}\\\\[2.84526pt]
w^{r}_{12}\end{pmatrix}=\begin{pmatrix}1&0&0&0\\\\[1.42262pt]
0&1&1/2&0\\\\[1.42262pt] 0&0&1/2&0\\\\[1.42262pt]
0&0&0&1\end{pmatrix}\begin{pmatrix}w^{H}_{10}\\\\[2.84526pt]
w^{H}_{2\phantom{0}}\\\\[2.84526pt] w^{H}_{13}\\\\[2.84526pt]
w^{H}_{12}\end{pmatrix}$

The elemental prolongation operator $\mathcal{P}_{E_{2}}$ is the adjoint of
$\mathcal{R}_{E_{2}}$. For elements that do not have hanging nodes, the
elemental prolongation and restriction operators are the identity. Putting
together all the elemental operators we can construct the global restriction
and prolongation operators, denoted by $\mathcal{R}$ and $\mathcal{P}$,
respectively. Their algebraic representations are still denoted with the same
symbols $\mathcal{R}\in\mathbb{R}^{n^{r}\times n^{H}}$ and
$\mathcal{P}\in\mathbb{R}^{n^{H}\times n^{r}}$.

###  3.2 DMP conditions on non-conforming meshes

We refer to an abstract variational problem:

$\begin{array}[]{l}\text{Find }p\in U\,\text{ such that }\\\\[5.69054pt]
{\displaystyle a(p,q)=f(q)\quad\forall{q}\in V},\end{array}$ (19)

for which we assume that a (continuous) maximum principle holds. The solution
of the problem above is equivalent to the solution of the linear system

$\mathbf{A}^{r}\mathbf{u}^{r}=\mathbf{f}^{r},$ (20)

where

$\begin{array}[]{ l c l
l}[\mathbf{A}^{r}]_{ij}&=&a(N^{r}_{j},N^{r}_{i}),&\quad\text{for}\,\,i,j\in
J^{r},\\\\[5.69054pt]
[\mathbf{f}^{r}]_{i}&=&f(N^{r}_{i}),&\quad\text{for}\,\,i\in
J^{r}.\end{array}$

In the process of assembling FE matrices, only elemental shape functions are
available. Hence, it is possibile to easily assemble

$\begin{array}[]{ l c l
l}[\mathbf{A}^{H}]_{ij}&=&a(N^{H}_{j},N^{H}_{i}),&\quad\text{for}\,\,i,j\in
J^{H},\\\\[5.69054pt]
[\mathbf{f}^{H}]_{i}&=&f(N^{H}_{i}),&\quad\text{for}\,\,i\in
J^{H}.\end{array}$

The stiffness matrix and right-hand-side associated with the conforming
discretization space can be obtained by applying the following transformations

$\mathbf{A}^{r}=\mathcal{R}\mathbf{A}^{H}\mathcal{P},\quad\mathbf{f}^{r}=\mathcal{R}\mathbf{f}^{H}.$
(21)

Once the solution $\mathbf{p}^{r}$ is found, it may be interpolated on the
non-conforming space $W^{H}$ applying

$\mathbf{p}^{H}=\mathcal{P}\mathbf{p}^{r}.$

In order to ensure the DMP for problem (19), the conditions of the Theorem 4
has to be verified for the matrix in (20). Even if the matrix $\mathbf{A}^{H}$
would fulfill such conditions, the operation (21) may introduce positive
extra-diagonal entries.

###  3.3 Finite element assembly on non-conforming meshes

The FE assembly of the problem (20) is realized on each element $E$
constructing the matrix $\mathbf{A}^{H}_{E}$ and the vector
$\mathbf{f}^{H}_{E}$. The restriction operations (21) are hence realized per
element, resulting in

$\begin{array}[]{l c
c}\mathbf{A}_{E}^{r}&=&\mathcal{R}_{E}\mathbf{A}_{E}^{H}\mathcal{P}_{E},\\\\[5.69054pt]
\mathbf{f}_{E}^{r}&=&\mathcal{R}_{E}\mathbf{f}_{E}^{H}.\end{array}$

The operation above, may introduce positive extra-diagonal entries in
$\mathbf{A}_{E}^{r}$, even if the matrix $\mathbf{A}_{E}^{H}$ has none. Hence,
one needs to check if the assumptions of Theorem 4 are satisfied for the
matrices $\mathbf{A}_{E}^{r}$ and in case add the elemental contribution of
the algebraic diffusion operator. In Figure 1, we report the steps to assemble
a stiffness matrix that ensures the DMP. In order to ensure the DMP, we check
that all local stiffness matrices have positive diagonal entries, and negative
extra-diagonal entries. In this way, the stiffness matrix, before setting the
Dirichlet boundary conditions, will have the same properties, and, finally,
the imposition of Dirichlet boundary conditions makes $\mathbf{A}^{r}$ an
M-matrix.

Data: Bilinear form $a(\,\cdot\,,\,\cdot\,)$, linear functional
$F(\,\cdot\,)$, mesh $\mathcal{T}$

Result: The associated FE matrix $\mathbf{A}^{r}$ and right-hand-side
$\mathbf{f}^{r}$, with $\mathbf{A}^{r}$ which satisfies the DMP condition.

for  _$E\leftarrow\mathcal{T}$_ do

Assemble $\mathbf{A}^{H}_{E}$ and $\mathbf{f}^{H}_{E}$;

$\mathbf{A}_{E}^{r}=\mathcal{R}_{E}\mathbf{A}_{E}^{H}\mathcal{P}_{E}$ and
$\mathbf{f}_{E}^{r}=\mathcal{R}_{E}\mathbf{f}_{E}^{H}$;

if  _$\mathbf{A}_{E}^{r}$ has positive extra-diagonal entries_ then

Compute $\mathbf{S}_{E}^{A^{r}_{E}}$;

$\mathbf{A}_{E}^{r}\leftarrow\mathbf{A}_{E}^{r}+\mathbf{S}_{E}^{A^{r}_{E}}$;

end if

Modify rows of $\mathbf{A}^{r}$ and $\mathbf{f}^{r}$ to impose Dirichlet
boundary conditions;

end for

Algorithm 1 Algorithm to assemble problem (19) on a non-conforming mesh and
ensuring the DMP.

####  3.3.1 Failure of Discrete Maximum Principle conditions

In order to justify the use of the Algorithm 1, we present the assembly of two
matrices associated with a flow problem and a transport operator. We consider
a square element, whose nodes are numbered as follow $\mathbf{x}_{1}=(0,0)$,
$\mathbf{x}_{2}=(1,0)$, $\mathbf{x}_{3}=(0,1)$, $\mathbf{x}_{4}=(1,1)$. The
node $\mathbf{x}_{1}=(0,0)$ is a hanging node and the two regular nodes
connected to it are $\mathbf{x}_{0}=(-1,0)$ and $\mathbf{x}_{2}$. The
restriction operator that goes from a space where the index set is
$J^{H}_{E}=(1,2,3,4)$ to a space where the index set is $J^{r}_{E}=(0,2,3,4)$
is

$\mathcal{R}_{E}=\begin{pmatrix}1/2&0&0&0\\\\[1.42262pt]
1/2&1&0&0\\\\[1.42262pt] 0&0&1&0\\\\[1.42262pt] 0&0&0&1\end{pmatrix}.$

##### Diffusion operator

We consider the weak form associated

$\int_{E}\nabla N^{H}_{i}\cdot\nabla N^{H}_{i}\mbox{d}\Omega,$

and the associated local stiffness matrix

$\mathbf{A}_{E}^{H}=\begin{pmatrix}\phantom{-}2/3&-1/6&-1/6&-1/3\\\\[1.42262pt]
-1/6&\phantom{-}2/3&-1/3&-1/6\\\\[1.42262pt]
-1/6&-1/3&\phantom{-}2/3&-1/6\\\\[1.42262pt]
-1/3&-1/6&-1/6&\phantom{-}2/3\end{pmatrix}.$ (22)

The matrix $\mathbf{A}_{E}^{H}$ has positive diagonal entries and negative
extra-diagonal entries, as it usually happens for diffusion operators. After
applying the transformation (22), we obtain a local stiffness matrix evaluated
with respect to the regular basis that reads

$\mathbf{A}_{E}^{r}=\begin{pmatrix}\phantom{-}2/3\phantom{0}&\phantom{-}1/12&-1/3&-5/12\\\\[1.42262pt]
\phantom{-}1/12&\phantom{-}1/6\phantom{0}&-1/6&-1/12\\\\[1.42262pt]
-1/3\phantom{0}&-1/6\phantom{0}&\phantom{-}2/3&-1/6\phantom{0}\\\\[1.42262pt]
-5/12&-1/12&-1/6&\phantom{-}2/3\phantom{0}\end{pmatrix}.$

Hence, in order to have a local stiffness matrix that ensures the DMP
conditions, the following algebraic diffusion operator has be added to
$\mathbf{A}_{E}^{r}$.

$\mathbf{S}_{E}^{A_{E}^{r}}=\begin{pmatrix}\phantom{-}1/12&-1/12&0&0\\\\[1.42262pt]
-1/12&\phantom{-}1/12&0&0\\\\[1.42262pt] 0&0&0&0\\\\[1.42262pt]
0&0&0&0\end{pmatrix}.$

##### Transport operator

We consider now the weak form associated to the transport operator

$\int_{E}N^{H}_{j}\dfrac{\partial N^{H}_{i}}{\partial x_{1}}$

which corresponds to having chosen a velocity vector of the form
$\underline{u}=(1,0)$. The elemental stiffness matrix assembled on the non-
conforming basis functions and the associated algebraic diffusion operator are

$\mathbf{A}_{E}^{H}=\frac{1}{12}\begin{pmatrix}-2&-2&-1&-1\\\\[1.42262pt]
\phantom{-}2&\phantom{-}2&\phantom{-}1&\phantom{-}1\\\\[1.42262pt]
-1&-1&-2&-2\\\\[1.42262pt]
\phantom{-}1&\phantom{-}1&\phantom{-}2&\phantom{-}2\end{pmatrix},\quad\mathbf{S}_{E}^{A_{E}^{H}}=\frac{1}{12}\begin{pmatrix}\phantom{-}3&-2&\phantom{-}0&-1\\\\[1.42262pt]
-2&\phantom{-}4&-1&-1\\\\[1.42262pt]
\phantom{-}0&-1&\phantom{-}3&-2\\\\[1.42262pt]
-1&-1&-2&\phantom{-}4\end{pmatrix}.$

The matrix $\mathbf{A}_{E}^{H}$ has negative diagonal entries and positive
extra-diagonal entries. Moreover, it is not diagonally dominant. On the other
hand, the matrix $\mathbf{A}_{E}^{H}+\mathbf{S}_{E}^{H}$ respects the
conditions of Theorem 4, but the computation of
$\mathcal{R}_{E}(\mathbf{A}_{E}^{H}+\mathbf{S}_{E}^{H})\mathcal{P}_{E}$ gives
an elemental stiffness matrix which does not respect such conditions. On the
other hand, the computation of
$\mathbf{A}_{E}^{r}=\mathcal{R}_{E}\mathbf{A}_{E}^{H}\mathcal{P}_{E}$ and the
corresponding algebraic diffusion operator are

$\mathbf{A}_{E}^{r}=\frac{1}{24}\begin{pmatrix}-1&-3&-1&-1\\\\[1.42262pt]
\phantom{-}1&\phantom{-}3&\phantom{-}1&\phantom{-}1\\\\[1.42262pt]
-1&-3&-4&-4\\\\[1.42262pt]
\phantom{-}1&\phantom{-}3&\phantom{-}4&\phantom{-}4\end{pmatrix},\quad\mathbf{S}_{E}^{A_{E}^{r}}=\frac{1}{24}\begin{pmatrix}\phantom{-}2&-1&\phantom{-}0&-1\\\\[1.42262pt]
-1&\phantom{-}5&-1&-3\\\\[1.42262pt]
\phantom{-}0&-1&\phantom{-}5&-4\\\\[1.42262pt]
-1&-3&-4&\phantom{-}8\end{pmatrix}.$

The sum of $\mathbf{A}_{E}^{r}$ and $\mathbf{S}_{E}^{A_{E}^{r}}$ gives an
M-matrix.

##  4 Conservation properties of continuous finite element methods

In this section, we derive the global and local conservation properties of FE
discretizations on non-conforming meshes. We refer to a diffusion problem,
such as the flow problem, but they can be equally proven for all the cases
considered in [45]. The idea at the base of the proofs is that the basis
functions form a partition of unity, i.e., the sum of the evaluation of all
the basis functions at a point $\mathbf{x}$ is one. Equivalently, we can
observe that the function one belongs to the space of test functions or that
the sum of the rows of the stiffness matrix is zero.

While the proof of the global conservation properties is identical to the one
for conforming meshes, the main difference in the proof of the local
conservation properties is in the definition of a suitable space of test
functions. The use of this definition allows us to show that the FE method is
conservative across any interface that separates two subsets $\Omega$, without
the requirement that such an interface coincides with the sides of some
elements of the mesh.

##### Global conservation properties

The function spaces defined in equation (4) are spanned by the functions
$N^{r}_{i}$. In order to show the global conservation, we introduce the flux
$q_{D}(\Omega)$ along the Dirichlet boundary $\Gamma_{D}$, and restate problem
(7) as follows:

$\begin{array}[]{l}\text{Find}\,\,({u}^{h},q_{D}(\Omega))\in{U}^{h}\times{G}^{h}\,\,\text{such
that}\\\\[5.69054pt] {\displaystyle
d({p}^{h},{v}^{h})+(q_{D}(\Omega),{v}^{h})_{D}=f({v}^{h})\quad\forall{v}^{h}\in{S}^{h}}.\end{array}$
(23)

Setting ${v}^{h}=1$ into the above equation, we obtain

$\int_{\Gamma_{D}}q_{D}(\Omega)\,\text{d}{\Gamma}+\int_{\Gamma_{N}}h\,\text{d}\Gamma=0.$

This shows that that $q_{D}(\Omega)$ is the flux associated to the Dirichlet
boundary conditions that makes the method globally conservative.

Formulation (23) can be split into two sub-problems:
$\text{Find}\,\,({u}^{h},q_{D}(\Omega))\in{U}^{h}\times{G}^{h}\,\,\text{such
that}$

$\displaystyle d({p}^{h},{v}^{h})-f({v}^{h})$ $\displaystyle=0$
$\displaystyle\forall{v}^{h}\in{V}^{h},$ (24)
$\displaystyle(q_{D}(\Omega),{v}^{h})_{D}$
$\displaystyle=d({p}^{h},{v}^{h})-f({v}^{h})$
$\displaystyle\forall{v}^{h}\in{G}^{h},$ (25)

which admit the following algebraic representation:

$\begin{pmatrix}\mathbf{D}_{II}&0\\\
-\mathbf{D}_{ID}&\mathbf{B}_{DD}\end{pmatrix}\begin{pmatrix}\mathbf{p}_{I}\\\
\mathbf{h}_{D}\end{pmatrix}=\begin{pmatrix}\mathbf{f}_{I}-\mathbf{D}_{ID}\mathbf{g}_{D}\\\
-\mathbf{f}_{D}\end{pmatrix},$ (26)

In the linear system above, we have

$\begin{array}[]{l c l c
l}[\mathbf{B}_{DD}]_{ij}&=&(N^{r}_{j},N^{r}_{i})_{\Gamma_{D}}&&i,j\in D,\\\
[\mathbf{D}_{II}]_{ij}&=&d(N^{r}_{j},N^{r}_{i})&&i,j\in I,\\\
[\mathbf{D}_{ID}]_{ij}&=&d(N^{r}_{j},N^{r}_{i})&&i\in I,\,j\in D,\\\
[\mathbf{f}_{I}]_{i}&=&f(N^{r}_{i})&&i\in I.\\\ \end{array}$

The matrix $\mathbf{B}_{DD}$ is the boundary mass on $\Gamma_{D}$. From (25)
or, equivalently, from the second row of (26), we notice that the computation
of the boundary flux associated to $\Gamma_{D}$ requires the inversion of the
boundary mass matrix $\mathbf{B}_{DD}$, which can be thought as a post-
processing evaluation.

##### Local conservation properties

Let $\Omega_{1}\subset\Omega$ and $\Omega_{2}\subset\Omega$ be two non-
overlapping subdomains such that
$\overline{\Omega}_{1}\cup\overline{\Omega}_{2}=\overline{\Omega}$. We denote
the boundary of the subdomain $\Omega_{1}$ by $\Gamma_{1}$, the boundary of
the subdomain $\Omega_{2}$ by $\Gamma_{2}$, and the interface between the
$\Omega_{1}$ and $\Omega_{2}$ (i.e., $\Gamma_{1}\cap\Gamma_{2}$) by $L$ . The
subdomain $\Omega_{1}$ could be interior to $\Omega$ or intersect the boundary
of $\Omega$ in either a finite number of isolated points or a set with a
finite measure.

For the sake of simplicity, we assume that $\Gamma_{1}$ only overlaps with the
Neumann boundary $\Gamma_{N}$ and denote by
$\Gamma_{1}^{N}=\Gamma_{1}\cap\Gamma_{N}$ their intersection such that
$\Gamma_{1}=L\cup\Gamma_{1}^{N}$. In the same way, the boundary
$\Gamma_{2}\setminus L$ is split into the Neumann and the Dirichlet portions
and defined as $\Gamma_{2}=L\cup\Gamma_{2}^{N}\cup\Gamma_{D}$.

We refer to $J^{r}_{L}$ as the index set of the conforming basis functions
whose support is non-null on the interface ${L}$, i.e.,

$J^{r}_{L}=\\{i\,:\,|\text{supp}(N^{r}_{i})\cap L|>0\\}.$

and introduce the function space
$G_{L}^{h}=\mathrm{span}\\{N^{r}_{i}\\}_{\\{i\in J^{r}_{L}\\}}$. We also
define the bilinear form and the linear functional restricted to the subdomain
$\Omega_{a}$ for $a=\\{1,2\\}$ as follows:

$d_{a}({p},{q})=\int_{\Omega_{a}}k\,\nabla p\cdot\nabla
q\,\text{d}\Omega\quad\text{and}\quad
f_{a}(q)=\int_{\Gamma_{a}^{N}}h\,{q}\,\text{d}\Gamma,$

with $a=\\{1,2\\}$.

In order to derive the local conservation properties, we denote the
approximated flux across the boundary $L$ with respect to $\Omega_{1}$ by
$q_{L}(\Omega_{1})$.

The boundary flux $q_{L}(\Omega_{1})$ is the solution of following
formulation:

$\begin{array}[]{l}\text{Find}\,({p}^{h},q_{L}(\Omega_{1}))\in U^{h}\times
G_{L}^{h}\,\text{ such that }\\\\[5.69054pt]
({v}^{h},q_{L}(\Omega_{1}))_{L}=d_{1}({p}^{h},{v}^{h})-f_{1}({v}^{h})\quad\forall{v}^{h}\in
S^{h},\end{array}$ (27)

which can be split into two sub-problems:

$\displaystyle\text{Find}\,({p}^{h},q_{L}(\Omega_{1}))\in U^{h}\times
G_{L}^{h}\,\text{ such that }$ $\displaystyle
d_{1}({p}^{h},{v}^{h})-f_{1}({v}^{h})$ $\displaystyle=0$
$\displaystyle\forall{v}^{h}\in S^{h}-G_{L}^{h},$
$\displaystyle({v}^{h},q_{L}(\Omega_{1}))_{L}$
$\displaystyle=d_{1}({p}^{h},{v}^{h})-f_{1}({v}^{h})$
$\displaystyle\forall{v}^{h}\in G_{L}^{h}.\quad\quad$

In the same way, one can define the flux for the complementary domain
$q_{L}(\Omega_{2})$ across the boundary $L$ as the solution of the problem:

$\begin{array}[]{l}\text{Find}\,({p}^{h},q_{L}(\Omega_{2}))\in U^{h}\times
G_{L}^{h}\,\text{ such that }\\\\[5.69054pt]
({v}^{h},q_{L}(\Omega_{2}))_{L}=d_{2}({p}^{h},{v}^{h})-f_{2}({v}^{h})-({v}^{h},q_{D}(\Omega))_{\Gamma_{D}}\quad\forall{v}^{h}\in
S^{h},\end{array}$ (28)

We establish the local conservation by first proving that the interface
fluxes, $q_{L}(\Omega_{1})$ and $q_{L}(\Omega_{2})$ are conservative with
respect to corresponding subdomains $\Omega_{1}$ and $\Omega_{2}$. To this aim
we select any ${v}^{h}\in S^{h}$ such that ${v}^{h}_{\Omega_{a}}=1$ with
$a\in\\{1,2\\}$, and obtain:

$\int_{L}q_{L}(\Omega_{1})\,\text{d}\Gamma+\int_{\Gamma_{1}^{N}}h\,\text{d}\Gamma=0$

and

$\int_{L}q_{L}(\Omega_{2})\,\text{d}\Gamma+\int_{\Gamma_{2}^{N}}h\,\text{d}\Gamma+\int_{\Gamma_{D}}q_{D}^{h}(\Omega)\,\text{d}\Gamma=0.$

It only remains to prove the uniqueness of the boundary flux over $L$. Summing
together equation (27) and equation (28) and restricting ${v}^{h}$ to
$G_{L}^{h}\subset S^{h}$ leads to:

$({v}^{h},q_{L}(\Omega_{1})+q_{L}(\Omega_{2}))_{L}=0$

which states that the two fluxes equilibrate on $L$ and completes the proof.

Remark 1 If the vertices of the interface $L$ correspond to the vertices the
mesh $\mathcal{T}$, the local conservation is attained point-wise:

$q_{L}(\Omega_{1})+q_{L}(\Omega_{2})=0,$

and the coefficients of the boundary fluxes can be obtained by computing the
inverse of a boundary mass matrix restricted to $L$.

Remark 2 As the sum of the rows of a discrete diffusion operator is zero, due
to the properties listed in Definition 5, the global and local conservation
properties hold also for the stabilized formulations.

##  5 Numerical examples

In this section, we analyze the DMP and conservation properties of FE
discretizations on non-conforming meshes. The applications will be the coupled
flow and transport problems in fractured porous media based on an equi-
dimensional representation of fractures. For the flow problem, we focus both
on the conservation properties of the computed fluxes and the positivity of
the solution. For the transport problem, the main focus is to show the
positivity of the solution. Before discussing in detail the numerical results,
we briefly describe the AMR strategy employed to generate the non-conforming
adapted meshes for domains with fractures.

We consider three numerical examples. The first example is a regular fracture
network [56]. Being the fractures of this example axis-aligned, we can employ
a structured mesh that resolves the interfaces between the fractures and the
embedding matrix. The second example is a two-dimensional adaptation of the
benchmark proposed in [61]. The third example describes a realistic fracture
network [56].

We employ the first and the second benchmarks to estimate the conservative
boundary fluxes, and the third test case to show the key role that
stabilization plays to obtain solutions that satisfy the DMP. In order to
evaluate the accuracy of the boundary fluxes $q_{L}(\Omega_{i})$ with
$i\in\\{1,2\\}$, we compute the relative $L^{2}$ error:

$e_{q}=\dfrac{\|q_{L}(\Omega_{i})-q_{R}(\Omega_{j})\|_{L}}{\|q_{R}(\Omega_{i})\|_{L}},$
(29)

where $q_{R}(\Omega_{i})$ denotes the reference flux, usually obtained on a
very fine mesh. Moreover, we evaluate the total flux $Q_{L}(\Omega_{j})$ along
the boundary $L$ with respect to the sub-domains $\Omega_{1}$ and $\Omega_{2}$
as follows:

$Q_{L}(\Omega_{i})=\int_{L}q_{L}(\Omega_{i})d\Gamma.$ (30)

All numerical tests have been implemented in “Parrot”, an application
implemented in the FE framework MOOSE. Parrot allows for the simulation of
several problems in heterogeneous materials, such as flow and transport,
Biot’s equations in the time and time-frequency domain. The solution of the
linear system arising from the FE discretization is performed using the
parallel direct solver MUMPS.

###  5.1 Adaptive mesh refinement for
heterogenous media

In this work, we employ the strategy proposed in [54] to create adapted meshes
for fracture networks. The idea of the method is to start from an initial
uniform mesh. Then, given the distribution of embedded fractures
$\mathcal{F}$, one may iterate over all elements and mark for refinement the
ones that have non-empty intersections with at least one fracture. The uniform
refinement of such elements is a step of the AMR procedure. At each step, we
also mark for refinements the elements close to the refined one to keep the
meshes $1$-irregular. Iterating the refinement step several times allows
creating sequences of meshes, which become progressively finer close to the
interface between the fractures and the matrix, where the discretization error
is usually larger. As interfaces are not explicitly resolved, material
properties may be non-constant over the mesh elements close to the interfaces.
Hence, in the process of assembling the stiffness matrix, the material
properties are evaluated at the quadrature points.

We point out that the employed mesh generation is a flexible procedure that
allows us to obtain meshes for any distribution of the fractures and with some
desired features. For example, one may refine only elements that have a non-
empty overlap with the boundary of fractures or a region around the fracture
which does not necessarily coincide with its aperture.

###  5.2 Resolved and unresolved meshes

Meshes generated with the AMR strategy do not resolve the interfaces between
matrix and fractures. Hence, we will refer to them as _unresolved_ , or
_adapted_ , meshes. This kind of meshes is in contrast to the _resolved_
meshes, i.e., meshes that explicitly resolve the interfaces. For resolved
meshes, elements close to the interface have edges that lie on the interfaces.

Discretization methods based on resolved meshes are usually referred to as
fitted methods, while methods based on unresolved methods are referred to as
unfitted methods. Despite the lack of regularity of interface problems, fitted
methods, in general, provide optimal convergence rates compared to the
polynomial order. On the other hand, standard unfitted FE methods do not
provide optimal convergence rates but allow for larger flexibility in the
meshing of complicated domains and interfaces. For the elements on which
material properties are not constant, the standard convergence proofs used for
FE are not applicable.

We refer to the resolved meshes as $\mathcal{R}^{\text{be},\text{fe}}$, where
be is the number of vertical subdivisions in the mesh and fe is the number of
element crossing the fracture. We point out that large values of fe with
smaller values of be give meshes with elongated elements in the fracture,
which may affect the accuracy of the solution.

Unresolved meshes are denoted by $\mathcal{U}^{\text{be}}_{\text{amr}}$, where
amr is the number of AMR steps applied to the uniform background mesh
consisting of be initial subdivisions.

For each example, we report the mesh characteristics, together with the number
of elements $n^{E}$, and the number of elements per fracture aperture $w$. and
for the first two examples, for which a resolved mesh can be created, we
perform a comparison of the solutions between the resolved and the unresolved
meshes.

###  5.3 Regular fracture network

The geometry of this test consists in a regular fracture network embedded in
square domain with $L_{1}=L_{2}=1\,[\text{m}]$. Material properties are listed
in Table 1. A flow problem on this geometry has been originally presented in
[62]. A comparison between different solution methods for both hybrid- and
equi-dimensional representations of the fracture have been presented in [56].
A hybrid-dimensional representation for the transport has been proposed in
[63].

We consider several resolutions for the resolved and the unresolved meshes.
For the resolved meshes, we choose
$\text{be}\in\\{80,160,320,640,1280,2560,4000,5120\\}$ and
$\text{fe}\in\\{1,2,4,8,16\\}$. In order to study the accuracy of the boundary
fluxes, we use a reference solution computed on a resolved mesh with $5120$
elements in the background and $64$ elements in an area that is the double of
the fracture aperture. This means that $32$ elements are inside the fracture
and $16$ elements per side are close to each fracture side. We remind that the
total number of elements is $(\text{be}+3\text{fe})^{2}$, while the total
number of nodes is $(\text{be}+3\text{fe}+1)^{2}$. The characteristics of the
unresolved meshes are reported in Table 12. Pressure and concentration
profiles are analized along the segments $\text{A}\text{A}^{\prime}$,
$(x_{1}=0.5$ [m]$)$, $\text{B}\text{B}^{\prime}$, $(x_{2}=0.7$ [m]$)$, and
$\text{C}\text{C}^{\prime}$, $(x_{2}=0.5$ [m]$)$, and
$\text{D}\text{D}^{\prime}$, $(x_{2}=0.75$ [m]$)$ as reported in Figure 3.

####  5.3.1 Flow problem and local conservative flux calculations

For the flow problem, $\Gamma_{D}$ is the right side at $x_{1}=L_{1}$, while
$\Gamma_{N}$ is composed of the other three sides. We set $g=1$ [m] and $h=-1$
[m$\times$s-1] on the left side of $\Gamma_{N}$. On the top and bottom sides,
homogeneous Neumann conditions are imposed. Figure 3 shows the spatial
distribution of the pressure computed on the mesh $\mathcal{U}^{320}_{8}$.

Table 1: Material properties employed for the benchmark “Regular fracture network”. Property | Symbol | Value | Unit
---|---|---|---
Fracture thickness | $\delta$ | $1\times 10^{-4}$ | m
Matrix permeability | $k_{m}$ | $1\times 10^{4}$ | m$\times$s-1
Fracture permeability | $k_{f}$ | $1$ | m$\times$s-1
Matrix porosity | $\phi_{m}$ | $1$ |
Fracture porosity | $\phi_{f}$ | $1$ |
Figure 3: Pressure distribution for the test case $\mathcal{U}^{320}_{8}$. We
also report the extrema of the three segments along which the relevant
properties are computed. Figure 4: Comparison of pressure values along two
lines for the resolved test cases. Figure 5: Comparison of pressure values
along two lines for the unresolved test cases.

We compare the pressure on some resolved and unresolved meshes computed along
with the segments $\text{A}\text{A}^{\prime}$ and $\text{B}\text{B}^{\prime}$
and depict all the results in Figures 4 and 5, respectively. We observe that
for resolved meshes no visible differences are present along both the lines.
On the other hand, the accuracy of the solutions on unresolved meshes is
mostly affected by $w$, the number of elements per fracture width. For this
case, solutions obtained with meshes with $w>8$ present identical pressure
distributions. In Figure 6, we report a comparison between the solutions
computed on the unresolved and the resolved meshes. We observe that the two
solutions coincide with one of the solutions presented in [56], which has been
computed with a mimetic finite difference (MFD) method on a very fine mesh,
consisting of $1\,175\,056$ elements and $10$ elements per each fracture in
its normal direction.

Figure 6: Comparison of pressure values along two lines between resolved and
unresolved test cases. The MFD solution computed in [56] is reported in black.

To show the conservation properties of the FE on adapted meshes, we consider
two subdivisions of $\Omega$. The first one is the previously-introduced
$\Omega_{m}$ and $\Omega_{f}$, for which we compute the total flux along with
the interface $\Gamma$. As a second one, we consider a subdivision in
$\Omega_{1}$ and $\Omega_{2}$, which are divided by the segment
$L:=\text{B}\text{B}^{\prime}$. In Table 2, we report the value of the total
fluxes $Q_{L}(S)$ with $S\in\\{\Omega_{1},\Omega_{2},\Omega_{m},\Omega_{f}\\}$
obtained by using Equation (30) for some of the resolved and unresolved
meshes. As already proved in Section 4, all the numerical results confirm that
the FE method on adapted meshes attains local conservation properties. Indeed,
we see that the value of the flux computed on one side of an interface is
matched by the negative value on the other side of the domain. The total
fluxes computed across the interface $\Gamma$ satisfy the same property.

Table 2: Total fluxes, Equation 30, for the resolved and unresolved test cases. Mesh | $Q_{\Gamma}(\Omega_{m})$ | $Q_{\Gamma}(\Omega_{f})$ | $Q_{\text{B}\text{B}^{\prime}}(\Omega_{1})$ | $Q_{\text{B}\text{B}^{\prime}}(\Omega_{2})$
---|---|---|---|---
$\mathcal{R}^{320,4}$ | 0.66219 | -0.66219 | -0.11775 | 0.11775
$\mathcal{R}^{640,4}$ | 0.66219 | -0.66219 | -0.11775 | 0.11775
$\mathcal{R}^{1280,4}$ | 0.66218 | -0.66218 | -0.11775 | 0.11775
$\mathcal{R}^{4000,2}$ | 0.66218 | -0.66218 | -0.11775 | 0.11775
$\mathcal{U}^{80}_{9}$ | 0.65692 | -0.65692 | -0.11704 | 0.11704
$\mathcal{U}^{160}_{8}$ | 0.656915 | -0.656915 | -0.11691 | 0.11691
$\mathcal{U}^{320}_{8}$ | 0.662382 | -0.662382 | -0.11803 | 0.11803
$\mathcal{U}^{640}_{7}$ | 0.662381 | -0.662381 | -0.11799 | 0.11799
$\mathcal{U}^{1280}_{7}$ | 0.662415 | -0.662415 | -0.11794 | 0.11794

(a) Boundary flux across the region of the boundary
$\text{B}\text{B}^{\prime}$ intersecting fractures. (b) Boundary flux across
the portion of the boundary $\text{B}\text{B}^{\prime}$ lying in the embedding
matrix.

Figure 7: Boundary fluxes across boundary $\text{B}\text{B}^{\prime}$. We
report $q_{\text{B}\text{B}^{\prime}}(\Omega_{2})$ for the unresolved mesh,
and $q_{\text{B}\text{B}^{\prime}}(\Omega_{1})$ the resolved test cases. The
black line is the reference solution.

In Figure 7, we compare the fluxes computed on the unresolved mesh
$\mathcal{U}^{80}_{11}$ and the resolved mesh $\mathcal{R}^{1280,16}$ against
the fluxes computed for the reference solution. Such mesh resolutions have
been computed as they are the coarsest ones for each case that can correctly
reproduce the flux. For reproducing the flux correctly, unresolved meshes with
a coarse background are sufficient but they need a large number of AMR steps,
requiring $w>16$. On the other hand, fluxes computed on the resolved mesh
require a much finer background. Such results suggest that resolved meshes are
not enough to correctly reproduce the flux, but they also require a certain
number of elements close to the interfaces. This fact is also highlighted by
the errors reported in Table 3 and Table 4. We observe that the error is
mostly dominated by the value of $w$ for unresolved discretizations. For
resolved meshes, although the error decreases more regularly, it remains
larger even if a larger number of background and fracture elements are
employed.

Table 3: $L^{2}$ relative error, Equation 29, for the boundary flux $q_{\text{B}\text{B}^{\prime}}(\Omega_{1})$ computed on the resolved meshes. fe: be: | 80 | 160 | 320 | 640 | 1280 | 2560 | 5120
---|---|---|---|---|---|---|---
1 | 0.99439 | 0.98905 | 0.97863 | 0.95874 | 0.92264 | 0.86367 | 0.78035
2 | 0.56599 | 0.56301 | 0.56064 | 0.55603 | 0.54743 | 0.53292 | 0.51144
4 | 0.38873 | 0.38697 | 0.38609 | 0.38437 | 0.3811 | 0.37542 | 0.3667
8 | 0.25898 | 0.2578 | 0.25745 | 0.25679 | 0.25551 | 0.25325 | 0.24969
16 | 0.16424 | 0.16339 | 0.16323 | 0.16294 | 0.16239 | 0.16138 | 0.15975
amr: be: | 80 | 160 | 320 | 640 | 1280 | 2560
---|---|---|---|---|---|---
2 |  |  |  |  |  | 0.50383
3 |  |  |  |  | 0.50383 | 0.3782
4 |  |  |  | 0.50383 | 0.3782 | 0.24934
5 |  |  | 0.50383 | 0.3782 | 0.24934 | 0.1466
6 |  | 0.50383 | 0.3782 | 0.24934 | 0.1466 | 0.061532
7 | 0.50384 | 0.3782 | 0.24934 | 0.1466 | 0.061532 |
8 | 0.3782 | 0.24934 | 0.1466 | 0.061532 |  |
9 | 0.24934 | 0.1466 | 0.061533 |  |  |
10 | 0.1466 | 0.061533 |  |  |  |
11 | 0.061535 |  |  |  |  |
Table 4: $L^{2}$ relative error, Equation 29, for the boundary flux
$q_{\text{B}\text{B}^{\prime}}(\Omega_{1})$ computed on the unresolved meshes.

####  5.3.2 Transport problem and DMP

Figure 8: Concentration values along two lines for the resolved test cases.
Figure 9: Concentration values along two lines for the unresolved test cases.
Figure 10: Comparison of concentration values along two lines between resolved
and unresolved test cases at time $t=0.01\,\text{[}s]$. Figure 11: Comparison
of concentration values along two lines between resolved and unresolved test
cases at time $t=0.1\,\text{[}s]$. Figure 12: Comparison of concentration
values along two lines between resolved and unresolved test cases at time
$t=0.5\,\text{[}s]$. Figure 13: Concentration distribution at time
$t_{1}=0.01\,\text{[s]}$, $t_{2}=0.10\,\text{[s]}$ , and
$t=T_{\rm{fin}}=0.50\,\text{[s]}$ for the test case $\mathcal{U}_{7}^{1280}$.

For the transport problem, $\Gamma_{D}$, which coincides with $\Gamma_{in}$,
is the left edge of $\Omega$. The boundary condition is set as $g=1.0$ [m-3].
We consider an overall simulation time equal to
$T_{\rm{fin}}=0.50~{}\rm{[s]}$, and set the time step equal to $\Delta
t=0.025~{}\rm{[s]}$.

In Figures 8 and 9, we report concentration profiles along the segments
$\text{C}\text{C}^{\prime}$ and $\text{D}\text{D}^{\prime}$ for some resolved
and unresolved meshes, respectively. While no relevant differences are found
for the unresolved meshes, the numerical results computed on resolved meshes
are very sensitive to the number of subdivisions of the background mesh $be$.
In particular, similar solutions are obtained for meshes with at least
$be=4000$, while all the meshes with lower resolution produce results with
notable differences, especially along the segment $\text{C}\text{C}^{\prime}$.
The results obtained with the unresolved meshes highlight again the importance
of not only resolving the interfaces but also to refine the region close to
them. We point out that the high aspect ratio of the elements located inside
the fracture aperture may affect the stability of the numerical method and
reduce the accuracy of the approximated solutions. Indeed, the numerical
diffusion added in the direction of the longest side causes the unphysical
undershoots visible where the segment $\text{C}\text{C}^{\prime}$ crosses by
the vertical fracture at $x=0.5$.

In Figures 10, 11, and 12, we compare one of the finest resolved test cases
with the numerical solution computed on the unresolved mesh with the largest
number of element layers along the fracture aperture. We analyze the
concentration along lines $\text{C}\text{C}^{\prime}$ and
$\text{D}\text{D}^{\prime}$ at different time instants: $t=0.01$, $t=0.1$ and
$t=T_{\rm{fin}}=0.50$ [s] and observe similar solutions for both the types of
meshes. The good match found between the unresolved and the resolved results
confirms that while the solution computed on the unresolved mesh is mostly
affected by the number AMR, the resolved test cases require a much finer
background to improve the aspect ratio of the elements located inside the
fracture.

In Figure 13, we show the distribution of the concentration at the same times
selected for the profiles depicted in Figure 12. We observe that
concentrations are negligible inside the embedding matrix and that the
physical quantity is mainly transported inside the high permeable fracture
network.

###  5.4 Single 2D

For the second example, we set $L_{1}=L_{2}=100\,[\text{m}]$. Such a square
domain is crossed by an oblique fracture with a thickness
$\delta=0.01\,[\text{m}]$. The subdomain $\Omega_{m}$ consists of two
subdomains, $\Omega_{m1}$ and $\Omega_{m2}$, with different permeabilities,
$k_{m1}$ and $k_{m2}$, and porosities ${\phi_{m1}}$ and ${\phi_{m2}}$. We
report details about the geometrical set up in Table 5, where we denote by
$x_{i}^{s}$ and $x_{i}^{e}$ the coordinates of starting and ending points,
respectively. The material properties are summarized in Table 6.

Table 5: Geometrical sizes employed in the problem “Single fracture network”.
Here $s$ and $e$ refer to the starting and the end position, respectively.
Rock Matrix
---
|  $x_{1}^{s}$ [m] |  $x_{1}^{e}$ [m] |  $x_{2}^{s}$ [m] |  $x_{2}^{s}$ [m]
$\mathrm{\Omega}_{m1}$ |  -50 |  50 |  -50 |  -40
$\mathrm{\Omega}_{m2}$ |  -50 |  50 |  -40 |  -50
$\mathrm{\Gamma_{in}}$ |  0 |  0 |  40 |  50
$\mathrm{\Gamma_{out}}$ |  50 |  50 |  0 |  10
$\mathrm{\Omega}_{A}$ |  -50 |  0 |  -50 |  50
$\mathrm{\Omega}_{B}$ |  0 |  50 |  -50 |  50
Fracture
|  $x_{1}^{s}$ [m] |  $x_{1}^{e}$ [m] |  $x_{2}^{s}$ [m] |  $x_{2}^{s}$ [m]
$\Omega_{f}$ |  -50 |  50 |  -30 |  30
Segments
|  $x_{1}^{s}$ [m] |  $x_{1}^{e}$ [m] |  $x_{2}^{s}$ [m] |  $x_{2}^{s}$ [m]
$\text{A}\text{A}^{\prime}$ |  -50 |  50 |  50 |  -50
$\text{B}\text{B}^{\prime}$ |  -50 |  50 |  30 |  -30
Table 6: Material properties employed in the problem “Single fracture network”. Property | Symbol | Value | Unit
---|---|---|---
Fracture thickness | $\delta$ | 0.01 | m
Matrix conductivity ($\Omega_{m1}$) | $k_{m1}$ | $10^{-5}$ | $\text{m}\times\text{s}$
Matrix conductivity ($\Omega_{m2}$) | $k_{m2}$ | $10^{-6}$ | $\text{m}\times\text{s}$
Fracture conductivity | $k_{f}$ | $10^{-1}$ | $\text{m}\times\text{s}$
Matrix porosity ($\Omega_{m1}$) | ${\phi_{m1}}$ | $0.2$ |
Matrix porosity ($\Omega_{m2}$) | ${\phi_{m2}}$ | $0.25$ |
Fracture porosity | ${\phi_{f}}$ | $0.4$ |

We compare numerical results computed on several resolved and unresolved
meshes. The characteristics of the considered meshes are reported in Table 11
and Table 13, respectively. In Figure 14, we report the resolved mesh
$\mathcal{R}^{114,2}$ and the unresolved mesh $\mathcal{T}^{100,8}$. As the
fracture is oblique, mesh elements of the resolved mesh are skewed, while the
unresolved meshes allow to keep the same aspect ratio for all the elements. We
analyze pressure and concentrations along the segments
$\text{A}\text{A}^{\prime}$ and $\text{B}\text{B}^{\prime}$ depicted in Figure
15, and whose coordinates are reported in Table 5. For the flow problem, we
estimate the approximated flux across the line $L:=\text{C}\text{C}^{\prime}$
located at $x_{1}=0$ [m].

Table 7: Total fluxes, Equation 30, for the resolved and unresolved test cases. All the fluxes values have been multiplied by a scale factor $1.0e^{6}$. Mesh |  $Q_{\text{C}\text{C}^{\prime}}(\Omega_{1})$ |  $Q_{\text{C}\text{C}^{\prime}}(\Omega_{2})$
---|---|---
$\mathcal{R}^{114,2}$ |  $1.86675$ |  $-1.86675$
$\mathcal{R}^{228,4}$ |  $1.85977$ |  $-1.85977$
$\mathcal{R}^{456,8}$ |  $1.85632$ |  $-1.85632$
$\mathcal{R}^{921,16}$ |  $1.85457$ |  $-1.85457$
$\mathcal{U}^{100}_{8}$ |  $1.87065$ |  $-1.87065$
$\mathcal{U}^{200}_{6}$ |  $1.86277$ |  $-1.86277$
$\mathcal{U}^{200}_{7}$ |  $1.86277$ |  $-1.86277$
$\mathcal{U}^{200}_{8}$ |  $1.86174$ |  $-1.86174$
$\mathcal{U}^{400}_{6}$ |  $1.85779$ |  $-1.85779$
$\mathcal{U}^{400}_{7}$ |  $1.85747$ |  $-1.85747$
$\mathcal{U}^{400}_{8}$ |  $1.85534$ |  $-1.85534$
$\mathcal{U}^{800}_{6}$ |  $1.85534$ |  $-1.85534$
$\mathcal{U}^{800}_{7}$ |  $1.85517$ |  $-1.85517$
$\mathcal{U}^{800}_{8}$ |  $1.85505$ |  $-1.85505$
$\mathcal{U}^{1600}_{6}$ |  $1.85411$ |  $-1.85411$
$\mathcal{U}^{1600}_{7}$ |  $1.85399$ |  $-1.85399$

(a) Resolved (b) Unresolved

Figure 14: (a) Example of a resolved mesh ($\mathcal{R}^{114,2}$). (b) Example
of an unresolved mesh ($\mathcal{U}^{100,8}$).

####  5.4.1 Flow problem and local conservation of the flux

For the flow problem, the Dirichlet boundary $\Gamma_{D}$ consists of two
segments $\Gamma_{in}$ and $\Gamma_{out}$, as reported in Table 5, while
$\Gamma_{N}$ is composed of the other remaining sides. We set $g_{in}=4$ [m]
and $g_{out}=1$ [m] on $\Gamma_{in}$ and $\Gamma_{out}$, respectively.
Homogeneous Neumann conditions are imposed on $\Gamma_{N}$. Figure 15 shows
the spatial distribution of the pressure obtained for the test case
$\mathcal{U}^{400}_{6}$.

In Figure 16, we compare the pressure computed along the lines
$\text{A}\text{A}^{\prime}$ and $\text{B}\text{B}^{\prime}$ on some resolved
and unresolved meshes. In particular, we compare coarser mesh solutions with
results computed on meshes fine enough to guarantee mesh-independence. We
observe that solutions present a similar behavior along the segment
$\text{A}\text{A}^{\prime}$, while the results computed along the line
$\text{B}\text{B}^{\prime}$ are more sensitive to the mesh resolution. As
already observed for the benchmark “Regular fracture network”, the resolved
meshes require a certain number of elements in the background, while the
solution computed on the unresolved meshes is mainly affected by the number of
elements close to the matrix-fracture interface.

Figure 15: Pressure distribution for the unresolved test case
$\mathcal{U}_{8}^{320}$. We also report the extrema of the three segments
along which the relevant properties are computed.

In Figure 17, we report the flux along the segment $\text{C}\text{C}^{\prime}$
computed on the meshes $\mathcal{R}^{921,16}$ and $\mathcal{U}^{1600}_{7}$.
Results for both types of meshes present the same profile but with opposite
sign, showing that the two fluxes are matched point-wise and not only in an
integral form. On the other hand, differences between the results obtained
with the two types of mesh are visible close to the two sides of the fracture.
The relative position between the segment $\text{C}\text{C}^{\prime}$ and the
oblique fracture require a large number of either mesh subdivisions, $be$, or
adaptive refinement steps, AMR, to reproduce the steep profile computed for
the resolved test cases. The results reported in Figure 18) for some
unresolved meshes confirm this observation. Indeed, by progressively
increasing the resolution in the fracture, the flux
$q_{\text{C}\text{C}^{\prime}}(\Omega_{2})$ in the central part of the segment
$\text{C}\text{C}^{\prime}$ tends to correctly reproduce the steep
characteristics of the flux.

For the unresolved test cases, we estimate the error by adopting the flux
computed on $\mathcal{U}^{1600}_{7}$ as a reference solution. The results
reported in Table 8 show that, for a fixed number of AMR steps, the error
$e_{q}$ progressively reduces by increasing the number of mesh subdivisions
$be$. In particular, the estimated converge rate presents an order larger than
$1.5$. Similar behavior can be observed fixing the number of mesh subdivisions
and analyzing the error behavior with respect to the number of AMR steps.

In Table 7 we report the value of the total fluxes $Q_{L}(S)$ with
$S\in\\{\Omega_{1},\Omega_{2}\\}$ computed for meshes listed in Tables 11 and
13. Again, all the numerical results confirm that the local conservation is
satisfied for both resolved and unresolved meshes.

Figure 16: Comparison of pressure values along two lines between resolved and unresolved test cases. Table 8: $L_{2}$ relative error, Equation 29, for the boundary flux $q_{\text{C}\text{C}^{\prime}}(\Omega_{1})$ computed on the unresolved meshes. be: amr: |  6 |  7 |  8
---|---|---|---
$100$ |  \- - |  \- - |  0.004862
$200$ |  0.009018 |  0.004894 |  0.003498
$400$ |  0.004911 |  0.003508 |  0.002235
$800$ |  0.003512 |  0.002238 |  0.001351
$1600$ |  0.002239 |  \- - |  \- -
Figure 17: Boundary fluxes, $q_{\text{C}\text{C}^{\prime}}(\Omega_{1})$ and
$q_{\text{C}\text{C}^{\prime}}(\Omega_{2})$, across boundary
$\text{C}\text{C}^{\prime}$ for the test cases $\mathcal{R}^{921,16}$ and
$\mathcal{U}^{1600}_{7}$. Figure 18: Boundary flux
$q_{\text{C}\text{C}^{\prime}}(\Omega_{2})$ across the region of the boundary
$\text{C}\text{C}^{\prime}$ intersecting the fracture for the test cases
$\mathcal{U}^{100,8}$, $\mathcal{U}^{200,8}$ and $\mathcal{U}^{1600}_{7}$.

#### Transport problem and DMP

Figure 19: Comparison of concentration values along two lines between two
selected resolved and unresolved test cases at time $t=T=10^{9}\,\text{[s]}$.
Figure 20: Concentration distribution at time $t_{1}=10^{8}\,[s]$,
$t_{2}=5\,10\,^{8}[s]$ , and $t=T=\,10^{9}\,[s]$ for the test case
$\mathcal{U}^{800}_{6}$.

For the transport problem, $\Gamma_{D}$ coincides with $\Gamma_{in}$, on which
we impose $g=0.01$ [m-3]. The rest of the boundary is $\Gamma_{out}$. We set
$T_{\rm{fin}}=10^{9}$ [s] and $\Delta t=10^{2}$ [s]. We compare the
concentrations obtained on some resolved and unresolved meshes by analyzing
their value at the final time $t=T_{\rm{fin}}$ along the lines
$\text{A}\text{A}^{\prime}$ and $\text{B}\text{B}^{\prime}$. The results are
reported in Figure 19.

Although for the flow problem results coincide between the two types of
meshes, the concentration profiles computed on $\mathcal{R}^{921,16}$ and
$\mathcal{U}^{800}_{6}$ present different behaviours close to the point
$\text{A}^{\prime}$ inside $\Omega_{m}$ and close to the point
$\text{B}^{\prime}$ inside $\Omega_{f}$. Such differences are related to the
skewness of the elements of the resolved meshes which introduce a larger
stabilization to preserve the positivity of the solution but reduce its
accuracy. Moreover, the aspect ratio of the bilinear elements deteriorates
when increasing the number of subdivision $be$. This observation is supported
by the good agreement found between the concentrations computed on the
unresolved meshes and the one computed on $\mathcal{R}^{114,2}$.

In Figure 20, we show the distribution of the concentration $c$ at the time
instants $t=10^{8}$, $t=5\,10\,^{8}$ , and $t=T_{\rm{fin}}=\,10^{9}$ [s].
Here, the values of the concentrations are always comprised between $0$ and
$0.01$ for the entire simulated time, and verify the ability of the AFC scheme
to ensure the DMP and to preserve the positivity of the solution.

###  5.5 Realistic fracture network

This example has been presented in [56], where a comparison between different
discretizations based on a one-dimensional representation of the fractures has
been performed. For this numerical example, we have $L_{1}=700$ [m] and
$L_{2}=600$ [m] and the set $\mathcal{F}$ is composed of $63$ fractures,
ranging from isolated fractures up to groups of tens of fractures each. The
aperture of the fractures is $10^{-2}$ [m]. Indeed, such a geometry represents
an outcrop in the Sotra island, near Bergen in Norway. The matrix permeability
is $k_{m}=10^{-14}$ [m$\times$s-1], while the fracture permeability is to
$k_{f}=10^{-8}$ m$\times$s-1. We consider no-flow boundary condition on top
and bottom, pressure $1013250$ [m] on the left, and pressure $0$ [m] on the
right of the boundary of the domain. Due to the high geometrical complexity of
the fracture network, meshes which resolve the interface $\Gamma$ are
difficult to be generated. On the other hand, AMR method enables a fully
automatic generation of meshes for domains that contain a large number of
highly complex heterogeneities.

As we have discussed in Section 2 and we will show in the numerical results,
the DMP may not hold for FE discretizations of diffusion problems on adapted
meshes. We will show that the proposed stabilization instead ensures the DMP
and provides a suitable discretization method.

The characteristics of the considered meshes are reported in Table 14. The
superscript be refers to the number of subdivisions along the horizontal
direction. The number of intervals along the vertical direction is
$6\text{be}/7$.

Table 9: Maximum and the minimum values for the pressure for the benchmark “Realistic fracture networks”. Mesh | $\min{p}^{h}$ | $\max{p}^{h}$
---|---|---
$\mathcal{U}^{7}_{7}$ | -3255.8181 | 1013250
$\mathcal{U}^{7}_{8}$ | -2.7606e-14 | 1014210.2276
$\mathcal{U}^{7}_{9}$ | -100.4391 | 1013250
$\mathcal{U}^{7}_{10}$ | -10.6924 | 1013250
$\mathcal{U}^{14}_{6}$ | -3255.8181 | 1013250
$\mathcal{U}^{14}_{7}$ | -2.7606e-14 | 1014210.2276
$\mathcal{U}^{14}_{8}$ | -100.4391 | 1013250
$\mathcal{U}^{14}_{9}$ | -10.6924 | 1013250
$\mathcal{U}^{28}_{5}$ | -3255.8451 | 1013250
$\mathcal{U}^{28}_{6}$ | -2.7607e-14 | 1014210.112
$\mathcal{U}^{28}_{7}$ | -100.4418 | 1013250
$\mathcal{U}^{28}_{8}$ | -10.6928 | 1013250

(a) Step 0 (b) Step 1

(c) Step 2 (d) Step 4

Figure 21: Sequences of meshes generated by iterating several times the AMR
procedure.

In Table 9, we report the minimum and the maximum for a non-stabilized FE
discretization implemented on the corresponding grid. We observe that for all
the cases, the DMP is violated, as the solution provided either a minimum
smaller than $0$ or a maximum larger than the value imposed on the left side
of the domain. On the other hand, the stabilized discretization provides a
solution that attains its maximum and minimum value at the boundary. In Figure
22, we report the solutions obtained with the stabilized and non-stabilized
methods on the mesh $\mathcal{U}^{28}_{8}$. On the non-stabilized solution
(a), we superimpose the contour lines of the stabilized solution (b) to
highlight the differences between the two solutions. We observe that the
largest differences are found in the central part of the domain.

(a) Unresolved Non Stabilized (b) Unresolved Stabilized

Figure 22: Pressure distribution for the unresolved test case
$\mathcal{U}^{28}_{8}$. We also report the extrema of the three segments along
which the relevant properties are computed.

In Table 10, we report the relative $L^{2}$ errors (i.e., replacing
$q_{L}(\Omega_{i})$ with $p^{h}$ in Equation 29) for the two methods with
respect to a reference solution computed on a mesh $\mathcal{U}^{56}_{11}$
with a non-stabilized method. The magnitude of the errors computed with the
stabilized and non-stabilized methods are comparable, with the stabilized one
providing, in general, a bit larger errors but ensuring the DMP. The addition
of the algebraic stabilization is in general non-consistent, unless non-linear
stabilization methods are considered [64].

We report in Figure 23 the pressure solution along two different lines:
$\text{A}\text{A}^{\prime}$, $y=500$ [m], and $\text{B}\text{B}^{\prime}$,
$x=625$ [m], for the meshes with $\text{be}=28$. We note that except for the
case with $\text{amr}=5$, the solutions behave similarly, showing again that
the behavior of the solutions is mostly dominated by the number of elements in
the fracture.

Finally, in Figure 24, we report a comparison between our stabilized method
and the methods based on the hybrid-dimensional representation of the domain
presented in [56]. Our stabilized method based on AMR provides a possible
alternative to hybrid methods for computing solutions with a two-dimensional
representation of the fractures. As we see, all the characteristics of the
solutions along these two lines are correctly described.

Table 10: $L_{2}$ relative error, Equation 29, for the boundary flux $q_{\text{B}\text{B}^{\prime}}(\Omega_{1})$ computed on the resolved meshes. | Non-stabilized
---|---
amr be | 7 | 14 | 28
5 |  |  | 0.31745
6 |  | 0.31745 | 0.024124
7 | 0.31745 | 0.024133 | 0.020829
8 | 0.024133 | 0.020821 | 0.0016197
9 | 0.020821 | 0.0016153 |
10 | 0.0016153 |  |
| Stabilized
amr be | 7 | 14 | 28
5 |  |  | 0.31604
6 |  | 0.31635 | 0.021175
7 | 0.31635 | 0.020986 | 0.027017
8 | 0.020986 | 0.026737 | 0.0078999
9 | 0.026737 | 0.0074733 |
10 | 0.0074733 |  |

Figure 23: Comparison of pressure values along two lines between stabilized
and non stabilized unresolved test cases.  Figure 24: Comparison of pressure
values along two lines between the unresolved test case and the numerical
solutions based on the hybrid-dimensional approaches proposed in [56].

##  6 Conclusions

We have presented a continuous finite element method for the numerical
simulation of the coupled flow and transport problems in fractured porous
media based on an equi-dimensional representation of the fractures and an
adaptive mesh refinement. Adaptive mesh refinement allows to automatically and
efficiently adapt an initially uniform mesh to any kind of heterogeneities.
The resulting non-conforming mesh does not resolve the matrix-fracture
interfaces which can be approximated with the desired accuracy.

For the proposed finite element discretization on non-conforming meshes, we
have extended standard results on the discrete maximum principle, proposed a
stabilization strategy based on algebraic diffusion operators, and shown that
global and local conservation properties hold. The stabilization strategy is
based on the algebraic flux correction method which allows us to ensure the
discrete maximum principle and to remove the unphysical oscillations in the
solution.

We have presented a numerical comparison with a standard finite element method
based conforming meshes which resolve the interfaces. Such numerical examples
highlight the main advantages in terms of accuracy and solution properties
deriving from the use of meshes with non-skewed elements. In particular, the
use of adapted meshes allows to avoid large aspect ratios for the elements
inside the fracture aperture, which represents the main drawback for equi-
dimensional models. On the other hand, the numerical results illustrate that
the number of elements inside and around the fracture has a significant
influence on the accuracy of the solution. As the last contribution, we have
shown that the proposed discretization method provides accurate results also
for realistic cases with complex fracture networks.

It is worth to point out that the proposed approach allows describing the
realistic equidimensional structure of fracture networks without requiring
their explicitly representation.

In the future we plan to employ such a discretization method to validate the
approaches based on a lower-dimensional representation of fractures, which are
usually employed to perform simulation in fractured media.

The study of the accuracy of the proposed method for realistic three-
dimensional cases should also be performed. As the proposed method does not
involve the coupling of different discretizations for the matrix and the
background, the linear systems to be solved are characterized by positive
definite M-matrices and, hence, are suitable for the efficient solution of the
problem with multigrid methods. This is another fundamental characteristic of
the proposed method. In this context, the resulting hierarchy of meshes could
also be considered in the construction of the levels of a multigrid solver.
Moreover, we marked for refinement the elements which presented a non-empty
overlap with at least one fractures. This idea comes on the a-priori
information that the discretization error in the simulation of heterogenous
media is usually localized close to the interfaces. Obviously, the development
of suitable error estimators would also provide a strategy to control such
numerical errors.

## Acknowledgments

M. G. C. Nestola acknowledges gratefully the support the PASC Project FASTER.
Marco Favino acknowledges gratefully the support of the Swiss National Science
Foundation (SNSF) through the grant PZ00P2$\\_$180112\. All methods and
routines used in this study are implemented within the open-source software
library _Parrot_ : https://github.com/favinom/parrot/. _Parrot_ ’s lead
developers are Maria GC Nestola<EMAIL_ADDRESS>and Marco Favino
(marco.favino@unil.ch).

##  7 Appendix A. Mesh Characteristics

Mesh characteristics of the resolved and the unresolved meshes employed in
Sec, 5.

Table 11: Resolved meshes for the benchmark “Single fracture network”. Mesh | $n^{E}$ | $n^{r}$ | $w$
---|---|---|---
$\mathcal{R}^{114,2}$ | 6 612 | 6 785 | 2
$\mathcal{R}^{228,4}$ | 26 448 | 27 793 | 4
$\mathcal{R}^{456,8}$ | 105 792 | 106 481 | 8
$\mathcal{R}^{921,16}$ | 423 168 | 424 545 | 16
Table 12: Unresolved meshes for the benchmark “Regular fracture network”. Mesh | ${E}_{tot}$ | ${N}_{tot}$ | $N_{tot}^{r}$ | $w$
---|---|---|---|---
$\mathcal{U}^{80}_{7}$ | 219 256 | 254 851 | 184 067 | 1.02
$\mathcal{U}^{80}_{8}$ | 434 224 | 505 641 | 363 225 | 2.05
$\mathcal{U}^{80}_{9}$ | 1 720 768 | 1863685 | 1578389 | 4.10
$\mathcal{U}^{80}_{10}$ | 4 300 600 | 4586823 | 4014951 | 8.19
$\mathcal{U}^{80}_{11}$ | 12 900 400 | 13473253 | 12328181 | 16.38
$\mathcal{U}^{160}_{6}$ | 236 848 | 272341 | 202069 | 1.02
$\mathcal{U}^{160}_{7}$ | 451 816 | 523131 | 381227 | 2.05
$\mathcal{U}^{160}_{8}$ | 1 736 896 | 1879729 | 1594897 | 4.10
$\mathcal{U}^{160}_{9}$ | 431 6728 | 4602867 | 4031459 | 8.19
$\mathcal{U}^{160}_{10}$ | 12 916 528 | 13489297 | 12344689 | 16.38
$\mathcal{U}^{320}_{5}$ | 310 360 | 345631 | 276431 | 1.02
$\mathcal{U}^{320}_{6}$ | 525 328 | 596421 | 455589 | 2.05
$\mathcal{U}^{320}_{7}$ | 1 807 264 | 1949893 | 1666085 | 4.10
$\mathcal{U}^{320}_{8}$ | 4 387 096 | 4673031 | 4102647 | 8.19
$\mathcal{U}^{320}_{9}$ | 12 986 896 | 13559461 | 12415877 | 16.38
$\mathcal{U}^{640}_{4}$ | 610 912 | 645721 | 578713 | 1.02
$\mathcal{U}^{640}_{5}$ | 825 880 | 896511 | 757871 | 2.05
$\mathcal{U}^{640}_{6}$ | 2101 312 | 2243497 | 1961833 | 4.10
$\mathcal{U}^{640}_{7}$ | 4 681 144 | 4966635 | 4398395 | 8.19
$\mathcal{U}^{640}_{8}$ | 13 280 944 | 13853065 | 12711625 | 16.38
$\mathcal{U}^{1280}_{3}$ | 1 826 344 | 1860211 | 1797635 | 1.02
$\mathcal{U}^{1280}_{4}$ | 2 041 312 | 2111001 | 1976793 | 2.05
$\mathcal{U}^{1280}_{5}$ | 3 303 520 | 3444781 | 316750 | 4.10
$\mathcal{U}^{1280}_{6}$ | 5 883 352 | 6167919 | 5604063 | 8.19
$\mathcal{U}^{1280}_{7}$ | 14 483 152 | 15054349 | 13917293 | 16.38
$\mathcal{U}^{2560}_{2}$ | 6 714 736 | 6746701 | 6693037 | 1.02
$\mathcal{U}^{2560}_{3}$ |  6 929 704 | 6997491 | 6872195 | 2.05
$\mathcal{U}^{2560}_{4}$ |  8 165 248 | 8304625 | 8036209 | 4.10
$\mathcal{U}^{2560}_{5}$ | 10 745 080 | 11027763 | 10472771 | 8.19
$\mathcal{U}^{2560}_{6}$ | 19 344 880 | 19914193 | 18786001 | 16.38
Table 13: Unresolved meshes for the benchmark “single fracture network”. Mesh | $n^{E}$ | $n^{H}$ | $n^{r}$ | $w$
---|---|---|---|---
$\mathcal{U}^{100}_{8}$ | 254 740 | 295 757 | 214 173 | 2.56
$\mathcal{U}^{200}_{6}$ | 162 628 | 141 640 | 121 903 | 1.28
$\mathcal{U}^{200}_{7}$ | 283 720 | 324 775 | 243 509 | 2.56
$\mathcal{U}^{200}_{8}$ | 698 200 | 780 233 | 617 035 | 5.12
$\mathcal{U}^{400}_{6}$ | 442 813 | 401 680 | 362 185 | 2.56
$\mathcal{U}^{400}_{7}$ | 815 920 | 898 029 | 735 471 | 5.12
$\mathcal{U}^{400}_{8}$ | 2 037 040 | 2 201 075 | 1 874 681 | 10.24
$\mathcal{U}^{800}_{6}$ | 1 291 360 | 1 373 625 | 1 212 347 | 5.12
$\mathcal{U}^{800}_{7}$ | 2 512 480 | 2 678 304 | 2 351 557 | 10.24
$\mathcal{U}^{800}_{8}$ | 6 890 080 | 7 218 127 | 6 565 333 | 20.48
$\mathcal{U}^{1600}_{6}$ | 4 423 360 | 4 587 867 | 4 265 313 | 10.24
$\mathcal{U}^{1600}_{7}$ | 8 800 960 | 9 129 323 | 8 479 089 | 20.48
Table 14: Unresolved meshes for the benchmark “realistic fracture networks”. Mesh | $n^{E}$ | $n^{H}$ | $n^{r}$ | $w$
---|---|---|---|---
$\mathcal{U}^{7}_{7}$ | 62232 | 75439 | 49251 | 0.01
$\mathcal{U}^{7}_{8}$ | 133128 | 162364 | 104142 | 0.03
$\mathcal{U}^{7}_{9}$ | 278223 | 340375 | 216354 | 0.06
$\mathcal{U}^{7}_{10}$ | 573276 | 701910 | 444958 | 0.12
$\mathcal{U}^{14}_{6}$ | 62232 | 75439 | 49251 | 0.01
$\mathcal{U}^{14}_{7}$ | 133128 | 162364 | 104142 | 0.03
$\mathcal{U}^{14}_{8}$ | 278223 | 340375 | 216354 | 0.06
$\mathcal{U}^{14}_{9}$ | 573276 | 701910 | 444958 | 0.12
$\mathcal{U}^{28}_{5}$ | 62286 | 75492 | 49321 | 0.01
$\mathcal{U}^{28}_{6}$ | 133182 | 162417 | 104212 | 0.03
$\mathcal{U}^{28}_{7}$ | 278277 | 340428 | 216424 | 0.06
$\mathcal{U}^{28}_{8}$ | 573330 | 701963 | 445028 | 0.12

## References

  * [1] J. W. Tester, B. Anderson, A. Batchelor, D. Blackwell, R. DiPippo, E. Drake, J. Garnish, B. Livesay, M. C. Moore, K. Nichols, et al., The future of geothermal energy: Impact of enhanced geothermal systems (egs) on the united states in the 21st century, Massachusetts Institute of Technology 209.
  * [2] M. W. McClure, R. N. Horne, Correlations between formation properties and induced seismicity during high pressure injection into granitic rock, Engineering Geology 175 (2014) 74–80.  doi:http://dx.doi.org/10.1016/j.enggeo.2014.03.015.
  * [3] C. E. Bond, R. Wightman, P. S. Ringrose, The influence of fracture anisotropy on co2 flow, Geophysical Research Letters 40 (7) (2003) 1284–1289.  doi:10.1002/grl.50313.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/grl.50313

  * [4] E. Bonnet, O. Bour, N. E. Odling, P. Davy, I. Main, P. Cowie, B. Berkowitz, Scaling of fracture systems in geological media, Reviews of Geophysics 39 (3) (2001) 347–383.  doi:10.1029/1999RG000074.
URL http:https://doi.org/10.1029/1999RG000074

  * [5] A. Rasmuson, I. Neretnieks, Radionuclide transport in fast channels in crystalline rock, Water Resources Research 22 (8) (1986) 1247–1256.
  * [6] F. Amann, V. Gischig, K. Evans, J. Doetsch, R. Jalali, B. Valley, H. Krietsch, N. Dutler, L. Villiger, B. Brixel, et al., The seismo-hydromechanical behavior during deep geothermal reservoir stimulations: open questions tackled in a decameter-scale in situ stimulation experiment, Solid Earth 9 (1) (2018) 115–137.
  * [7] I. Berre, F. Doster, E. Keilegavlen, Flow in fractured porous media: A review of conceptual models and discretization approaches, Transport in Porous Media 130 (1) (2019) 215–236.
  * [8] B. Eikemo, K.-A. Lie, G. T. Eigestad, H. K. Dahle, Discontinuous galerkin methods for advective transport in single-continuum models of fractured media, Advances in water resources 32 (4) (2009) 493–506.
  * [9] R. Liu, B. Li, Y. Jiang, N. Huang, Mathematical expressions for estimating equivalent permeability of rock fracture networks, Hydrogeology Journal 24 (7) (2016) 1623–1649.
  * [10] E. T. Chung, Y. Efendiev, T. Leung, M. Vasilyeva, Coupling of multiscale and multi-continuum approaches, GEM-International Journal on Geomathematics 8 (1) (2017) 9–41.
  * [11] M. Karimi-Fard, L. J. Durlofsky, K. Aziz, et al., An efficient discrete fracture model applicable for general purpose reservoir simulators, in: SPE Reservoir Simulation Symposium, Society of Petroleum Engineers, 2003.
  * [12] T. H. Sandve, I. Berre, J. M. Nordbotten, An efficient multi-point flux approximation method for discrete fracture–matrix simulations, Journal of Computational Physics 231 (9) (2012) 3784–3800.
  * [13] R. Ahmed, M. G. Edwards, S. Lamine, B. A. Huisman, M. Pal, Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model, Journal of Computational Physics 284 (2015) 462–489.
  * [14] P. Angot, F. Boyer, F. Hubert, Asymptotic and numerical modelling of flows in fractured porous media, ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique 43 (2) (2009) 239–275.
  * [15] A. Fumagalli, E. Keilegavlen, S. Scialò, Conforming, non-conforming and non-matching discretization couplings in discrete fracture network simulations, Journal of Computational Physics 376 (2019) 694–712.
  * [16] P. F. Antonietti, C. Facciolà, M. Verani, Mixed-primal discontinuous galerkin approximation of flows in fractured porous media on polygonal and polyhedral grids.
  * [17] M. Köppel, V. Martin, J. E. Roberts, A stabilized lagrange multiplier finite-element method for flow in porous media with fractures, GEM-International Journal on Geomathematics 10 (1) (2019) 7.
  * [18] P. Schädle, P. Zulian, D. Vogler, S. R. Bhopalam, M. G. Nestola, A. Ebigbo, R. Krause, M. O. Saar, 3d non-conforming mesh model for flow in fractured porous media using lagrange multipliers, Computers & Geosciences 132 (2019) 42 – 55.  doi:https://doi.org/10.1016/j.cageo.2019.06.014.
URL http://www.sciencedirect.com/science/article/pii/S0098300418312032

  * [19] N. Frih, V. Martin, J. E. Roberts, A. Saâda, Modeling fractures as interfaces with nonmatching grids, Computational Geosciences 16 (4) (2012) 1043–1060.  doi:10.1007/s10596-012-9302-6.
URL http://link.springer.com/10.1007/s10596-012-9302-6

  * [20] W. M. Boon, J. M. Nordbotten, I. Yotov, Robust discretization of flow in fractured porous media, SIAM Journal on Numerical Analysis 56 (4) (2018) 2203–2233.
  * [21] E. Burman, P. Hansbo, M. G. Larson, A simple finite element method for elliptic bulk problems with embedded surfaces, Computational Geosciences 23 (1) (2019) 189–199.
  * [22] L. Li, S. H. Lee, et al., Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media, SPE Reservoir Evaluation & Engineering 11 (04) (2008) 750–758.
  * [23] L. H. Odsæter, T. Kvamsdal, M. G. Larson, A simple embedded discrete fracture–matrix model for a coupled flow and transport problem in porous media, Computer Methods in Applied Mechanics and Engineering 343 (2019) 572–601.
  * [24] L. Formaggia, A. Fumagalli, A. Scotti, P. Ruffo, A reduced model for darcy’s problem in networks of fractures, ESAIM: Mathematical Modelling and Numerical Analysis 48 (4) (2014) 1089–1116.
  * [25] N. Schwenck, B. Flemisch, R. Helmig, B. I. Wohlmuth, Dimensionally reduced flow models in fractured porous media: crossings and boundaries, Computational Geosciences 19 (6) (2015) 1219–1230.
  * [26] B. Flemisch, A. Fumagalli, A. Scotti, A review of the xfem-based approximation of flow in fractured porous media, in: Advances in Discretization Methods, Springer, 2016, pp. 47–76.
  * [27] E. Burman, P. Hansbo, M. G. Larson, A cut finite element method for a model of pressure in fractured media, arXiv preprint arXiv:2006.01521.
  * [28] E. Burman, P. Hansbo, M. G. Larson, K. Larsson, Cut finite elements for convection in fractured domains, Computers & Fluids 179 (2019) 726–734.
  * [29] P. F. Antonietti, C. Facciolà, M. Verani, Polytopic discontinuous galerkin methods for the numerical modelling of flow in porous media with networks of intersecting fractures, arXiv preprint arXiv:2002.06420.
  * [30] P. F. Antonietti, L. Formaggia, A. Scotti, M. Verani, N. Verzott, Mimetic finite difference approximation of flows in fractured porous media, ESAIM: Mathematical Modelling and Numerical Analysis 50 (3) (2016) 809–832.
  * [31] L. Formaggia, A. Scotti, F. Sottocasa, Analysis of a mimetic finite difference approximation of flows in fractured porous media, ESAIM: Mathematical Modelling and Numerical Analysis 52 (2) (2018) 595–630.
  * [32] A. Fumagalli, E. Keilegavlen, Dual virtual element methods for discrete fracture matrix models, Oil & Gas Science and Technology–Revue d’IFP Energies nouvelles 74 (2019) 41.
  * [33] W. M. Boon, J. M. Nordbotten, I. Yotov, Robust discretization of flow in fractured porous media, SIAM Journal on Numerical Analysis 56 (4) (2018) 2203–2233.
  * [34] M. G. Larson, A. J. Niklasson, A conservative flux for the continuous galerkin method based on discontinuous enrichment, Calcolo 41 (2) (2004) 65–76.
  * [35] V. Martin, J. Jaffré, J. E. Roberts, Modeling fractures and barriers as interfaces for flow in porous media, SIAM Journal on Scientific Computing 26 (5) (2005) 1667–1691.
  * [36] C. D’Angelo, A. Scotti, A mixed finite element method for darcy flow in fractured porous media with non-matching grids, ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique 46 (2) (2012) 465–489.
  * [37] F. Chave, D. A. Di Pietro, L. Formaggia, A hybrid high-order method for darcy flows in fractured porous media, SIAM Journal on Scientific Computing 40 (2) (2018) A1063–A1094.
  * [38] J. Douglas Jr, J. E. Roberts, Mixed finite element methods for second order elliptic problems, Mat. Apl. Comput 1 (1) (1982) 91–103.
  * [39] B. Riviè, M. F. Wheeler, K. Banaś, et al., Part ii. discontinuous galerkin method applied to a single phase flow in porous media, Computational Geosciences 4 (4) (2000) 337–349.
  * [40] S. Sun, J. Liu, A locally conservative finite element method based on piecewise constant enrichment of the continuous galerkin method, SIAM Journal on Scientific Computing 31 (4) (2009) 2528–2548.
  * [41] J. Wheeler, Simulation of heat transfer from a warm pipe buried in permafrost, 74th, in: National Meeting of the American Institute of Chemical Engineers, New Orleans,(March, 1973).
  * [42] G. F. Carey, Some further properties of the superconvergent flux projection, Communications in numerical methods in engineering 18 (4) (2002) 241–250.
  * [43] B. Cockburn, J. Gopalakrishnan, H. Wang, Locally conservative fluxes for the continuous galerkin method, SIAM Journal on Numerical Analysis 45 (4) (2007) 1742–1776.
  * [44] L. H. Odsæter, M. F. Wheeler, T. Kvamsdal, M. G. Larson, Postprocessing of non-conservative flux for compatibility with transport in heterogeneous media, Computer Methods in Applied Mechanics and Engineering 315 (2017) 799–830.
  * [45] T. J. Hughes, G. Engel, L. Mazzei, M. G. Larson, The continuous galerkin method is locally conservative, Journal of Computational Physics 163 (2) (2000) 467–488.
  * [46] P. G. Ciarlet, P.-A. Raviart, Maximum principle and uniform convergence for the finite element method, Computer Methods in Applied Mechanics and Engineering 2 (1) (1973) 17–31.
  * [47] B. Fischer, A. Ramage, D. Silvester, A. Wathen, On parameter choice and iterative convergence for stabilised discretisations of advection–diffusion problems, Computer methods in applied mechanics and engineering 179 (1-2) (1999) 179–195.
  * [48] T. J. Hughes, A multidimentional upwind scheme with no crosswind diffusion, Finite Element Methods for Convection Dominated Flows, AMD 34.
  * [49] F. Brezzi, L. P. Franca, A. Russo, Further considerations on residual-free bubbles for advective-diffusive equations, Computer Methods in Applied Mechanics and Engineering 166 (1-2) (1998) 25–33.
  * [50] V. John, E. Schmeyer, Finite element methods for time-dependent convection–diffusion–reaction equations with small diffusion, Computer methods in applied mechanics and engineering 198 (3-4) (2008) 475–494.
  * [51] D. Kuzmin, Explicit and implicit fem-fct algorithms with flux linearization, Journal of Computational Physics 228 (7) (2009) 2517–2534.
  * [52] M. Köppel, V. Martin, J. E. Roberts, A stabilized lagrange multiplier finite-element method for flow in porous media with fractures, GEM-International Journal on Geomathematics 10 (1) (2019) 7.
  * [53] Z. Xu, Y. Yang, The hybrid dimensional representation of permeability tensor: A reinterpretation of the discrete fracture model and its extension on nonconforming meshes, Journal of Computational Physics (2020) 109523.
  * [54] M. Favino, J. Hunziker, E. Caspari, B. Quintal, K. Holliger, R. Krause, Fully-automated adaptive mesh refinement for media embedding complex heterogeneities: application to poroelastic fluid pressure diffusion, Computational Geosciences (2019) 1–20.
  * [55] J. Bonilla, S. Badia, Monotonicity-preserving finite element schemes with adaptive mesh refinement for hyperbolic problems, Journal of Computational Physics (2020) 109522.
  * [56] B. Flemisch, I. Berre, W. Boon, A. Fumagalli, N. Schwenck, A. Scotti, I. Stefansson, A. Tatomir, Benchmarks for single-phase flow in fractured porous media, Advances in Water Resources 111 (2018) 239–258.
  * [57] I. Berre, W. M. Boon, B. Flemisch, A. Fumagalli, D. Gläser, E. Keilegavlen, A. Scotti, I. Stefansson, A. Tatomir, K. Brenner, et al., Verification benchmarks for single-phase flow in three-dimensional fractured porous media, arXiv preprint arXiv:2002.07005.
  * [58] D. Kuzmin, Algebraic flux correction i, in: Flux-corrected transport, Springer, 2012, pp. 145–192.
  * [59] Y. Saito, T. Mitsui, Stability analysis of numerical schemes for stochastic differential equations, SIAM Journal on Numerical Analysis 33 (6) (1996) 2254–2267.
  * [60] R. S. Varga, On a discrete maximum principle, SIAM Journal on Numerical Analysis 3 (2) (1966) 355–359.
  * [61] I. Berre, W. M. Boon, B. Flemisch, A. Fumagalli, D. Gläser, E. Keilegavlen, A. Scotti, I. Stefansson, A. Tatomir, K. Brenner, et al., Verification benchmarks for single-phase flow in three-dimensional fractured porous media, arXiv preprint arXiv:2002.07005.
  * [62] S. Geiger, M. Dentz, I. Neuweiler, et al., A novel multi-rate dual-porosity model for improved simulation of fractured and multiporosity reservoirs, SPE journal 18 (04) (2013) 670–684.
  * [63] L. H. Odsæter, T. Kvamsdal, M. G. Larson, A simple embedded discrete fracture–matrix model for a coupled flow and transport problem in porous media, Computer Methods in Applied Mechanics and Engineering 343 (2019) 572–601.
  * [64] D. Kuzmin, M. Shashkov, D. Svyatskiy, A constrained finite element method satisfying the discrete maximum principle for anisotropic diffusion problems, Journal of Computational Physics 228 (9) (2009) 3448 – 3463.  doi:https://doi.org/10.1016/j.jcp.2009.01.031.
URL http://www.sciencedirect.com/science/article/pii/S0021999109000552