Modelling and discretization of flow in porous media with thin, full-tensor
permeability inclusions
M. Starnoni1,2, I. Berre1, E. Keilegavlen1, & J.M. Nordbotten1
1Department of Mathematics, University of Bergen, Bergen, Norway
2Department of Environment, Land and Infrastructure Engineering, Politecnico
di Torino, Torino, Italy

## Abstract

When modelling fluid flow in fractured reservoirs, it is common to represent
the fractures as lower-dimensional inclusions embedded in the host medium.
Existing discretizations of flow in porous media with thin inclusions assume
that the principal directions of the inclusion permeability tensor are aligned
with the inclusion orientation. While this modelling assumption works well
with tensile fractures, it may fail in the context of faults, where the damage
zone surrounding the main slip surface may introduce anisotropy that is not
aligned with the main fault orientation. In this paper, we introduce a
generalized dimensional reduced model which preserves full-tensor permeability
effects also in the out-of-plane direction of the inclusion. The governing
equations of flow for the lower-dimensional objects are obtained through
vertical averaging. We present a framework for discretization of the resulting
mixed-dimensional problem, aimed at easy adaptation of existing simulation
tools. We give numerical examples that show the failure of existing
formulations when applied to anisotropic faulted porous media, and go on to
show the convergence of our method in both 2D and 3D.

Key points

  * •

Existing local discretizations of flow in fractured porous media fail in
modelling out-of plane anisotropic properties of thin inclusions

  * •

We present a new framework to modelling and discretizing flow in porous media
with thin, full-tensor permeability inclusions

  * •

We show convergence of our method in both 2D and 3D faulted porous media

Keywords discretization, faults, permeability, mixed-dimensional, flow, porous
media

##  1 Introduction

Modeling and simulation of flow in porous media with faults, fractures, and
other thin inclusions representing discontinuities is central to a wide range
of subsurface engineering applications, including geothermal energy
exploitation (Bödvarsson and Tsang,, 1982), shale gas extraction (Cao et al.,,
2016), carbon sequestration (Johnson et al.,, 2009), and energy storage
(Nagelhout and Roest,, 1997).

The inclusions are characterized by a high aspect ratio, and permeability
significantly different from that of the host medium; hence, they severely
affect flow patterns. This poses a challenge for traditional simulation
models, which are based on upscaling of fine-scale details into an equivalent
permeability (Oda,, 1985; Farmer,, 2002; Liu et al.,, 2016; Sævik et al.,,
2013). We instead focus on an alternative approach, which explicitly
represents the inclusions in the mathematical and simulation models and
thereby to a large degree avoids challenges relating to parameter upscaling.
To avoid elongated cells at the inclusion in the computational grid, it is
common to represent the inclusions as co-dimension one objects embedded in the
host medium (Boon et al.,, 2018; Nordbotten et al.,, 2019). The intersection
of inclusions further gives rise to line and point intersections of co-
dimension two and three. Each of these objects (matrix, inclusions, and
intersection points and lines) are represented as independent subdomains
separated by interfaces. We refer to this representation of the geometry as
mixed-dimensional.

Governing equations for fluid flow in lower-dimensional representation of the
inclusion can be derived by integration in the direction orthogonal to the
inclusion. This leads to a decomposition of the governing equations into an
in-plane component that represents flow within the inclusion, and an out-of-
plane component that couples flow between the inclusion and the host medium.
While the in-plane flow has been modeled with both linear and non-linear, as
well as both isotropic and non-isotropic flow models (Martin et al.,, 2005;
Reichenberger et al.,, 2006; Brenner et al.,, 2017, 2018), existing models for
the coupling term are limited by an assumption on orthogonal flow between
inclusion and host. Reduced order models for flow were also developed for
aquifers, leading to the same set of equations, see for instance Bear, (1979),
Yortsos, (1995), and Nordbotten and Celia, (2011). These existing models will
be denoted as ”local” in the following, meaning that each partial differential
equation (PDE) contains only quantities associated with the subdomain where
the PDE is defined.

Local models generally work well when the inclusion is a joint (tensile
fracture). However, inclusions with a more complex geological history may have
significantly more complex flow properties in the out-of-plane direction. For
instance, the damage zone in the vicinity of faults may exhibit shear
fractures, slip surfaces, and/or deformation bands, as summarized in Fossen et
al., (2007). These features introduce secondary permeability anisotropy in the
damage zone as they tend to have preferred orientations, as shown by both
field studies (Fossen et al.,, 2005; Johansen and Fossen,, 2008) and core
analysis (Hesthammer et al.,, 2000). This leads to preferential flow
directions that are neither parallel nor orthogonal to the main plane. This
type of flow cannot be represented by existing models that employ dimension
reduction. To the Authors’ best knowledge, the only attempt to modeling faults
and their surrounding damage zones in a mixed-dimensional framework can be
found in Fumagalli and Scotti, (2019). However, they still apply local
formulations to model the damage zones as lower-dimensional objects which are
connected on one side to the fault and on the other side to the rock matrix,
hence conceptually seeing the whole fault zone as a multilayer object. An
alternative approach would be to implement the fault core as a
transmissibility multiplier and the damage zone by modifying the grid
permeability in the cells adjacent to the model faults, as illustrated in
Wilson et al., (2020). In the following, we will consistently refer to the
thin inclusions as faults, notwithstanding that all methods presented herein
can be applied to models of fractures and other thin inclusions, however, we
expect that the methods proposed are of more importance for faults.

The contribution of this paper is two-fold: First, we present a generalized
dimensional reduced model that can preserve full-tensor permeability effects
also in the out-of-plane direction of the fault. The resulting reduced
equations have a form similar to that of traditional models, however the more
general coupling structure leads to additional terms both in the in-plane and
out-of-plane equations. These terms, as well as our whole novel formulation,
will be denoted as ”semi-local” in the following, emphasizing the fact that
the new PDEs will contain quantities that, while physically in the same
location, from a modeling perspective reside outside the subdomain where the
PDE is defined, specifically the internal boundary between the subdomain and
its higher dimensional neighbor.

Multiple discretization schemes have been proposed for the local
dimensionally-reduced models, including methods based on finite volumes
(Helmig et al.,, 1997; Karimi-Fard et al.,, 2003; Sandve et al.,, 2012), mixed
finite elements (Martin et al.,, 2005; Boon et al.,, 2018), virtual elements
(Fumagalli and Keilegavlen,, 2019) and mimetic methods (Formaggia et al.,,
2018). A comparison study for all these discretizations of flow in fractured
media can be found in Flemisch et al., (2018) and Berre et al., (2020) for 2D
and 3D flow, respectively. However, the additional terms arising in our
formulation bring the semi-local model outside the scope of previously
proposed discretization methods. The second contribution of the paper is
therefore the derivation of discretization schemes for semi-local models. We
achieve this in two stages: First, based on the unified framework for
discretization of mixed-dimensional problems with local interface laws
presented in Nordbotten et al., (2019), we present conditions under which any
standard discretization scheme for fixed-dimensional problems can be extended
to mixed-dimensional problems with semi-local interface laws. Second, we
present a concrete discretization approach based on finite volume methods.

The paper is organized as follows. In Sec. 2, the mathematical model is
presented, first for a domain with a single fault, and then for a general
faults configuration. Thereafter, in Sec. 3, the unified discretization is
formulated. After presenting simulation results in Sec. 4, concluding remarks
are given in Sec. 5.

##  2 Flow modelling in faulted porous media

In this section, the mathematical model for flow in faulted porous media is
presented, first for a porous domain containing a single fault (Sections 2.1
and 2.2), and then for a general network of faults (Section 2.3). For the
general case, we also provide the weak formulation of the interface problem
(Sections 2.4-2.5), which will be useful from the perspective of
implementation.

\begin{overpic}[width=325.215pt]{upscaling} \put(25.0,20.0){\small$\Psi_{3}$}
\put(30.0,5.0){\small$\Psi_{2}$} \put(30.0,35.0){\small$\Psi_{1}$}
\put(80.0,23.0){\small$\Omega_{3}$} \put(85.0,5.0){\small$\Omega_{2}$}
\put(85.0,35.0){\small$\Omega_{1}$} \put(23.0,37.0){\small${\bm{n}}_{3}$}
\put(23.0,5.0){\small${\bm{n}}_{3}$}
\put(10.0,30.0){\small$\partial_{\Psi_{3}}\Psi_{1}$}
\put(10.0,12.0){\small$\partial_{\Psi_{3}}\Psi_{2}$} \put(7.0,21.0){\small$a$}
\end{overpic} Figure 1: Representation of the fault as a thin three-
dimensional domain $\Psi_{3}$ (left) and as a two-dimensionl manifold
$\Omega_{3}$ (right). The boundary of $\Psi_{j}$ adjacent to $\Psi_{3}$ is
denoted by $\partial_{\Psi_{3}}\Psi_{j}$, for $j=1,2$, while ${\bm{n}}_{i}$ is
the normal vector which is always pointing outwards from $\Psi_{i}$, for
$i=1,2,3$.

###  2.1 Domain with a single fault

We start by considering two three-dimensional porous media $\Psi_{1}$ and
$\Psi_{2}$, each of them with its Neumann and Dirichlet boundaries
$\partial_{N}$ and $\partial_{D}$, respectively. The two three-dimensional
domains are separated by a fault $\Psi_{3}$, which is a thin, almost two-
dimensional object of thickness $a$ (in the following $a$ will be denoted as
the aperture), but which is currently represented as three-dimensional. We
note that $\Psi_{3}$ need not be planar, i.e. $a$ need not be constant. We
denote by $\partial_{\Psi_{3}}\Psi_{j}$, for $j=1,2$, the boundary of
$\Psi_{j}$ adjacent to $\Psi_{3}$. Furthermore, let ${\bm{n}}_{i}$ be the
normal vector which is always pointing outwards from $\Psi_{i}$. It thus
follows that ${\bm{n}}_{3}=-{\bm{n}}_{j}$ on $\partial_{\Psi_{3}}\Psi_{j}$. A
representation of the fault as a thin three-dimensional domain $\Psi_{3}$ is
illustrated in the left of Fig. 1. Darcy flow in the three-dimensional porous
medium is then governed by the following set of equations ($i=1,2,3$):

$\displaystyle\nabla\cdot{\bm{q}}_{i}+f_{i}=0\quad$ $\displaystyle
on\quad\Psi_{i}$ (1) $\displaystyle{\bm{q}}_{i}=-{\bm{K}}_{i}\nabla
p_{i}\quad$ $\displaystyle on\quad\Psi_{i}$ (2)
$\displaystyle\lambda_{3,j}={\bm{q}}_{3}\cdot{\bm{n}}_{3}=-{\bm{q}}_{j}\cdot{\bm{n}}_{j}=-\lambda_{j,3}\quad\quad(j=1,2)\quad$
$\displaystyle on\quad\partial_{\Psi_{3}}\Psi_{j}$ (3)
$\displaystyle{\bm{q}}_{i}\cdot{\bm{n}}_{i}=g_{i}\quad$ $\displaystyle
on\quad\partial_{N}\Psi_{i}$ (4) $\displaystyle\text{tr }p_{i}=0\quad$
$\displaystyle on\quad\partial_{D}\Psi_{i}$ (5)

Here, $p$ is pressure, $\mathbf{q}$ is the Darcy flux, $f$ is a source, and
$\bm{K}$ is a second-order tensor representing the absolute permeability
divided by fluid viscosity. Equation (1) represents mass conservation, while
equation (2) is Darcy’s law. Equation (3) represents flux continuity
conditions on $\partial_{\Psi_{3}}\Psi_{j}$, where $\lambda_{3,j}$ represents
flow from $\Psi_{3}$ to $\Psi_{j}$, thus by flux continuity it follows that
$\lambda_{3,j}=-\lambda_{j,3}$. Finally, equations (4)-(5) are boundary
conditions on $\partial_{N}\Psi_{i}$ and $\partial_{D}\Psi_{i}$, repectively.

Figure 2: Illustration of possible structures of the permeability of a fault
embedded in a porous domain indicated by the principal axis of the
permeability tensor: (a) orthogonal permeability, (b) homogeneous full-
permeability structure, (c) different structure on each half of the fault.

Before deriving the governing equations for the lower-dimensional manifold, we
discuss the decomposition of the permeability tensor within the fault.
Existing local laws for faults as embedded thin inclusions assume that the
principal directions of the local permeability tensor are aligned with the
fault orientation, as illustrated in Fig. 2.a. Hence, more general
orientations of the principal permeability directions, shown in Fig. 2.b-2.c,
cannot be represented by existing models. To be concrete, we let the
permeability on $\Psi_{3}$ have the following decomposition in terms of a
coordinate system aligned with the fault orientation:

${\bm{K}}_{3}=\begin{bmatrix}{\bm{K}}_{3,\parallel}&{\bm{k}}_{3,t}\\\
{\bm{k}}^{T}_{3,t}&k_{3,\bot}\end{bmatrix}$ (6)

Here, ${\bm{K}}_{3,\parallel}$ is a $2\times 2$ second-order tensor
representing the within-fault permeability and $k_{3,\bot}$ is a scalar
representing the normal permeability. The off-diagonal term ${\bm{k}}_{3,t}$
is a two-vector representing the symmetric off-diagonal components of
${\bm{K}}_{3}$; for local interface laws, these off-diagonal components are
assumed to be negligible, i.e. ${\bm{k}}_{3,t}=0$ (Nordbotten and Celia,,
2011; Berre et al.,, 2020). The inclusion of this anisotropic term leads to
significant complications in the modeling and discretization, and is the main
topic of this work. With this structure of the fault permeability, the Darcy
flux for the fault can be decomposed as
${\bm{q}}_{3}=[{\bm{q}}_{3,\parallel},q_{3,\bot}]$, where the 2-vector
tangential component ${\bm{q}}_{3,\parallel}$ and the scalar normal component
$q_{3,\bot}$ have the following form:

$\displaystyle{\bm{q}}_{3,\parallel}=-{\bm{K}}_{3,\parallel}\nabla_{\parallel}p_{3}-{\bm{k}}_{3,t}\nabla_{\bot}p_{3},$
(7) $\displaystyle
q_{3,\bot}=-{\bm{k}}_{3,t}\cdot\nabla_{\parallel}p_{3}-k_{3,\bot}\nabla_{\bot}p_{3}.$
(8)

Here, $\nabla_{\parallel}$ and $\nabla_{\bot}=\dfrac{\partial}{\partial n}$
represent the in-plane and out-of-plane components of the gradient for the
fault, respectively.

\begin{overpic}[width=216.81pt]{sketch_interfaces}
\put(32.0,50.0){\small$\Omega_{1}$} \put(32.0,25.0){\small$\Omega_{2}$}
\put(27.0,40.0){\small$\Omega_{3}$} \put(102.0,37.0){\small$\Omega_{3}$}
\put(102.0,29.0){\small$\Gamma_{2,3}$} \put(102.0,45.0){\small$\Gamma_{1,3}$}
\put(85.0,56.0){\small$\partial_{\Omega_{3}}\Omega_{1}$}
\put(85.0,18.0){\small$\partial_{\Omega_{3}}\Omega_{2}$}
\put(75.0,65.0){\small$\Omega_{1}$} \put(75.0,5.0){\small$\Omega_{2}$}
\end{overpic} Figure 3: Illustration of the mixed-dimensional geometry.
$\Omega_{3}$ is connected to the higher dimensional neighbors $\Omega_{j}$
through the interfaces $\Gamma_{j,3}$, for $j=1,2$. Note that $\Omega_{3}$,
$\Gamma_{j,3}$ and $\partial_{\Omega_{3}}\Omega_{j}$ are all coinciding in
physical space.

###  2.2 Model reduction

To proceed, we apply integration over the perpendicular direction to achieve a
dimension reduction of the fault. This replaces $\Psi_{3}$ with a lower-
dimensional domain $\Omega_{3}$ (see right of Fig. 1). Note that we use $\Psi$
to represent the equi-dimensional geometry, that is all $\Psi_{j}$ are 3D, and
$\Omega$ to denote the mixed-dimensional geometry. We also introduce two
interfaces $\Gamma_{j,3}$ on each side $j=1,2$ of $\Omega_{3}$, as illustrated
in Fig. 3. The interfaces physically represent the half zone comprised between
the fault and either side of the surrounding matrix. In a mixed-dimensional
setting, they have no perpendicular extent, and serve as connectors between
two objects of different dimensions. Note that, due to the dimension reduction
of the model, $\Omega_{3}$, $\Gamma_{1,3}$, $\Gamma_{2,3}$,
$\partial_{\Omega_{3}}\Omega_{1}$ and $\partial_{\Omega_{3}}\Omega_{2}$ are
all coinciding in physical space. Furthermore, we define the integrated Darcy
flux ${\bm{q}}_{3}^{(2)}$ and the average pressure $p^{(2)}_{3}$, respectively
as

${\bm{q}}_{3}^{(2)}=\int_{-a/2}^{a/2}{\bm{q}}^{(3)}_{3,\parallel}dn,\quad\quad
p^{(2)}_{3}=\dfrac{1}{a}\int_{-a/2}^{a/2}p_{3}^{(3)}dn.$ (9)

Here, we use subscripts to index the domains, and superscripts (when necessary
for clarity) to emphasize the effective topological dimension of the domain,
e.g. $p_{3}^{(3)}$ and $p_{3}^{(2)}$ are the pressures within the fault in the
3D (on $\Psi_{3}$) and 2D (on $\Omega_{3}$) representations, respectively.
When passing to a mixed-dimensional representation of the geometry, i.e. when
integrating eqs. (1) and (7) along the perpendicular direction, the out-of-
plane component of the gradient is converted into a jump operator as follows:

$\int_{-a/2}^{a/2}\nabla_{\bot}p_{3}^{(3)}dn=(\text{tr }p_{1}-\text{tr
}p_{2}).$ (10)

The governing equations for the fault are then obtained from equations (1),
(7), (4) and (5) by integrating in the perpendicular direction. Moreover,
since the fault is assumed to be thin, we assume that the permeability is
constant across the perpendicular direction. Together with the definitions
above, this results in

$\displaystyle\nabla_{3}\cdot{\bm{q}}_{3}^{(2)}-(\lambda_{1,3}+\lambda_{2,3})+f_{3}^{(2)}=0\quad$
$\displaystyle on\quad\Omega_{3}$ (11)
$\displaystyle{\bm{q}}_{3}^{(2)}=-a{\bm{K}}_{3,\parallel}\nabla_{3}p_{3}^{(2)}+{\bm{\mu}}_{1,3}+{\bm{\mu}}_{2,3}\quad$
$\displaystyle on\quad\Omega_{3}$ (12)
$\displaystyle{\bm{q}}_{3}^{(2)}\cdot{\bm{n}}_{3}^{(2)}=g_{3}^{(2)}\quad$
$\displaystyle on\quad\partial_{N}\Omega_{3}$ (13) $\displaystyle\text{tr
}p_{3}^{(2)}=0\quad$ $\displaystyle on\quad\partial_{D}\Omega_{3}$ (14)

where we have also introduced the integrated source term and boundary flux

$f_{3}^{(2)}=\int_{-a/2}^{a/2}f_{3}^{(3)}dn,\quad\quad
g_{3}^{(2)}=\int_{-a/2}^{a/2}g_{3}^{(3)}dn.$ (15)

We emphasize that the differential operator $\nabla_{3}$ in eqs. (11)-(12)
operates on the manifold $\Omega_{3}$. Compared to traditional upscaled
models, see for instance Nordbotten et al., (2019), additional terms
${\bm{\mu}}_{j,3}$ appear in equation (12), analogous to the flux terms
$\lambda_{j,3}$ in equation (11). This two-vector term, which is not present
in previous work, represents the within-fault flux induced by pressure
differences between the fault and the surrounding matrix, and is defined for
either side of the fault as

${\bm{\mu}}_{j,3}=\epsilon_{j,3}{\bm{k}}_{3,t}(p_{3}^{(2)}-\text{tr }p_{j}),$
(16)

where the permutation variable $\epsilon_{j,3}$ is positive if the coordinate
systems of $\Omega_{3}$ and $\partial_{\Omega_{3}}\Omega_{j}$ coincide, and
negative otherwise.

To complete the model, we derive a constitutive law for $\lambda_{j,3}$. This
is obtained by integrating equation (8) in the perpendicular direction, that
is

$\int_{-a/2}^{a/2}q_{3,\bot}^{(3)}dn=-\int_{-a/2}^{a/2}{\bm{k}}_{3,t}\cdot\nabla_{\parallel}p_{3}^{(3)}dn-\int_{-a/2}^{a/2}k_{3,\bot}\nabla_{\bot}p_{3}^{(3)}dn.$
(17)

The left hand side of equation (17) is approximated using the trapeizodal
rule, that is

$\int_{-a/2}^{a/2}q_{3,\bot}^{(3)}dn\approx\dfrac{a}{2}(\epsilon_{1,3}\lambda_{1,3}+\epsilon_{2,3}\lambda_{2,3}),$
(18)

where continuity of the flux across the boundary between the fault and the
surrounding matrix is applied. The first term at the right hand side of
equation (17) is approximated as

$\int_{-a/2}^{a/2}{\bm{k}}_{3,t}\cdot\nabla_{\parallel}p_{3}^{(3)}dn={\bm{k}}_{3,t}\cdot\int_{-a/2}^{a/2}\nabla_{\parallel}p_{3}^{(3)}dn\approx{\color[rgb]{1,0,0}a}{\bm{k}}_{3,t}\cdot\nabla_{3}p_{3}^{(2)}.$
(19)

Finally, the second term at the right hand side of (17) is resolved using the
jump operator defined in equation (10) as follows:

$\int_{-a/2}^{a/2}k_{3,\bot}\nabla_{\bot}p_{3}^{(3)}dn=\epsilon_{1,3}k_{3,\bot}(p_{3}^{(2)}-\text{tr
}p_{1})+\epsilon_{2,3}k_{3,\bot}(p_{3}^{(2)}-\text{tr }p_{2}).$ (20)

By incorporating eqs. (18), (19) and (20) into equation (17), we identify the
flux $\lambda_{j,3}$ having the following form:

$\lambda_{j,3}=-k_{3,\bot}\dfrac{2(p_{3}^{(2)}-\text{tr
}p_{j})}{a}-\epsilon_{j,3}{\bm{k}}_{3,t}\cdot\nabla_{3}p_{3}^{(2)}.$ (21)

Here, the first term on the right-hand side represents the local component of
the constitutive law, while the second part is the semi-local contribution
that induces a flux across $\Gamma_{j,3}$ due to the pressure gradient within
the lower-dimensional manifold $\Omega_{3}$.

Inspecting equations (16) and (21), we see that both the normal permeability
$k_{3,\bot}$ and the off-diagonal permeability ${\bm{k}}_{3,t}$ are in the
reduced model naturally interpreted as properties of the interface
$\Gamma_{j,3}$. In the continuation, we will thus generalize the model as
derived above, and index these quantities with the interface, i.e.
${\bm{k}}_{3,t}\rightarrow{\bm{k}}_{3,j,t}$ and $k_{3,\bot}\rightarrow
k_{3,j,\bot}$ are assigned independently to either side of the fault.

In summary, omitting superscripts for the sake of clarity, we can write the
mixed-dimensional equations (1)-(5), (11)-(14), (16), and (21) in a unified
way, that is for $i=\\{1,2,3\\}$

$\displaystyle\nabla_{i}\cdot{\bm{q}}_{i}-\sum_{j\in{\hat{S}}_{i}}\lambda_{j,i}+f_{i}=0\quad$
$\displaystyle on\quad\Omega_{i}$ (22)
$\displaystyle{\bm{q}}_{i}=-{\bm{\kappa}}_{i,\parallel}\nabla_{i}p_{i}+\sum_{j\in{\hat{S}}_{i}}\epsilon_{j,i}{\bm{\kappa}}_{i,j,t}(p_{i}-\text{tr
}p_{j})\quad$ $\displaystyle on\quad\Omega_{i}$ (23)
$\displaystyle{\bm{q}}_{i}\cdot{\bm{n}}_{i}=\lambda_{i,3}\quad\quad(i\neq
3)\quad$ $\displaystyle on\quad\partial_{\Omega_{3}}\Omega_{i}$ (24)
$\displaystyle\lambda_{j,3}=-\kappa_{3,j,\bot}(p_{3}-\text{tr
}p_{j})-{\bm{\kappa}}_{3,j,t}\cdot\nabla_{3}p_{3}\quad\quad(j=1,2)\quad$
$\displaystyle on\quad\Gamma_{j,3}$ (25)
$\displaystyle{\bm{q}}_{i}\cdot{\bm{n}}_{i}=g_{i}\quad$ $\displaystyle
on\quad\partial_{N}\Omega_{i}$ (26) $\displaystyle\text{tr }p_{i}=0\quad$
$\displaystyle on\quad\partial_{D}\Omega_{i}$ (27)

where ${\hat{S}}_{i}$ is the set of neighbors of $\Omega_{i}$ of higher
dimension, e.g. ${\hat{S}}_{3}=\\{\Omega_{1},\Omega_{2}\\}$. Equations
(22)-(27) are complemented with the natural convention that there is no four-
dimensional domain in the model, thus ${\hat{S}}_{i}=\emptyset$ for $i=1,2$,
and one clearly has for these three-dimensional domains also that $a=1$,
${\bm{K}}_{\parallel}=\bm{K}$ and $\nabla_{i}=\nabla$.

We remark that due to the model reduction, the within-fault permeability
${\bm{K}}_{3,\parallel}$ and the normal permeability $k_{3,j,\bot}$ scale with
the aperture $a$ and its inverse, respectively, while the off-diagonal
permeability ${\bm{k}}_{3,j,t}$ remains as in the equi-dimensional model. In
order to present equations (22)-(27) without reference to this small
parameter, these scalings have been incorporated directly into the material
constants. Thus, the mixed-dimensional permeability ${\bm{\kappa}}_{3}$ is
related to the equi-dimensional ${\bm{K}}_{3}$ as follows

${\bm{\kappa}}_{3}=\begin{bmatrix}{\bm{\kappa}}_{3,\parallel}&{\bm{\kappa}}_{3,j,t}\\\
{\bm{\kappa}}^{T}_{3,j,t}&\kappa_{3,j,\bot}\end{bmatrix}=\begin{bmatrix}a{\bm{K}}_{3,\parallel}&{\bm{k}}_{3,j,t}\\\
{\bm{k}}^{T}_{3,j,t}&2a^{-1}k_{3,j,\bot}\end{bmatrix}.$ (28)

We point out that, when one reduces multiple dimensions at once, these
scalings get exponents corresponding to the number of dimensions below the
ambient dimension. We also emphasize that the normal and off-diagonal
permeabilities are in principle not a property of the fault itself, but
instead a property which belongs to the internal interface $\Gamma_{j,3}$
between the fault and either side of the higher-dimensional neighbors. This
represents an important extension of the existing local laws for fractured
porous media, making the model also applicable to faulted porous media, since
it allows for capturing the anisotropic character of the fault damage zone.
Moreover, since different values of ${\bm{k}}_{3,j,t}$ and $k_{3,j,\bot}$ can
be assigned to each side of the fault, our model can represent different
permeability structures on each side of the fault.

A schematic illustration of the different quantities and their domain of
definition for the local and semi-local formulations is shown in Fig. 4.

\begin{overpic}[width=216.81pt]{local_vs_semi_local} \put(-35.0,70.0){Semi-
local} \put(110.0,70.0){Local} \put(50.0,15.0){\small${\bm{q}}_{j},p_{j}$}
\put(-45.0,55.0){\small$\lambda_{j,3}\sim(\nabla_{3}p_{3},p_{j},p_{3})$}
\put(-45.0,17.0){\small${\bm{q}}_{3}\sim(\nabla_{3}p_{3},p_{3},p_{1},p_{2})$}
\put(103.0,60.0){\small$\lambda_{j,3}\sim(p_{j},p_{3})$}
\put(105.0,32.0){\small${\bm{q}}_{3}\sim\nabla_{3}p_{3}$} \end{overpic} Figure
4: Illustration of the quantities associated with the local and semi-local
formulations.

###  2.3 Domain with a general network of faults

Following the theory by Boon et al., (2018), equations (22)-(27) can be
generalized also to faults intersections, both the one-dimensional (1D) line
intersections between two faults and the zero-dimensional (0D) point
intersections of three faults (see Fig. 5 for an illustration of the mixed-
dimensional geometry). To this end, we use subscripts to index each domain
(matrix, fault, or intersection) by number as in the previous section, and let
$I$ denote the index set of all domains. Superscripts for the topological
dimension associated with each individual domain will be consistently omitted,
keeping in mind that the dimension is always a property of the domain, i.e.
$d=d_{i}$. Hence, we can write for all $i\in I$ the equations

$\displaystyle\nabla_{i}\cdot{\bm{q}}_{i}-\sum_{j\in{\hat{S}}_{i}}\lambda_{j,i}+f_{i}=0\quad$
$\displaystyle on\quad\Omega_{i}$ (29)
$\displaystyle{\bm{q}}_{i}=-{\bm{\kappa}}_{i,\parallel}\nabla_{i}p_{i}+\sum_{j\in{\hat{S}}_{i}}\epsilon_{j,i}{\bm{\kappa}}_{i,j,t}(p_{i}-\text{tr
}p_{j})\quad$ $\displaystyle on\quad\Omega_{i}$ (30)
$\displaystyle{\bm{q}}_{i}\cdot{\bm{n}}_{i}=\lambda_{i,j}\quad\quad(j\in{\check{S}}_{i})\quad$
$\displaystyle on\quad\partial_{\Omega_{j}}\Omega_{i}$ (31)
$\displaystyle\lambda_{j,i}=-\kappa_{i,j,\bot}(p_{i}-\text{tr
}p_{j})-\epsilon_{j,i}{\bm{\kappa}}_{i,j,t}\cdot\nabla_{i}p_{i}\quad\quad(j\in{\hat{S}}_{i})\quad$
$\displaystyle on\quad\Gamma_{j,i}$ (32)
$\displaystyle{\bm{q}}_{i}\cdot{\bm{n}}_{i}=g_{i}\quad$ $\displaystyle
on\quad\partial_{N}\Omega_{i}$ (33) $\displaystyle\text{tr }p_{i}=0\quad$
$\displaystyle on\quad\partial_{D}\Omega_{i}$ (34)

where ${\check{S}}_{i}$ is the set of neighbors of $\Omega_{i}$ of lower
dimension, e.g. ${\check{S}}_{1}=\\{\Omega_{2},\Omega_{3},\Omega_{4}\\}$. It
is easy to show that as long as the mixed-dimensional permeabilities are
diagonally dominant in the sense of

$\kappa_{i,j,\bot}\det{\bm{\kappa}}_{i,\parallel}>{\bm{\kappa}}_{i,j,t}\cdot{\bm{\kappa}}_{i,j,t},$
(35)

then the coefficients are globally positive definite, and equations (29)-(34)
are well-posed as long as $\partial_{D}\Omega_{i}$ has non-zero measure for at
least one domain (Boon et al.,, 2020).

\begin{overpic}[width=173.44534pt]{sketch_domain_network_inclusions_separate}
\put(52.0,95.0){\small$\Omega_{1}$} \put(5.0,52.0){\small$\Omega_{2}$}
\put(40.0,47.0){\small$\Omega_{3}$} \put(60.0,45.0){\small$\Omega_{4}$}
\put(25.0,20.0){\small$\Omega_{5}$} \put(50.0,26.0){\small$\Omega_{7}$}
\put(36.0,30.0){\small$\Omega_{6}$} \put(36.0,2.0){\small$\Omega_{8}$}
\end{overpic} Figure 5: Illustration of a 3D domain $\Omega_{1}$ containing
three faults $\Omega_{j}$ ($j=2,3,4$) with their respective three 1D line
intersections $\Omega_{k}$ ($k=5,6,7$) and one 0D point intersection
$\Omega_{8}$.

###  2.4 Mixed-dimensional formulation of the fault-matrix flows

While equations (29)-(34) constitute a full semi-local model, they are stated
in a form which is not immediately amenable for discretization. This and the
following subsection explore the model in more detail, with a goal of
rewriting the equations in a form that can be handled by standard
discretization schemes with only minimal adaptations. A discretization
approach based on this reformulation is then given in Section 3.

In order to simplify the exposition, we will introduce a mixed-dimensional
notation following Nordbotten et al., (2019). In particular, we will denote
the collection of pressure functions as
$\mathfrak{p}=\left(p_{1},...,p_{|I|}\right)$, and similarly the collection of
all fluxes (both in domains and across boundaries) as
$\mathfrak{q}=\left({\bm{q}}_{1},...,{\bm{q}}_{|I|},\lambda_{1,1},...,\lambda_{|I|,|J|}\right)$.
It is sometimes convenient to refer explicitly to only the domain or boundary
fluxes, and we will therefore sometimes abuse notation and simply write
$\mathfrak{q}=(q,\lambda)$. We refer to these as mixed-dimensional functions,
and consistently denote them with calligraphic font. We adopt the natural
convention that when evaluating a mixed-dimensional function at a point, say
$x\in\Omega_{i}$, then we simply evaluate the function on that domain, so that
$\mathfrak{p}(x)=p_{i}(x)$. In a similar sense, we denote the disjoint union
of domains as
$\mathfrak{F}=\left(\amalg_{i}\Omega_{i}\right)\sqcup\left(\amalg_{j,i}\Gamma_{j,i}\right)$.

With this notion of mixed-dimensional functions, the extension of the
divergence and gradient operators to the mixed-dimensional setting is natural.
First, we extend the concept of continuous functions by requiring that for
$\mathfrak{q}$ to be continuous, then it must hold that, for all
$\Gamma_{j,i}$, ${\bm{q}}_{i}\cdot{\bm{n}}_{i}=\lambda_{i,j}$. Then, for any
point $x\in\Omega_{i}$ we define

$\left(\mathfrak{D}\cdot\mathfrak{q}\right)(x)=\left[\nabla_{i}\cdot{\bm{q}}_{i}-\sum_{j\in{\hat{S}}_{i}}\lambda_{j,i}\right]_{x}\quad\quad\text{and}\quad\quad\left(\mathbb{D}\mathfrak{p}\right)(x)=\left[\nabla_{i}p_{i}\right]_{x},$
(36)

while for any point on an interface $x\in\Gamma_{j,i}$ we define

$\left(\mathbb{D}\mathfrak{p}\right)(x)=\left[p_{i}-\text{tr
}p_{j}\right]_{x}.$ (37)

Now we can write equations (29) - (34) simply as:

$\displaystyle\mathfrak{D}\cdot\mathfrak{q}+\mathfrak{f}=0\quad$
$\displaystyle on\quad\mathfrak{F}$ (38)
$\displaystyle\mathfrak{q}=-\mathfrak{K}\mathbb{D}p\quad$ $\displaystyle
on\quad\mathfrak{F}$ (39)
$\displaystyle\mathfrak{q}\cdot\mathfrak{n}=\mathfrak{g}\quad$ $\displaystyle
on\quad\partial_{N}\mathfrak{F}$ (40) $\displaystyle\text{tr
}\mathfrak{p}=0\quad$ $\displaystyle on\quad\partial_{D}\mathfrak{F}$ (41)

where we have also introduced the collection of sources
$\mathfrak{f}=\left(f_{1},...,f_{|I|}\right)$ and the collection of boundary
fluxes $\mathfrak{g}=\left(g_{1},...,g_{|I|}\right)$. Here, the material
coefficients are now all part of the mixed-dimensional permeability
$\mathfrak{K}$, which is defined such as that for any mixed-dimensional
gradient $\mathfrak{u}=\mathbb{D}p=(u,\mu)$, it holds that for any point
$x\in\Omega_{i}$:

$\left(\mathfrak{K}\mathfrak{u}\right)(x)={\bm{\kappa}}_{i,\parallel}u_{i}-\sum_{j\in{\hat{S}}_{i}}\epsilon_{j,i}{\bm{\kappa}}_{i,j,t}\mu_{j,i},$
(42)

while for any point on an interface $x\in\Gamma_{j,i}$, it holds that

$\left(\mathfrak{K}\mathfrak{u}\right)(x)=\kappa_{i,j,\bot}\mu_{j,i}+\epsilon_{j,i}{\bm{\kappa}}_{i,j,t}\cdot
u_{i}.$ (43)

It is then also sometimes convenient to write equation (39) in matrix form,
that is for $\mathfrak{q}=(q,\lambda)$ and
$\mathfrak{u}=\mathbb{D}\mathfrak{p}=(u,\mu)$, one has:

$\begin{Bmatrix}q\\\
\lambda\end{Bmatrix}=-\begin{bmatrix}\mathfrak{K}_{\Omega\Omega}&\mathfrak{K}_{\Omega\Gamma}\\\
\mathfrak{K}_{\Gamma\Omega}&\mathfrak{K}_{\Gamma\Gamma}\end{bmatrix}\begin{Bmatrix}u\\\
\mu\end{Bmatrix}.$ (44)

Equation (44) highlights the contribution from the semi-local terms in the
mixed-dimensional version of Darcy’s law.

###  2.5 Weak formulation as an interface system

The semi-local terms in equations (29)-(34) lead to coupling terms between
domains that are local in physical space, but non-local in the mixed-
dimensional representation of the geometry. A critical example are the fault
and its sides, which, from the perspective of implementation, we would prefer
to only interact via the interfaces $\Gamma_{j,i}$, and not directly, as is
the case for the last term in equation (30).

Thus we are motivated to consider a reformulation of the governing equations
before considering numerical discretizations. We proceed by first performing
an LU decomposition of equation (44) as follows:

$\mathfrak{K}_{U}\begin{Bmatrix}q\\\
\lambda\end{Bmatrix}=-\mathfrak{K}_{L}\begin{Bmatrix}u\\\ \mu\end{Bmatrix},$
(45)

where $\mathfrak{K}_{U}$ and $\mathfrak{K}_{L}$ are defined, respectively, as:

$\mathfrak{K}_{U}=\begin{bmatrix}I&\mathfrak{K}_{\Omega\Gamma}\mathfrak{K}_{\Gamma\Gamma}^{-1}\\\
0&I\end{bmatrix}\quad\quad\text{and}\quad\quad\mathfrak{K}_{L}=\begin{bmatrix}A_{\Omega}&0\\\
\mathfrak{K}_{\Gamma\Omega}&\mathfrak{K}_{\Gamma\Gamma}\end{bmatrix},$ (46)

and $A_{\Omega}$ is the Schur-complement defined as

$A_{\Omega}=\mathfrak{K}_{\Omega\Omega}-\mathfrak{K}_{\Omega\Gamma}\mathfrak{K}_{\Gamma\Gamma}^{-1}\mathfrak{K}_{\Gamma\Omega}.$
(47)

Note that, since $\mathfrak{K}_{\Gamma\Gamma}$ consists only of scalar values
$(\kappa_{i,j,\bot})$, this reformulation only depends on the trivial
inversion of scalars.

In the following it will be helpful to discuss the components of the mixed-
dimensional gradient and divergence, and we therefore additionally define the
”full jump” $\mathbb{d}\mathfrak{q}$ such that for any point $x\in\Omega_{i}$
it holds that

$\left(\mathbb{d}\mathfrak{q}\right)(x)=\left[-\sum_{j\in{\hat{S}}_{i}}\lambda_{j,i}\right]_{x},$
(48)

while the ”half jump” $\mathbb{d}^{\star}\mathfrak{p}$ is simply the
restriction of $\mathbb{D}\mathfrak{p}$ to $\Gamma_{j,i}$. We then write (with
the natural extension of $\nabla$ and $\nabla\cdot$):

$\mathfrak{D}\cdot\mathfrak{q}=\nabla\cdot
q+\mathbb{d}\lambda\quad\quad\text{and}\quad\quad\mathbb{D}\mathfrak{p}=\left(\nabla
p,\mathbb{d}^{\star}\mathfrak{p}\right).$ (49)

We now proceed by (formally) eliminating internal domain variables, in order
to obtain a problem only posed on interfaces. We note that equations (38) and
(39) can now be written as the first order system:

$\displaystyle\mathfrak{D}\cdot\mathfrak{q}=\mathfrak{f}$ (50)
$\displaystyle\mathfrak{K}_{U}\mathfrak{q}=-\mathfrak{K}_{L}\mathfrak{D}\mathfrak{p}$
(51)

where use of equation (45) has been made. By writing out equation (50) in
local notation for each $\Omega_{i}$ and by stating equation (51) explicitly
as two equations, we obtain the following set of equations:

$\displaystyle\nabla\cdot q=f-\mathbb{d}\lambda$ (52) $\displaystyle
q+A_{\Omega}\nabla
p=-\mathfrak{K}_{\Omega\Gamma}\mathfrak{K}_{\Gamma\Gamma}^{-1}\lambda$ (53)
$\displaystyle\lambda=-\left(\mathfrak{K}_{\Gamma\Omega}\nabla
p+\mathfrak{K}_{\Gamma\Gamma}\mathbb{d}^{\star}\mathfrak{p}\right)$ (54)

This reveals that equations (52) and (53) form a locally well-posed system (of
standard Darcy type) on each $\Omega_{i}$, and we can therefore consider
$p=p(\lambda)$ for any given $\lambda$.

We formalize this concept by introducing the (continuous) solution operators
for the standard elliptic value problem on $\Omega_{i}$,
$\mathcal{S}_{\Omega_{i}}^{K}$, defined as:

$\left(\upsilon,\nabla\upsilon,\text{tr
}\upsilon,F\right)=\mathcal{S}_{\Omega_{i}}^{K}\left(f,\chi,b,\upsilon_{0}\right),$
(55)

where $\upsilon$ is the solution to

$\displaystyle\nabla\cdot\varphi=f-F\quad\quad$
$\displaystyle\text{on}\quad\Omega_{i}$ (56)
$\displaystyle\varphi=-K\left(\nabla\upsilon+\chi\right)\quad\quad$
$\displaystyle\text{on}\quad\Omega_{i}$ (57) $\displaystyle\varphi\cdot
n=b\quad\quad$
$\displaystyle\text{on}\quad\partial\Omega_{i}\setminus\partial\Omega$ (58)
$\displaystyle\upsilon=0\quad\quad$
$\displaystyle\text{on}\quad\partial\Omega_{i}\cap\partial\Omega$ (59)
$\displaystyle\dfrac{1}{|\Omega_{i}|}\int_{\Omega_{i}}\upsilon=\upsilon_{0}\quad\quad$
$\displaystyle\text{if}\quad\partial\Omega_{i}\cap\partial\Omega\neq\oslash$
(60)

where $\partial\Omega$ is the global boundary and
$F=\dfrac{1}{|\Omega_{i}|}\left(\int_{\Omega_{i}}f-\int_{\Omega_{i}}b\right)$
if $\partial\Omega_{i}\cap\partial\Omega\neq\oslash$, and zero otherwise.
Using this solution operator, we see that the solution to equations (52) and
(53) can be stated as functions of $\lambda$ (and a set of number of numbers
$p_{0}$ corresponding to the domains where
$\partial\Omega_{i}\cap\partial\Omega\neq\oslash$) as:

$\left(p,\nabla p,\text{tr
}p,F\right)_{\Omega_{i}}\left(\lambda,p_{0}\right)=\mathcal{S}_{\Omega_{i}}^{A_{i}}\left(f_{i}-\left(\mathbb{d}\lambda\right)_{i},A_{i}^{-1}\left(\mathfrak{K}_{\Omega\Gamma}\mathfrak{K}_{\Gamma\Gamma}^{-1}\lambda\right)_{i},\lambda_{\check{I}_{i}},p_{0}\right).$
(61)

Inserting $p=p(\lambda,p_{0})$ etc. into equation (54), we have now
reformulated the fault-matrix problem into a pure interface problem. From the
perspective of implementation, we desire to consider the interface problem in
the weak sense, and we therefore multiply by test functions $w$ and integrate
to obtain the problem: Find $\lambda\in L^{2}(\Gamma)$ such that, for all
$w\in L^{2}(\Gamma_{j})$

$\left(\mathfrak{K}_{\Gamma\Gamma}^{-1}\lambda,w\right)_{\Gamma_{j,i}}+\left(\mathfrak{K}_{\Gamma\Gamma}^{-1}\mathfrak{K}_{\Gamma\Omega}\nabla
p(\lambda,p_{0}),w\right)_{\Gamma_{j,i}}+\left(\mathbb{d}^{\star}\mathfrak{p}(\lambda,p_{0}),w\right)_{\Gamma_{j,i}}=0$
(62)

and $F_{i}(\lambda,p_{0})=0$ if
$\partial\Omega_{i}\cap\partial\Omega\neq\oslash$. We point out that the inner
products in equation (62) are bounded from a formal perspective, since for
$\lambda\in L^{2}(\Gamma)$, then $p_{i}\in H^{1}(\Omega_{i})$, and both
$\mathfrak{K}_{\Gamma\Gamma}^{-1}\mathfrak{K}_{\Gamma\Omega}\nabla p$ and
$\text{tr }p$ will lie in (at least) $L^{2}(\Gamma_{j,i})$.

Finally, we emphasize that equations (61)-(62) are attractive from the
perspective of implementation, since the inner products appearing are easy to
evaluate, and the solution operators $\mathcal{S}_{\Omega_{i}}^{A_{i}}$ can be
approximated by any standard method, as we will detail in the next section.

##  3 Discretizations of flow for faulted porous media

The equations derived in Section 2.5, and in particular the interface problem
of equation (61), form the starting point for the discretization approach laid
out in this section. We present the general discretization framework in
Section 3.1, and discuss implementational aspects in Section 3.2.

###  3.1 Unified discretization

Equation (61) provides a solution operator for the arbitrary standard method
used to solve the elliptic boundary value problem (52)-(53) on $\Omega_{i}$.
To be concrete, we consider each domain $\Omega_{i}$ and its Neumann boundary
$\partial\Omega_{i}=\partial_{N}\Omega_{i}\cup\sum_{j\in\check{S}_{i}}\partial_{\Omega_{j}}\Omega_{i}$
as endowed with a numerical discretization. Then, the solution operator
$\mathcal{S}_{i}$ can be stated as

$\mathcal{S}_{i}:\left[N(\Omega_{i}),N^{d_{i}}(\Omega_{i}),N(\partial\Omega_{i})]\rightarrow[N(\Omega_{i}),N^{d_{i}}(\Omega_{i}),N(\partial\Omega_{i})\right],$
(63)

where $N(\Omega_{i})$, $N^{d_{i}}(\Omega_{i})$, and $N(\partial\Omega_{i})$
are the discrete representations of $L^{2}(\Omega_{i})$,
$\left(L^{2}(\Omega_{i})\right)^{d_{i}}$, and $L^{2}(\partial\Omega_{i})$,
respectively, and ${d_{i}}$ is the topological dimension of $\Omega_{i}$. In
particular, $\mathcal{S}_{i}$ takes as input sinks, vector sources, and
Neumann data and returns as output pressures, pressure gradients, and pressure
traces. Most discretization schemes for elliptic equations can provide such a
solution operator; we discuss the concrete implementation in the next
subsection.

To discretize the flux coupling term $\lambda_{j,i}$, we introduce a mortar-
like grid $\mathcal{T}_{j,i}$ on the interface $\Gamma_{j,i}$ on which the
boundary flux $\lambda_{j,i}$ will be defined. The flux variables are
represented as piecewise constant on the mortar grid $\mathcal{T}_{j,i}$, thus
$\lambda_{j,i}\in P_{0}(\mathcal{T}_{j,i})\subset\L^{2}(\Omega_{i})$. In order
to allow communications between subdomains, and thus explicitly relate the
degrees of freedom of the numerical methods $\mathcal{S}_{i}$ and the mortar
grids $\mathcal{T}_{j,i}$, we introduce projection operators, namely
$\Pi_{N(\Omega_{i})}$ and $\Pi_{L^{2}(\Omega_{i})}$. The former is the
compound operator projecting from the coupling variables on the mortar grids
to the subdomain degrees of freedom, that is

$\begin{split}\Pi_{N(\Omega_{i})}:&\left[L^{2}(\Omega_{i}),\left(L^{2}(\Omega_{i})\right)^{d_{i}},L^{2}\left(\Omega_{\check{S}_{i}}\right),L^{2}(\partial\Omega_{i})\right]\\\
&\rightarrow\left[N(\Omega_{i}),N^{d_{i}}(\Omega_{i}),N(\partial\Omega_{i})\right],\end{split}$
(64)

while the latter conversely moves from the numerical variables to the coupling
variables, that is

$\begin{split}\Pi_{L^{2}(\Omega_{i})}:&\left[N(\Omega_{i}),N^{d_{i}}(\Omega_{i}),N(\partial\Omega_{i})\right]\\\
&\rightarrow\left[L^{2}(\Omega_{i}),\left(L^{2}(\Omega_{i})\right)^{d_{i}},L^{2}\left(\Omega_{\check{S}_{i}}\right),L^{2}(\partial\Omega_{i})\right].\end{split}$
(65)

Now, following the variational formulation derived in Sec. 2.5, we exploit
equation (62) in order to provide discretization-independent framework for
faulted porous media. This takes the form: for given numerical discretizations
$\mathcal{S}_{i}$, find $\lambda_{j,i}\in P_{0}(\mathcal{T}_{j,i})$, for all
$i\in I$ and $j\in\hat{S}_{i}$ such that

$\begin{split}&\left(\mathbb{d}^{\star}\mathfrak{p},w_{j}\right)_{\Gamma_{j,i}}+\left(\mathfrak{K}^{-1}_{\Gamma\Gamma}\left(\lambda_{j,i}+\mathfrak{K}_{\Gamma\Omega}\cdot\nabla
p\right),w_{j}\right)_{\Gamma_{j,i}}=0\\\ &\quad\text{for all }w_{j}\in
P_{0}(\mathcal{T}_{j,i})\\\ \end{split}$ (66)

subject to discrete constraints (for all $i\in I$):

$\displaystyle[p_{i},u_{i},t_{j}]$
$\displaystyle=\Pi_{L^{2}\left(\Omega_{i}\right)}\mathcal{S}_{i}(\psi_{i}+a_{i},b_{i},c_{i})$
(67) $\displaystyle[a_{i},b_{i},c_{i}]$
$\displaystyle=\Pi_{N(\Omega_{i})}\left[-\sum_{j\in\hat{S}_{i}}\lambda_{j,i},-\sum_{j\in\hat{S}_{i}}A_{i}^{-1}\mathfrak{K}_{\Omega\Gamma}\mathfrak{K}_{\Gamma\Gamma}^{-1}\lambda_{j,i},\sum_{j\in\check{S}_{i}}\lambda_{i,j}\right]$
(68)

where the dummy variables $a_{i}$, $b_{i}$ and $c_{i}$ have the
interpretations of sources, forces, and fluxes due to interactions with other
domains, respectively. In contrast, the variables $p_{i}$, $u_{i}$, and
$t_{j}$ are the pressures, pressure gradients, and pressure traces after
projection onto the grids $\mathcal{T}_{j,i}$.

The interpretation of this scheme is as follows. Eq. (67) resolves the
internal differential equations in each subdomain, eq.(68) is the projection
of variables from the flux grids to the numerical boundary (and source) data,
while equation (66) simply states that the flux $\lambda_{j,i}$ between the
fault and its surrodundings should satisfy Darcy’s law. In the following
section, we present the strategy for implementation of this approach and give
details for a specific numerical scheme.

###  3.2 MPFA discretization

It is of interest to consider the requirements put on the subdomain solution
operators $\mathcal{S}_{i}$ in some more detail. From the variational
formulations stated above, we see that for a discretization on a generic
subdomain $\Omega_{i}$ to interact with the interface $\Gamma_{j}$, we need to
provide operators which:

  1. 1.

Handle Neumann boundary data on the form $\Pi_{N(\Omega_{i})}\lambda_{j}$, for
all interfaces $\Gamma_{j}$ where $\Omega_{i}$ is the higher-dimensional
neighbor.

  2. 2.

Handle source terms $\Pi_{N(\Omega_{i})}\lambda_{j}$ from interfaces
$\Gamma_{j}$ where $\Omega_{i}$ is the lower-dimensional neighbor.

  3. 3.

Provide a discrete operator $\text{tr }p_{i}$ so that
$\Pi_{L^{2}(\Omega_{i})}$ can project the pressure trace from
$\partial_{j}\Omega_{i}$ to $\Gamma_{j}$ where $\Omega_{i}$ is the higher-
dimensional neighbor.

  4. 4.

Provide a pressure $p_{i}$ so that $\Pi_{L^{2}(\Omega_{i})}$ can project the
pressure to all $\Gamma_{j}$ where $\Omega_{i}$ is the lower-dimensional
neighbor.

  5. 5.

Handle the divergence of vector source terms
$\Pi_{N(\Omega_{i})}(\nabla\cdot{\bm{\mu}}_{j,i})$ from interfaces
$\Gamma_{j}$ where $\Omega_{i}$ is the lower-dimensional neighbor.

  6. 6.

Provide a pressure gradient $u_{i}$ so that $\Pi_{L^{2}(\Omega_{i})}$ can
project the pressure gradient to all $\Gamma_{j}$ where $\Omega_{i}$ is the
lower-dimensional neighbor.

The four first requirements are readily available for any discretization
scheme for elliptic equations. Specifically, we have based our solution
operators on a cell-centered finite volume method termed the multi-point flux
approximation (MPFA) (Aavatsmark,, 2002; Nordbotten and Keilegavlen,, 2020).
Treatment of vector source terms (item 5) is not as natural in primal
discretization schemes such as finite elements, but is easy to include in most
flux-based discretization methods such as e.g. mixed finite elements. We have
employed the approach introduced in Starnoni et al., (2019), which treats the
vector source term as part of the discrete divergence operator, and thereby
provides an expression of the fluxes in terms of jumps in cell-centers vector
sources. Finally, the pressure gradients are discretized as piece wise
constant on each cell from an interpolation of the face cells fluxes (item 6).
We implemented our model in PorePy, an open-source software for simulation of
multiphysics processes in fractured porous media (Keilegavlen et al.,, 2021).

\begin{overpic}[width=108.405pt]{sketch_discretization}
\put(30.0,76.0){\small$\Omega_{h}$}
\put(50.0,60.0){\small$\partial_{j}\Omega_{h}$}
\put(75.0,30.0){\small$\Gamma_{j}$} \put(75.0,0.0){\small$\Omega_{l}$}
\put(60.0,12.0){\small$\Pi_{L^{2}(\Omega_{h})}$}
\put(60.0,42.0){\small$\Pi_{L^{2}(\Omega_{l})}$}
\put(20.0,12.0){\small$\Pi_{N(\Omega_{l})}$}
\put(20.0,42.0){\small$\Pi_{N(\Omega_{h})}$} \end{overpic} Figure 6:
Illustration of a coupling between subdomains. $\Omega_{h}$ and $\Omega_{l}$
are the higher and lower subdomains respectively, $\Gamma_{j}$ is the
interface between the two subdomains, $\partial_{j}\Omega_{h}$ is the portion
of the boundary of $\Omega_{h}$ as seen from $\Gamma_{j}$,
$\Pi_{N(\Omega_{k})}$ is the projection operator from coupling variables on
the mortar grid to each of the subdomains degrees of freedom ($k=h,l$), and
$\Pi_{L^{2}(\Omega_{k})}$ is the projection operator from numerical variables
to coupling variables.

To better understand the structure of the discrete coupling, it is instructive
to write out the coupled system for two subdomains $\Omega_{h}$ and
$\Omega_{l}$ separated by an interface $\Gamma_{j}$ (see Fig. 6). Let
$\overline{p}_{h}$ and $\overline{p}_{l}$, be the vectors of cell-center
pressures in $\Omega_{h}$ and $\Omega_{l}$ respectively, and let
$\overline{\lambda}_{j}$ be the vector of discrete mortar fluxes in
$\Gamma_{j}$. The discrete coupled system in absence of external sources can
then be represented on the generic form

$\begin{bmatrix}A_{h}&0&G_{h}\Pi_{N(\Omega_{h})}\\\
0&A_{l}&B_{l}\Pi_{N(\Omega_{l})}+J_{l}\Pi_{N(\Omega_{l})}T_{j}\\\
-\Pi_{L^{2}(\Omega_{h})}P_{h}&\Pi_{L^{2}(\Omega_{l})}P_{l}+T_{j}\Pi_{L^{2}(\Omega_{l})}R_{l}&D_{j}\\\
\end{bmatrix}\begin{bmatrix}\overline{p}_{h}\\\ \overline{p}_{l}\\\
\overline{\lambda}_{j}\end{bmatrix}=\begin{bmatrix}0\\\ 0\\\ 0\end{bmatrix}.$
(69)

The first two rows of the system (69) represent the discretised differential
equations in each subdomain, while the third row is the discretized Darcy’s
law in the direction perpendicular to the interface. Here, $A_{h}$ and $A_{l}$
are the fixed-dimensional discretizations on the subdomains, $G_{h}$ is the
discretization of Neumann boundary conditions on $\Omega_{h}$, $B_{l}$ is the
discretization of source terms in $\Omega_{l}$, $J_{l}$ is the discretization
of the vector source term on $\Omega_{l}$, $T_{j}$ is the discretized
$\mathfrak{K}_{\Omega\Gamma}\mathfrak{K}_{\Gamma\Gamma}^{-1}$ product on
$\Gamma_{j}$, and $\Pi_{N(\Omega_{h})}$ and $\Pi_{N(\Omega_{l})}$ are the
projection operators from coupling variables on the mortar grid to each of the
subdomains degrees of freedom. Furthermore, $P_{h}$ provides a discrete
representation of the pressure trace operator on $\Omega_{h}$, $P_{l}$ gives
the pressure unknowns on $\Omega_{l}$, $R_{l}$ gives the reconstruction of the
pressure gradient on $\Omega_{l}$, and $\Pi_{L^{2}(\Omega_{k})}$ is the
projection operator from numerical variables to coupling variables. Finally,
$D_{j}$ is the discretized inverse normal permeability on $\Gamma_{j}$.

We conclude by making two remarks: firstly, there is no direct coupling
between $\Omega_{h}$ and $\Omega_{l}$ and secondly, global boundary conditions
are left out of the system.

##  4 Numerical examples

We validate the semi-local model and our implementation by a suite of
numerical examples. First, we consider a case with a single fault, and show
how the semi-local model can capture the effects of anisotropic off-diangonal
permeabilities, while the local model fails to do so. Second, we probe the
robustness of our discretization on more complex geometries in 2D and 3D.

###  4.1 Comparison to the equi-dimensional model

In this first example, we compare our reduced model to an equi-dimensional
model. The aim is to highlight the enhanced modelling capabilities of our
formulation with respect to the standard local formulation. With reference to
this latter point, we present results of two test cases: the first one where
the fault has the same off-diagonal permeability on both sides (see Fig. 2.b),
and a second one where different permeability structures are assigned to each
side of the fault (see Fig. 2.c).

\begin{overpic}[width=151.76964pt]{sketch_equi_dim}
\put(5.0,54.0){\small$\Omega_{3}$} \put(63.0,66.0){\small$\Omega_{1}$}
\put(63.0,32.0){\small$\Omega_{2}$} \put(8.0,90.0){\small$h_{out}$}
\put(83.0,90.0){\small$h_{out}$} \put(48.0,6.0){\small$h_{in}$} \end{overpic}
Figure 7: Setup of Case 1.

####  4.1.1 Case 1: homogeneous permeability

We consider a 2D square domain of side $L=1~{}m$ cut by a horizontal fault of
aperture $a=1~{}cm$ located in the middle of the domain. In the mixed-
dimensional setting we therefore have two 2D domains $\Omega_{1}$ and
$\Omega_{2}$ and one 1D fault $\Omega_{3}$, as illustrated in Fig. 7. The
hydraulic conductivity is isotropic homogeneous for the 2D matrix, that is
${\bm{K}}_{j}=K_{j}\mathbf{I}$, with $j=1,2$, while for the fault we consider
the following equi-dimensional full tensor:

${\bm{K}}_{3}=\begin{bmatrix}K_{f,\parallel}&k_{f,t}\\\
k_{f,t}&k_{f,\bot}\end{bmatrix},$ (70)

For simplicity, we take $K_{1}=K_{2}=K_{m}$. Boundary conditions consist of an
applied difference in hydraulic head along the vertical direction and no-flow
conditions elsewhere. In particular, the inlet pressure $h_{in}$ is specified
on the portion of the bottom boundary where $0.25<x<0.75~{}m$, while the
outlet pressure $h_{out}$ is specified on the portion of the top boundary
where $x<0.25~{}\&~{}x>0.75~{}m$ (see Fig. 7). Data for the simulations are
reported in Table 1. We consider as reference solution the solution obtained
with an equi-dimensional model of $N=40k$ structured square cells (mesh size
$dx=5~{}mm$), where the fault is discretized with two rows of 200 elements
each. Then, for the reduced models, we consider triangular grids with
approximately $N=[40,160,700,3k,11k]$ (respectively $N_{f}=[4,8,16,32,64]$
cells for the fault), and report the average $L^{2}$ error in pressure along
the fault

$\varepsilon_{p}=\dfrac{\sqrt{\sum_{i}\Delta_{i}(p_{i}-p_{i,eq})^{2}}}{\sqrt{\sum_{i}\Delta_{i}p_{i,eq}^{2}}},$
(71)

where $\Delta_{i}$ is the size of the fault element in the reduced model, and
$p_{i,eq}$ is calculated from the equi-dimensional model as the mean value of
the two fault cells at each location $x_{i}$:

$p_{i,eq}(x_{i})=\sum_{j=y_{1},y_{2}}p_{ij},$ (72)

where $y_{j}=L/2\pm dx/2$.

Convergence results are shown in Fig. 8a. As Fig. 8a clearly shows, our
formulation presents about first-order convergence rate, while the local
formulation does not converge. This is due to the strong anisotropy of the
fault, which is not captured by the standard local formulation. As a result of
the anisotropy of the fault, the flow will take a preferential direction
towards one of the two inlets, therefore breakig the symmetry of the local
formulation. This is better observed in Fig. 8b showing the pressure
distribution along the fault for the three models. As Fig. 8b clearly shows,
the semi-local and the equi-dimensional models coincide, while the local
formulation exhibits an erroneous symmetric profile.

Table 1: Data for Case 1. Values of the fault hydraulic conductivity are given for the equi-dimensional model, i.e. before scaling. Parameter |  Description | Value
---|---|---
$K_{m}$ |  Matrix hydraulic conductivity | $1~{}m/s$
$K_{f,\parallel}$ |  Fault tangential hydraulic conductivity | $100~{}m/s$
$k_{f,\bot}$ |  Fault normal hydraulic conductivity | $100~{}m/s$
$k_{f,t}$ |  Fault off-diagonal hydraulic conductivity | $80~{}m/s$
$a$ |  Fault aperture | $0.01~{}m$
$L$ |  Side of the square domain | $1~{}m$
$h_{in}$ |  Hydraulic head at the bottom boundary | $10~{}m$
$h_{out}$ |  Hydraulic head at the top boundary | $1~{}m$

(a)

(b)

Figure 8: Case 1: (a) convergence of the average error in pressure within the
fault and (b) pressure distribution along the fault for different methods.

####  4.1.2 Case 2: dual permeability

As a further illustration of the enhanced modeling capabilities of the semi-
local model, we modify the setup used in the previous section to have
different permeability structures on the two sides of the fault. This is
relevant for modeling of geological faults, where the two sides of the fault
may undergo different damage processes. To that end, we divide the fault into
an upper and lower part (see Fig. 9) and assign different permeability
structures to the two sides, that is for $j=1,2$:

${\bm{K}}_{3,j}=\begin{bmatrix}K_{f,\parallel}&k_{f,j,t}\\\
k_{f,j,t}&k_{f,\bot}\end{bmatrix}.$ (73)

In particular, values of $K_{m}$, $K_{f,\parallel}$ and $k_{f,\bot}$ are the
same as those given in Table 1, while $k_{f,1,t}$ and $k_{f,2,t}$ take values
of $50$ and $80~{}m/s$, respectively. The aperture of the fault is set to
$a=2~{}cm$ and we use the same boundary conditions as in Case 1.

Convergence results for the local and semi-local models are shown in Figs.
10a-10b, with the reference solution again computed from an equi-dimensional
model with a grid with 40k cells. As in the previous case, the local model
fails to converge, while the semi-local model exhibits first order convergence
up to the last refinement step. Here, the mesh size is of the same order of
the fault aperture, thus further error reduction cannot be expected due to the
modeling error in the dimension reduction.

\begin{overpic}[width=325.215pt]{sketch_dual_kt}
\put(75.0,26.0){\small${K}_{1}$} \put(75.0,8.0){\small${K}_{2}$}
\put(75.0,15.0){\small${K}_{3,2}$} \put(75.0,19.5){\small${K}_{3,1}$}
\put(102.0,15.0){$\small{1~{}cm}$} \put(102.0,20.0){$\small{1~{}cm}$}
\end{overpic} Figure 9: Setup of Case 2.

(a)

(b)

Figure 10: Case 2: (a) convergence of the average error in pressure within the
fault and (b) pressure distribution along the fault for different methods.

###  4.2 Self-convergence

In this section, we test the robustness of the method on more challenging
fault configurations in 2D and 3D.

####  4.2.1 2D case

We consider the same test case as Case 1 in Boon et al., (2018). The domain is
a unit square including a network of five faults (Fig. 11a). Of these five
faults, one cuts the square domain into two 2D subdomains, denoted as
$\Omega_{1}$ and $\Omega_{2}$, respectively. The faults are numbered for
$j=3,..,7$ and are of two kinds: $\Omega_{3}$ and $\Omega_{4}$ are conductive,
that is $K_{3}=K_{4}=K_{f,1}$, while the other three are blocking, that is
$K_{5}=K_{6}=K_{7}=K_{f,2}$. The hydraulic conductivity is isotropic
homogeneous for the 2D matrix, with $K_{1}=K_{2}=K_{m}$, while for the faults
we consider an equi-dimensional full tensor with $k_{j,t}=0.1K_{j,\parallel}$,
for $j=3,..,7$. Boundary conditions consist of an applied difference in
hydraulic head along the vertical direction and no-flow conditions elsewhere.
Data for the simulations are reported in Table 2. We consider as reference
solution the solution obtained with approximately $N=133k$ cells for the 2D
domain and a total number of $N_{f}=510$ cells for the faults. Then we
consider grids with approximately $N=[300,1k,4k,17k,67k]$ (respectively
$N_{f}=[26,48,93,183,363]$), and report the average $L^{2}$ error in pressure
along the faults.

The convergence results, shown in Fig. 11b, indicate a rate of at least first
order. The test thus confirms the performance of our method also in cases that
involve faults that are intersecting and have low permeability. Both these
features are highly relevant in a geologic setting where fault may have
complex geometry and reduced permeability compared to the host rock.

Table 2: Data for the 2D self-convergence test. Values of the fault hydraulic conductivity are given for the equi-dimensional model, i.e. before scaling. Parameter |  Description | Value
---|---|---
$K_{m}$ |  Matrix hydraulic conductivity | $1~{}m/s$
$K_{f,1,\parallel}$ |  Fault tangential hydraulic conductivity | $100~{}m/s$
$k_{f,1,\bot}$ |  Fault normal hydraulic conductivity | $100~{}m/s$
$k_{f,1,t}$ |  Fault off-diagonal hydraulic conductivity | $10~{}m/s$
$K_{f,2,\parallel}$ |  Fault tangential hydraulic conductivity | $0.01~{}m/s$
$k_{f,2,\bot}$ |  Fault normal hydraulic conductivity | $0.01~{}m/s$
$k_{f,2,t}$ |  Fault off-diagonal hydraulic conductivity | $0.001~{}m/s$
$a$ |  Fault aperture | $0.01~{}m$
$h_{in}$ |  Hydraulic head at the top boundary | $1~{}m$
$h_{out}$ |  Hydraulic head at the bottom boundary | $0~{}m$

\begin{overpic}[width=151.76964pt]{sketch_2d_intersections}
\put(5.0,52.0){\small$\Omega_{1}$} \put(88.0,5.0){\small$\Omega_{2}$}
\put(55.0,80.0){\small$\Omega_{3}$} \put(63.0,66.0){\small$\Omega_{4}$}
\put(75.0,52.0){\small$\Omega_{5}$} \put(25.0,32.0){\small$\Omega_{6}$}
\put(78.0,15.0){\small$\Omega_{7}$} \end{overpic} (a)

(b)

Figure 11: 2D self-convergence test: (a) mixed-dimensional geometry and (b)
convergence of the average error in pressure within the faults.

####  4.2.2 3D case

As a final verification, we consider a 3D case with multiple intersecting
faults. The setup is based on Case 2 in the benchmark study described in Berre
et al., (2020). The domain is a unit cube including a network of 9 faults,
whose intersections divide the cubic domain into several subdomains, as
illustrated in Fig. 12a. These 3D subdomains are grouped into two regions,
where we assigned different permeabilities $K_{m,1}$ and $K_{m,2}$, both
homogeneous and isotropic (see Berre et al., (2020) for a visualization of
these two regions). For the faults we consider full tensors with tangential
permeability ${\bm{K}}_{j,\parallel}=K_{f,\parallel}\mathbf{I}_{\parallel}$,
normal permeability $k_{j,\bot}=k_{f,\bot}$, and off-diagonal permeability
${\bm{k}}_{j,t}=0.1k_{j,t}\mathbf{i}_{\parallel}$. Boundary conditions consist
of an imposed normal flux $q_{in}$ on the portion of the boundary where
$(x,y,z)<0.25~{}m$ and a constant hydraulic head $h_{out}$ on the portion of
the boundary where $(x,y,z)>0.875~{}m$. Data for the simulations are reported
in Table 3. We consider as reference solution the solution obtained with
approximately $N_{3}=85k$ cells for the 3D domain and a total number of
$N_{f}=8364$ cells for all faults. Then we consider
$N_{3}=[500,1k,2k,4k,10k,20k,40k]$ (respectively
$N_{f}=[148,282,384,814,1536,2298,3456]$) and report the average $L^{2}$ error
in pressure along the faults.

Convergence results are shown in Fig. 12b, indicating first order convergence
on average. This confirms the consistency of our implementation also for 3D
problems with complex fault geometries.

Table 3: Data for the 3D self-convergence test. Values of the fault hydraulic conductivity are given for the equi-dimensional model, i.e. before scaling. Parameter |  Description | Value
---|---|---
$K_{m,1}$ |  Matrix hydraulic conductivity | $1~{}m/s$
$K_{m,2}$ |  Matrix hydraulic conductivity | $0.1~{}m/s$
$K_{f,\parallel}$ |  Fault tangential hydraulic conductivity | $1e^{4}~{}m/s$
$k_{f,\bot}$ |  Fault normal hydraulic conductivity | $1e^{4}~{}m/s$
$k_{f,t}$ |  Fault off-diagonal hydraulic conductivity | $1e^{3}~{}m/s$
$a$ |  Fault aperture | $1e^{-4}~{}m$
$q_{in}$ |  Normal flux at the inflow boundary | $-1~{}m/s$
$h_{out}$ |  Hydraulic head at the outflow boundary | $1~{}m$

(a)

(b)

Figure 12: 3D self-convergence test: (a) mixed-dimensional geometry and (b)
convergence of the average error in pressure within the faults.

##  5 Conclusions

We presented an improved framework to modelling and discretizing flow in
generally anisotropic porous media with thin inclusions, within the context of
mixed-dimensional partial differential equations. Our model considers a full
permeability tensor for the inclusions, resulting in additional terms arising
in our formulation as compared to existing local discretizations. We expect
our model to be important for modeling of flow in faulted porous media,
however the methods proposed herein can be in any case applied to models of
fractures, in fact our full-permeability model naturally reduces to the
existing models of fracture-matrix flow when the off-diagonal components of
the inclusion permeability tensor are set to zero.

We provided numerical examples showing convergence of the method for both 2D
and 3D faulted porous media. In particular, we provided numerical evidence
that, as opposed to existing local discretizations, our model is capable of
simulating the anisotropic behaviour of the faults near damage zone.

We remark that, in the spirit of flux-mortars coupling schemes, our
formulation is independent of the discretization methods used to discretize
the flow equations in the porous matrix and the faults. However, we only
showed results obtained using a multi-point flux finite volume approach.
Nevertheless, the formulation also applies to other discretization methods,
e.g. mixed finite elements.

Acknowledgements
This work forms part of Norwegian Research Council project 250223. Data will
be made public on Zenodo at the time of publication.

## References

  * Aavatsmark, (2002) Aavatsmark, I. (2002).  An introduction to multipoint flux approximations for quadrilateral grids.  Computational Geosciences, 6(3-4):405–432.
  * Bear, (1979) Bear, J. (1979).  Hydraulics of groundwater new york.  Mc GrawHill Inc.
  * Berre et al., (2020) Berre, I., Boon, W. M., Flemisch, B., Fumagalli, A., Gläser, D., Keilegavlen, E., Scotti, A., Stefansson, I., Tatomir, A., Brenner, K., et al. (2020).  Verification benchmarks for single-phase flow in three-dimensional fractured porous media.  Advances in Water Resources, 147:103759.
  * Bödvarsson and Tsang, (1982) Bödvarsson, G. S. and Tsang, C. F. (1982).  Injection and thermal breakthrough in fractured geothermal reservoirs.  Journal of Geophysical Research: Solid Earth, 87(B2):1031–1048.
  * Boon et al., (2020) Boon, W. M., Nordbotten, J. M., and Vatne, J. E. (2020).  Functional analysis and exterior calculus on mixed-dimensional geometries.  Annali di Mathematica.
  * Boon et al., (2018) Boon, W. M., Nordbotten, J. M., and Yotov, I. (2018).  Robust discretization of flow in fractured porous media.  SIAM Journal on Numerical Analysis, 56(4):2203–2233.
  * Brenner et al., (2017) Brenner, K., Hennicker, J., Masson, R., and Samier, P. (2017).  Gradient discretization of hybrid-dimensional darcy flow in fractured porous media with discontinuous pressures at matrix–fracture interfaces.  IMA Journal of Numerical Analysis, 37(3):1551–1585.
  * Brenner et al., (2018) Brenner, K., Hennicker, J., Masson, R., and Samier, P. (2018).  Hybrid-dimensional modelling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions.  Journal of Computational Physics, 357:100–124.
  * Cao et al., (2016) Cao, P., Liu, J., and Leong, Y.-K. (2016).  A fully coupled multiscale shale deformation-gas transport model for the evaluation of shale gas extraction.  Fuel, 178:103–117.
  * Farmer, (2002) Farmer, C. (2002).  Upscaling: a review.  International journal for numerical methods in fluids, 40(1-2):63–78.
  * Flemisch et al., (2018) Flemisch, B., Berre, I., Boon, W., Fumagalli, A., Schwenck, N., Scotti, A., Stefansson, I., and Tatomir, A. (2018).  Benchmarks for single-phase flow in fractured porous media.  Advances in Water Resources, 111:239–258.
  * Formaggia et al., (2018) Formaggia, L., Scotti, A., and Sottocasa, F. (2018).  Analysis of a mimetic finite difference approximation of flows in fractured porous media.  ESAIM: Mathematical Modelling and Numerical Analysis, 52(2):595–630.
  * Fossen et al., (2005) Fossen, H., Johansen, T. E. S., Hesthammer, J., and Rotevatn, A. (2005).  Fault interaction in porous sandstone and implications for reservoir management; examples from southern utah.  AAPG bulletin, 89(12):1593–1606.
  * Fossen et al., (2007) Fossen, H., Schultz, R. A., Shipton, Z. K., and Mair, K. (2007).  Deformation bands in sandstone: a review.  Journal of the Geological Society, 164(4):755–769.
  * Fumagalli and Keilegavlen, (2019) Fumagalli, A. and Keilegavlen, E. (2019).  Dual virtual element methods for discrete fracture matrix models.  Oil & Gas Science and Technology–Revue d’IFP Energies nouvelles, 74:41.
  * Fumagalli and Scotti, (2019) Fumagalli, A. and Scotti, A. (2019).  A multi-layer reduced model for flow in porous media with a fault and surrounding damage zones.  arXiv preprint arXiv:1903.01117.
  * Helmig et al., (1997) Helmig, R. et al. (1997).  Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems. Springer-Verlag.
  * Hesthammer et al., (2000) Hesthammer, J., Johansen, T., and Watts, L. (2000).  Spatial relationships within fault damage zones in sandstone.  Marine and Petroleum Geology, 17(8):873–893.
  * Johansen and Fossen, (2008) Johansen, T. E. S. and Fossen, H. (2008).  Internal geometry of fault damage zones in interbedded siliciclastic sediments.  Geological Society, London, Special Publications, 299(1):35–56.
  * Johnson et al., (2009) Johnson, S., Morris, J., et al. (2009).  Hydraulic fracturing mechanisms in carbon sequestration applications.  In 43rd US Rock Mechanics Symposium & 4th US-Canada Rock Mechanics Symposium. American Rock Mechanics Association.
  * Karimi-Fard et al., (2003) Karimi-Fard, M., Durlofsky, L. J., Aziz, K., et al. (2003).  An efficient discrete fracture model applicable for general purpose reservoir simulators.  In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers.
  * Keilegavlen et al., (2021) Keilegavlen, E., Berge, R., Fumagalli, A., Starnoni, M., Stefansson, I., Varela, J., and Berre, I. (2021).  Porepy: an open-source software for simulation of multiphysics processes in fractured porous media.  Computational Geosciences, 25(1):243–265.
  * Liu et al., (2016) Liu, R., Li, B., Jiang, Y., and Huang, N. (2016).  Mathematical expressions for estimating equivalent permeability of rock fracture networks.  Hydrogeology Journal, 24(7):1623–1649.
  * Martin et al., (2005) Martin, V., Jaffré, J., and Roberts, J. E. (2005).  Modeling fractures and barriers as interfaces for flow in porous media.  SIAM Journal on Scientific Computing, 26(5):1667–1691.
  * Nagelhout and Roest, (1997) Nagelhout, A. and Roest, J. (1997).  Investigating fault slip in a model of an underground gas storage facility.  International Journal of Rock Mechanics and Mining Sciences, 34(3-4):212–e1.
  * Nordbotten et al., (2019) Nordbotten, J. M., Boon, W. M., Fumagalli, A., and Keilegavlen, E. (2019).  Unified approach to discretization of flow in fractured porous media.  Computational Geosciences, 23(2):225–237.
  * Nordbotten and Celia, (2011) Nordbotten, J. M. and Celia, M. A. (2011).  Geological storage of CO2: modeling approaches for large-scale simulation.  John Wiley & Sons.
  * Nordbotten and Keilegavlen, (2020) Nordbotten, J. M. and Keilegavlen, E. (2020).  An introduction to multi-point flux (mpfa) and stress (mpsa) finite volume methods for thermo-poroelasticity.  arXiv preprint arXiv:2001.01990.
  * Oda, (1985) Oda, M. (1985).  Permeability tensor for discontinuous rock masses.  Geotechnique, 35(4):483–495.
  * Reichenberger et al., (2006) Reichenberger, V., Jakobs, H., Bastian, P., and Helmig, R. (2006).  A mixed-dimensional finite volume method for two-phase flow in fractured porous media.  Advances in water resources, 29(7):1020–1036.
  * Sævik et al., (2013) Sævik, P. N., Berre, I., Jakobsen, M., and Lien, M. (2013).  A 3d computational study of effective medium methods applied to fractured media.  Transport in porous media, 100(1):115–142.
  * Sandve et al., (2012) Sandve, T. H., Berre, I., and Nordbotten, J. M. (2012).  An efficient multi-point flux approximation method for discrete fracture–matrix simulations.  Journal of Computational Physics, 231(9):3784–3800.
  * Starnoni et al., (2019) Starnoni, M., Berre, I., Keilegavlen, E., and Nordbotten, J. (2019).  Consistent mpfa discretization for flow in the presence of gravity.  Water Resources Research.
  * Wilson et al., (2020) Wilson, P., Smith, S., Povey, D., and Harris, S. (2020).  Ranking and selecting fault models using flow-indicator fault properties and simple streamline simulations.  Petroleum Geoscience.
  * Yortsos, (1995) Yortsos, Y. C. (1995).  A theoretical analysis of vertical flow equilibrium.  Transport in Porous Media, 18(2):107–129.