# Call for participation: Verification benchmarks for single-phase flow in
three-dimensional fractured porous media

Inga Berre  Department of Mathematics, University of Bergen, Allégaten 41,
5007 Bergen, Norway  Christian Michelsen Research, Bergen, Norway  Wietse Boon
Department of Hydromechanics and Modelling of Hydrosystems, University of
Stuttgart, Pfaffenwaldring 61, 70569 Stuttgart, Germany  Bernd Flemisch
Department of Hydromechanics and Modelling of Hydrosystems, University of
Stuttgart, Pfaffenwaldring 61, 70569 Stuttgart, Germany<EMAIL_ADDRESS>stuttgart.de  Alessio Fumagalli  Department of Mathematics, University of
Bergen, Allégaten 41, 5007 Bergen, Norway  Dennis Gläser  Department of
Hydromechanics and Modelling of Hydrosystems, University of Stuttgart,
Pfaffenwaldring 61, 70569 Stuttgart, Germany  Eirik Keilegavlen  Department of
Mathematics, University of Bergen, Allégaten 41, 5007 Bergen, Norway  Anna
Scotti  Laboratory for Modeling and Scientific Computing MOX, Politecnico di
Milano, p.za Leonardo da Vinci 32, 20133 Milano, Italy  Ivar Stefansson
Department of Mathematics, University of Bergen, Allégaten 41, 5007 Bergen,
Norway  Alexandru Tatomir  Department of Applied Geology, Geosciences Center,
University of Göttingen, Goldschmidtstrasse 3, 37077 Göttingen, Germany

###### Abstract

This call for participation proposes four benchmark tests to verify and
compare numerical schemes to solve single-phase flow in fractured porous
media. With this, the two-dimensional suite of benchmark tests presented by
Flemisch et al. [1] is extended to include three-dimensional problems.
Moreover, transport simulations are included as a means to compare
discretization methods for flow. With this publication, we invite researchers
to contribute to the study by providing results to the test cases based on
their applied discretization methods.

##  1 Introduction

In this call for participation, we propose four benchmark cases for the
verification of computational models for single-phase flow in three-
dimensional, fractured porous media. It can be seen as an extension of [1],
where only two-dimensional test cases were considered. The current proposal
also presents benchmarks where simulations of linear tracer transport are used
as a means to compare flow discretizations. However, the process of collecting
and comparing results is different. While the publication [1] contains both
the case descriptions and a comparison of several computational models, this
time we aim for a two-stage process. In particular, the present document
describes the benchmark cases and serves as a call for participation to the
interested scientific community. After a phase of collection and discussion of
the results among the participants, we aim for a second publication that
presents a comparison of the participating models. The case descriptions are
accompanied by data in form of geometry descriptions, simulation results and
plotting scripts, all available in the Git repository https://git.iws.uni-
stuttgart.de/benchmarks/fracture-flow-3d.git. While access to this repository
will be restricted to the benchmark participants during the phase of
collection and comparison of the results, all data will be made publicly
available upon submission of the comparison publication.

The rest of this document is organized as follows. In Section 2, we describe
the mathematical models for fluid flow and transport in fractured porous
media. Section 3 briefly describes the participating discretizations proposed
by the authors as well as the discretization of the transport problem. Four
test cases are introduced in Section 4 in form of domain and boundary
condition descriptions as well as pointers to the relevant data in the Git
repository. The final Section 5 provides instructions on how to participate in
the benchmark study.

##  2 Mathematical models

We introduce two models for flow and transport in fractured media. First, the
flow model is presented in the conventional, equi-dimensional setting,
allowing for a natural introduction to the physical parameters. From this
formulation, we derive the mixed-dimensional model through appropriate
reduction of the equations. Note that the mixed-dimensional model forms the
focus of this study. Finally, we introduce the equi- and mixed-dimensional
transport models.

###  2.1 The equi-dimensional flow model

We consider a steady-state incompressible single-phase flow through a porous
medium assumed to be described by Darcy’s law, resulting in the governing
system of equations

$\displaystyle\begin{aligned} \bm{u}+\mathbb{K}\nabla{h}&=0,\\\
\nabla\cdot\bm{u}&={q},\end{aligned}\quad\text{in }\Lambda.$ (1a) Here,
$\bm{u}$ denotes the macroscopic fluid velocity in
$\mathrm{m}\text{/}\mathrm{s}$ whereas $\mathbb{K}$ and ${h}$ stand for
hydraulic conductivity and pressure head measured in
$\mathrm{m}\text{/}\mathrm{s}$ and $\mathrm{m}\text{/}$, respectively; $q$
represents a source/sink term measured in $\text{/}\mathrm{s}$, and the
regular domain $\Lambda\subset\mathbb{R}^{3}$ will be called the equi-
dimensional domain. Coupled to (1a), we consider boundary conditions on the
boundary $\partial\Lambda$ of $\Lambda$, namely $\displaystyle\begin{aligned}
{h}|_{{\partial\Lambda_{{\footnotesize{h}}}}}&=\overline{{h}}&\quad\text{on
}{\partial\Lambda_{{\footnotesize{h}}}},\\\
\bm{u}\cdot\bm{n}|_{{\partial\Lambda_{{\footnotesize
u}}}}&=\overline{u}&\text{on }{\partial\Lambda_{{\footnotesize
u}}}.\end{aligned}$ (1b)

We assume
$\partial\Lambda={\partial\Lambda_{{\footnotesize{h}}}}\cup{\partial\Lambda_{{\footnotesize
u}}}$,
${\partial\Lambda_{{\footnotesize{h}}}}\cap{\partial\Lambda_{{\footnotesize
u}}}=\emptyset$, and $|{\partial\Lambda_{{\footnotesize{h}}}}|>0$. In (1b)
$\cdot|_{A}$ is a suitable trace operator on $A\subset\partial\Lambda$,
depending on the quantity at hand. Also, $\overline{{h}}$ indicates the
pressure head imposed on the boundary
${\partial\Lambda_{{\footnotesize{h}}}}$, while $\overline{u}$ is the
prescribed Darcy velocity normal to the boundary
${\partial\Lambda_{{\footnotesize u}}}$ with respect to the outer unit normal
vector $\bm{n}$.

Problem (1) can be recast in its primal formulation, obtaining the equations

$\displaystyle\begin{aligned}
-\nabla\cdot\mathbb{K}\nabla{h}&={q}&\quad\text{in }\Lambda,\\\
{h}|_{{\partial\Lambda_{{\footnotesize{h}}}}}&=\overline{{h}}&\text{on
}{\partial\Lambda_{{\footnotesize{h}}}},\\\
-\mathbb{K}\nabla{h}\cdot\bm{n}|_{{\partial\Lambda_{{\footnotesize
u}}}}&=\overline{u}&\text{on }{\partial\Lambda_{{\footnotesize
u}}}.\end{aligned}$ (2)

Problem (1) and (2) are equivalent, however, different numerical schemes are
based on either of the two formulations. Under regularity assumptions on
$\Lambda$ and the data, the previous problems admit a unique weak solution. We
refer to [2, 3, 4, 5] for more details.

We assume that $\Lambda$ contains several fractures, i.e., inclusions in the
domain. The fracture walls are considered planar and smooth, and the fractures
have two distinguishing features: (1) The fracture thickness, which we measure
by the aperture, denoted by $\varepsilon$, is small compared to the extension
of the fracture. (2) The fracture hydraulic conductivity may differ
significantly from that of the rest of $\Lambda$, implying that the fractures
may have significant impact on the flow in $\Lambda$.

We furthermore make the assumption that the principal directions of the local
hydraulic conductivities are aligned with the orientation of the fractures. In
particular, the hydraulic conductivity in the matrix ($\mathbb{K}_{3}$), the
fractures ($\mathbb{K}_{2}$), as well as in the intersections between two
fractures ($\mathbb{K}_{1}$) and at the crossings of intersections
($\mathbb{K}_{0}$), can be decomposed in the following way:

$\displaystyle\mathbb{K}_{3}$ $\displaystyle=K_{3}^{eq},$
$\displaystyle\mathbb{K}_{2}$
$\displaystyle=\begin{bmatrix}K_{2}^{eq}&\begin{matrix}0\\\ 0\end{matrix}\\\
\begin{matrix}0&0\end{matrix}&\kappa_{2}^{eq}\end{bmatrix},$
$\displaystyle\mathbb{K}_{1}$ $\displaystyle=\begin{bmatrix}K_{1}^{eq}&0&0\\\
0&\kappa_{1}^{eq}&0\\\ 0&0&\kappa_{1}^{eq}\\\ \end{bmatrix},$
$\displaystyle\mathbb{K}_{0}$
$\displaystyle=\begin{bmatrix}\kappa_{0}^{eq}&0&0\\\ 0&\kappa_{0}^{eq}&0\\\
0&0&\kappa_{0}^{eq}\\\ \end{bmatrix}.$ (3)

Here, $K_{d}^{eq}$ and $\kappa_{d}^{eq}$, for different values of $d$, denote
the tangential and normal hydraulic conductivity, respectively. Thus,
$K_{d}^{eq}$ is an elliptic $(d\times d)$-tensor function whereas
$\kappa_{d}^{eq}$ is a positive scalar function. Note that the subscript $d$
indicates that the features will be represented by $d$-dimensional objects in
the reduced model, as derived in the next section. The superscript $eq$, on
the other hand, indicates that these quantities are related to the equi-
dimensional model.

###  2.2 Mixed-dimensional flow model

The small aperture of the fractures justifies a reduction of dimensionality
procedure to $\Lambda$ where fractures and their intersections are
approximated by lower-dimensional objects. For more details on the derivation
refer to [6, 7, 8, 9, 10, 11, 12, 13, 14, 15].

Here, we use $\Omega$ to denote the mixed-dimensional decomposition of
$\Lambda$. Let $\Omega$ with outer boundary $\partial\Omega$ be composed of a
possibly unconnected $3$-dimensional domain $\Omega_{3}$ which represents the
rock matrix. Furthermore, $\Omega$ contains up to $3$ lower-dimensional, open
subdomains, namely the fracture planes $\Omega_{2}$, their intersection lines
$\Omega_{1}$ and intersection points $\Omega_{0}$. For compatibility, we
assume that $\Omega_{d}\not\subset\Omega_{d^{\prime}}$ for all $d^{\prime}>d$.
Finally, we introduce $\Gamma_{d}=\Omega_{d}\cap\partial\Omega_{d+1}$ as the
set of $d$-interfaces between inter-dimensional sub-domains $\Omega_{d}$ and
$\Omega_{d+1}$ endowed with a normal unit vector $\bm{n}$ pointing outward
from $\Omega_{d+1}$.

Remaining consistent with the notation convention above, data and unknowns
will also be annotated with a subscript related to the dimension. As a first
example, on a $d$-dimensional feature $\Omega_{d,i}\subseteq\Omega_{d}$ with
counting index $i$, let $\varepsilon_{d,i}$ denote the cross-sectional volume,
area, or length of the corresponding physical domain for $d=0,..,2$
respectively. It naturally has the unit of measure
${\mathrm{m}}^{3-d}\text{/}$ and is extended as non-dimensional unity in
$\Omega_{3}$. Moreover, we introduce for each $d$-feature with index $i$, a
typical length $a_{d,i}$ such that $\varepsilon_{d,i}=a_{d,i}^{3-d}$. In the
continuation, we will omit the subscript $i$ if no ambiguity arises.

We continue this subsection by first presenting the reduced model associated
with (1) in the two-dimensional fractures $\Omega_{2}$ followed by its
generalization for all $d=0,\dots,3$.

####  2.2.1 Two-dimensional fracture flow

The primary variables in this formulation are the velocity $\bm{u}_{3}=\bm{u}$
and hydraulic head ${h}_{3}={h}$ in the rock matrix $\Omega_{3}$, as well as
the integrated, tangential velocity $\bm{u}_{2}$ and average hydraulic head
${h}_{2}$ in the fracture. These are given pointwise for $x\in\Omega_{2}$ by

$\displaystyle\bm{u}_{2}(x)=\int_{\varepsilon_{2}(x)}\bm{u}_{\|}\quad\text{and}\quad{h}_{2}(x)=\dfrac{1}{\varepsilon_{2}(x)}\int_{\varepsilon_{2}(x)}{h}.$

Here, $\bm{u}_{\|}$ denotes the components of $\bm{u}$ tangential to
$\Omega_{2}$. The integrals are computed in the normal direction of the
fracture and thus the corresponding units of measurement are
${\mathrm{m}}^{2}\text{/}\mathrm{s}$ and $\mathrm{m}\text{/}$ for $\bm{u}_{2}$
and $h_{2}$, respectively.

Let us derive the reduced Darcy’s law by averaging and the mass balance
equation by integration over the direction normal to the fractures. Recall
that the vector $\bm{n}$ here refers to the normal unit vector pointing from
$\Omega_{3}$ into $\Omega_{2}$.

$\displaystyle\begin{aligned}
\frac{1}{\varepsilon_{2}}\bm{u}_{2}+K_{2}^{eq}\nabla_{2}{h}_{2}&=0\\\
\nabla_{2}\cdot\bm{u}_{2}-\left\llbracket{\bm{u}_{3}\cdot\bm{n}}\right\rrbracket&={q}_{2}\end{aligned}\quad\text{in
}\Omega_{2},$ (4a) with $\nabla_{2}$ the del-operator in the tangential
directions and ${q}_{2}$ representing the integrated source term, i.e.
${q}_{2}(s)=\int_{\varepsilon_{2}(s)}{q}$. Note that here, we have assumed
$K_{2}^{eq}$ to be constant in the direction normal to $\Omega_{2}$. The jump
operator is defined as
$\left\llbracket{\bm{u}_{3}\cdot\bm{n}}\right\rrbracket|_{\Omega_{d}}=\sum(\bm{u}_{3}\cdot\bm{n}|_{\Gamma_{2}})$,
therewith representing the mass exchange between fracture and matrix. In
particular, for each subdomain $\Omega_{2,i}\subseteq\Omega_{2}$, we sum over
all flux contributions over sections of $\Gamma_{2}$ which coincide
geometrically with $\Omega_{2,i}$. These fluxes are assumed to satisfy the
following Darcy-type law given by a finite difference between the hydraulic
head in $\Omega_{2}$ and on $\Gamma_{2}$:
$\displaystyle\bm{u}_{3}\cdot\bm{n}+\kappa_{2}^{eq}\frac{2}{a_{d}}({h}_{2}-{h}_{3})=0\quad\text{on
}\Gamma_{2}.$ (4b)

Note that to be mathematically precise, each term in this equation represents
an appropriate trace or projection of the corresponding variable onto
$\Gamma_{2}$.

####  2.2.2 Generalized flow model

Next, we generalize the equations described above to domains of all
dimensions, thus including the intersection lines and points. For that
purpose, we introduce the integrated velocity $\bm{u}_{d}$ for $d=1$ and
average hydraulic head ${h}_{d}$ with $d=0,1$ given pointwise for
$x\in\Omega_{d}$ by

$\displaystyle\bm{u}_{1}(x)=\int_{\varepsilon_{1}(x)}\bm{u}_{\|}\quad\text{and}\quad{h}_{d}(x)=\dfrac{1}{\varepsilon_{d}(x)}\int_{\varepsilon_{d}(x)}{h},\text{
for }d=0,1.$

Again, $\bm{u}_{\|}$ denotes the components of $\bm{u}$ tangential to
$\Omega_{1}$. The integrals are computed in the directions normal to the
intersection lines $\Omega_{1}$, and points $\Omega_{0}$. The corresponding
units of measurement are therefore ${\mathrm{m}}^{3}\text{/}\mathrm{s}$ and
$\mathrm{m}\text{/}$ for $\bm{u}_{1}$ and $h_{d}$, respectively. The analogues
of (4a) on these lower-dimensional manifolds are then given by

$\displaystyle\begin{aligned}
\frac{1}{\varepsilon_{1}}\bm{u}_{1}+K_{1}^{eq}\nabla_{1}{h}_{1}&=0\\\
\nabla_{1}\cdot\bm{u}_{1}-\left\llbracket{\bm{u}_{2}\cdot\bm{n}}\right\rrbracket&={q}_{1}\quad\text{in
}\Omega_{1},\\\
-\left\llbracket{\bm{u}_{1}\cdot\bm{n}}\right\rrbracket&={q}_{0}\quad\text{in
}\Omega_{0}.\end{aligned}$ (5)

Here, $\nabla_{1}$ denotes the del-operator, i.e. the derivative, in
$\Omega_{1}$. Moreover, the linear jump operator
$\left\llbracket{\cdot}\right\rrbracket$ is naturally generalized to
$\left\llbracket{\bm{u}_{d+1}\cdot\bm{n}}\right\rrbracket|_{\Omega_{d}}=\sum(\bm{u}_{d+1}\cdot\bm{n}|_{\Gamma_{d}})$,
where we for each subdomain $\Omega_{d,i}\subseteq\Omega_{d}$ sum over all
flux contributions over sections of $\Gamma_{d}$ which coincide geometrically
with $\Omega_{d,i}$. Finally $q_{1}$ and $q_{0}$ correspond to the integrated
source terms in the intersection lines and points, respectively.

Due to our choice of defining $\bm{u}_{d}$ as the integrated velocity, a
scaling with $\varepsilon_{d+1}$ appears in the equation governing the flux
across $\Gamma_{d}$:

$\displaystyle\frac{1}{\epsilon_{d+1}}\bm{u}_{d+1}\cdot\bm{n}+\kappa_{d}^{eq}\frac{2}{a_{d}}({h}_{d}-{h}_{d+1})=0\quad\text{on
}\Gamma_{d},\ d=0,1.$ (6)

Recalling that $\epsilon_{3}=1$, it now follows that the effective tangential
and normal hydraulic conductivities are given by:

$\displaystyle K_{d}$ $\displaystyle=\varepsilon_{d}K_{d}^{eq},$
$\displaystyle\text{in }\Omega_{d},\ d$ $\displaystyle=1,\dots,3$ (7a)
$\displaystyle\kappa_{d}$
$\displaystyle=\varepsilon_{d+1}\frac{2}{a_{d}}\kappa_{d}^{eq},$
$\displaystyle\text{on }\Gamma_{d},\ d$ $\displaystyle=0,\dots,2.$ (7b)

From these definitions, it is clear that the units of $K_{d}$ and $\kappa_{d}$
are ${\mathrm{m}}^{4-d}\text{/}\mathrm{s}$ and
${\mathrm{m}}^{2-d}\text{/}\mathrm{s}$, respectively.

Collecting the above equations, we obtain the generalization of system (4) to
subdomains of all dimensions. The resulting system consists of Darcy’s law in
both tangential and normal directions followed by the mass conservation
equations:

$\displaystyle\bm{u}_{d}+K_{d}\nabla_{d}{h}_{d}$ $\displaystyle=0,$ in
$\displaystyle\Omega_{d},\ d=1,\dots,3,$ (8a)
$\displaystyle\bm{u}_{d+1}\cdot\bm{n}+\kappa_{d}({h}_{d}-{h}_{d+1})$
$\displaystyle=0,$ on $\displaystyle\Gamma_{d},\ d=0,\dots,2,$ (8b)
$\displaystyle\nabla_{d}\cdot\bm{u}_{3}$ $\displaystyle={q}_{3},$ in
$\displaystyle\Omega_{3},$ (8c)
$\displaystyle\nabla_{d}\cdot\bm{u}_{d}-\left\llbracket{\bm{u}_{d+1}\cdot\bm{n}}\right\rrbracket$
$\displaystyle={q}_{d},$ in $\displaystyle\Omega_{d},\ d=1,2,$ (8d)
$\displaystyle-\left\llbracket{\bm{u}_{1}\cdot\bm{n}}\right\rrbracket$
$\displaystyle={q}_{0},$ in $\displaystyle\Omega_{0}.$ (8e)

The source term is given by ${q}_{3}$ for the rock matrix and
${q}_{d}(x)=\int_{\varepsilon_{d}(x)}{q}$ measured in
${\mathrm{m}}^{3-d}\text{/}\mathrm{s}$.

The system (8) is then compactly described by:

$\displaystyle\bm{u}_{d}+K_{d}\nabla_{d}{h}_{d}$ $\displaystyle=0,$
$\displaystyle\text{in }\Omega_{d},\ d$ $\displaystyle=1,\dots,3,$ (9a)
$\displaystyle\bm{u}_{d+1}\cdot\bm{n}+\kappa_{d}({h}_{d}-{h}_{d+1})$
$\displaystyle=0,$ $\displaystyle\text{on }\Gamma_{d},\ d$
$\displaystyle=0,\dots,2,$ (9b)
$\displaystyle\nabla_{d}\cdot\bm{u}_{d}-\left\llbracket{\bm{u}_{d+1}\cdot\bm{n}}\right\rrbracket$
$\displaystyle={q}_{d},$ $\displaystyle\text{in }\Omega_{d},\ d$
$\displaystyle=0,\dots,3.$ (9c)

in which the non-physical $\bm{u}_{4}$ and $\bm{u}_{0}$ are understood as
zero. The associated boundary conditions are inherited from the
equidimensional model with the addition of a no-flux condition at embedded
fracture endings:

$\displaystyle{h}_{d}$ $\displaystyle=\overline{{h}}$ $\displaystyle\text{on
}\partial\Omega_{d}\cap{\partial\Lambda_{{\footnotesize{h}}}},\ d$
$\displaystyle=0,\dots,3,$ (10a) $\displaystyle\bm{u}_{d}\cdot\bm{n}$
$\displaystyle=\varepsilon_{d}\overline{u}$ $\displaystyle\text{on
}\partial\Omega_{d}\cap{\partial\Lambda_{{\footnotesize u}}},\ d$
$\displaystyle=1,\dots,3,$ (10b) $\displaystyle\bm{u}_{d}\cdot\bm{n}$
$\displaystyle=0$ $\displaystyle\text{on
}\partial\Omega_{d}\backslash(\Gamma_{d-1}\cup\partial\Lambda),\ d$
$\displaystyle=1,\dots,3.$ (10c)

To finish the section, we present the primal formulation of the mixed-
dimensional fracture flow model. In analogy with (2), this formulation is
derived by substituting the Darcy’s laws (9a), (9b) into the conservation
equation (9c):

$\displaystyle-\nabla_{d}\cdot
K_{d}\nabla_{d}{h}_{d}+\left\llbracket{\kappa_{d}({h}_{d}-{h}_{d+1})}\right\rrbracket$
$\displaystyle={q}_{d},$ $\displaystyle\text{in }\Omega_{d},\ d$
$\displaystyle=0,\dots,3.$ (11)

In this case, we set the divergence term to zero if $d=0$ and the jump term to
zero if $d=3$. The boundary conditions are given by

$\displaystyle{h}_{d}$ $\displaystyle=\overline{{h}}$ $\displaystyle\text{on
}\partial\Omega_{d}\cap{\partial\Lambda_{{\footnotesize{h}}}},\ d$
$\displaystyle=0,\dots,3,$ (12a) $\displaystyle-
K_{d}\nabla_{d}{h}_{d}\cdot\bm{n}$ $\displaystyle=\varepsilon_{d}\overline{u}$
$\displaystyle\text{on }\partial\Omega_{d}\cap{\partial\Lambda_{{\footnotesize
u}}},\ d$ $\displaystyle=1,\dots,3,$ (12b) $\displaystyle-
K_{d}\nabla_{d}{h}_{d}\cdot\bm{n}$ $\displaystyle=0$ $\displaystyle\text{on
}\partial\Omega_{d}\backslash(\Gamma_{d-1}\cup\partial\Lambda),\ d$
$\displaystyle=1,\dots,3.$ (12c)

Many discretization schemes of the mixed-dimensional model described here,
ignore flow in one-dimensional fracture intersections and zero-dimensional
intersections thereof. Although these correspond to discretizing a simpler
model, this is perfectly in line with the proposed study. In particular, it
will be interesting to evaluate differences in the results that can be
attributed to whether these lower-dimensional intersections are included in
the model, and therewith the discretization, or not.

###  2.3 Equi-dimensional transport model

We now consider a scalar quantity $c$ with the unit of measure
${\mathrm{m}}^{-3}\text{/}$, which is transported through the porous medium
subject to the velocity field resulting from the flow model presented in the
previous sections. The purely advective transport of $c$ is described by the
conservation equation:

$\phi\frac{\partial
c}{\partial{t}}+\nabla\cdot\left(c\bm{u}\right)={q}_{c}\quad\text{in
}\Lambda,$ (13)

where $\phi$ is the porosity of the porous medium and ${q}_{c}$ is a
source/sink term for $c$ given in ${\mathrm{m}}^{-3}\text{/}\mathrm{s}$. We
define Dirichlet boundary conditions on those boundary segments where inflow
occurs, i.e.:

$c|_{\partial\Lambda_{\footnotesize c}}=\overline{c}\quad\text{on
}\partial\Lambda_{\footnotesize c},\quad\partial\Lambda_{\footnotesize
c}=\\{x\in\partial\Lambda:\,\,\bm{u}\cdot\bm{n}<0\\},$ (14)

with $\overline{c}$ being the value for $c$ prescribed on the boundary
$\partial\Lambda_{\footnotesize c}$.

###  2.4 Mixed-dimensional transport model

In analogy to Section 2.2, we choose the average value for $c$ as primary
variable, which is defined as $c_{3}=c$ in $\Omega_{3}$ and for the lower
dimensional objects (with $d\leq 2$) as

$c_{d}(s)=\dfrac{1}{\varepsilon_{d}(s)}\int_{\varepsilon_{d}(s)}c.$

Following the derivation of the mixed-dimensional flow model presented in
Section 2.2, the resulting mixed-dimensional transport model reads:

$\varepsilon_{d}\phi_{d}\frac{\partial
c_{d}}{\partial{t}}+\nabla_{d}\cdot\left(c_{d}\bm{u}_{d}\right)-\left\llbracket{\tilde{c}_{d+1}\left(\bm{u}_{d+1}\cdot\bm{n}\right)}\right\rrbracket={q}_{c,d}\quad\text{in
}\Omega_{d},\quad d=0,\dots,3.$ (15)

Note that for $d=0$, the divergence term is void. Here, the porosity is simply
$\phi_{d}=\phi^{eq}$, with unit of measure ${\mathrm{m}}^{-3}\text{/}$, and
$\tilde{c}_{d+1}$ is evaluated on the basis of a first order upwind scheme,
i.e.:

$\tilde{c}_{d+1}=\begin{cases}c_{d+1}&\text{if}\quad\bm{u}_{d+1}\cdot\bm{n}|_{\Gamma_{d}}>0\\\
c_{d}&\text{if}\quad\bm{u}_{d+1}\cdot\bm{n}|_{\Gamma_{d}}<0.\end{cases}$ (16)

As in the flow model, the jump operator represents the sum of the fluxes over
all contributions defined on sections of $\Gamma_{d}$ that coincide
geometrically with $\Omega_{d,i}$.

##  3 Discretization

The intent of the proposed benchmark study is to quantitatively evaluate
different discretization schemes for the mixed-dimensional flow models (9),
(10) and (11), (12). The transport models described above are employed as a
means to further evaluate and compare results for different flow
discretizations. In the repository, the authors have already uploaded results
to the test cases based on the following methods:

  1. 1.

UiB-TPFA: Two-point flux approximation finite volume method (cell-centered
pressure, conforming DFM grid, discontinuous over lower-dimensional features,
mass transfer in intersection lines and points).

  2. 2.

UiB-MPFA: Multi-point flux approximation finite volume method (cell-centered
pressure, conforming DFM grid, discontinuous over lower-dimensional features,
mass transfer in intersection lines and points).

  3. 3.

UiB-RT0: Lowest order Raviart-Thomas mixed finite elements (cell-centered
pressure and face flux, conforming DFM grid, discontinuous over lower-
dimensional features, mass transfer in intersection lines and points).

  4. 4.

UiB-MVEM: Lowest order mixed virtual element method (cell-centered pressure
and face flux, conforming DFM grid, discontinuous over lower-dimensional
features, mass transfer in intersection lines and points).

  5. 5.

USTUTT-MPFA: Multi-point flux approximation finite volume method (cell-
centered pressure, conforming DFM grid, discontinuous over lower-dimensional
features, _no_ mass transfer in intersection lines and points).

Methods 1-4 are all implemented following the mixed-dimensional approach
described in [16].

To evaluate the fluxes produced by each discretization quantitatively, they
are to be inserted into a standard cell-centered first-order-upwind scheme for
the transport, together with the implicit Euler method for temporal
discretization with a fixed time-step that will be prescribed for each test
case.

##  4 Benchmark cases

In the next subsections we introduce the proposed benchmark cases. For each
case, pressure and tracer concentration need to be computed along with several
associated macroscopic metrics. In Subsection 4.1 a problem containing a
single fracture problem is considered. Subsection 4.2 contains a synthetic
network composed by 9 fractures. The example in Subsection 4.3 considers the
geometrically challenging case of almost intersecting fractures, fractures
with small intersections, and other features that a fracture network may
exhibit. Finally, in Subsection 4.4 we study a selection of a few dozen
fractures from a real network. In all the cases the purpose is to validate the
performance of the proposed numerical schemes and model choices.

###  4.1 Case 1: single fracture

##### Description

Figure 1 illustrates the first proposed benchmark example for testing the
discrete fracture flow and transport models in three-dimensional space. The
geometry of the system is slightly modified from the previous works of [17]
and [18].

Figure 1: Conceptual model and geometrical description of the domain for test
case of Subsection 4.1.

The domain, $\Omega$, is a cube-shaped region ($0\text{\,}\mathrm{m}\text{/}$,
$100\text{\,}\mathrm{m}\text{/}$) x ($0\text{\,}\mathrm{m}\text{/}$,
$100\text{\,}\mathrm{m}\text{/}$) x ($0\text{\,}\mathrm{m}\text{/}$,
$100\text{\,}\mathrm{m}\text{/}$) which is crossed by a planar fracture,
$\Omega_{2}$, with a thickness of
$1\text{\times}{10}^{-2}\text{\,}\mathrm{m}\text{/}$. The fracture divides the
domain into the sub-domains $\Omega_{3,1}$ and $\Omega_{3,2}$. Furthermore, an
additional sub-domain $\Omega_{3,3}$ is defined in the lower part of the
domain, representing a heterogeneity within the rock matrix. Inflow into the
system occurs through a narrow band defined by
$\\{$0\text{\,}\mathrm{m}\text{/}$\\}\times($0\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$)\times($90\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$)$.
Similarly, the outlet is a narrow band defined by
$($0\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$)\times\\{$0\text{\,}\mathrm{m}\text{/}$\\}\times($0\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$)$.
Except for the inlet and outlet bands, on which Dirichlet boundary conditions
are set, all other boundaries are treated as no-flow boundaries. For the
hydraulic head, $h_{in}=$4\text{\,}\mathrm{m}\text{/}$$ and
$h_{out}=$1\text{\,}\mathrm{m}\text{/}$$ are set on the inlet and outlet
boundaries, respectively, while for the transport
$c_{in}=$1\text{\times}{10}^{-2}\text{\,}{\mathrm{m}}^{-3}\text{/}$$ is set at
the inlet. The overall simulation time is
$1\text{\times}{10}^{9}\text{\,}\mathrm{s}\text{/}$.

##### Parameters

Matrix hydraulic conductivity $K_{3,1}$, $K_{3,2}$ |  $1\text{\times}{10}^{-6}$$\bm{I}$ | $\mathrm{m}\text{/}\mathrm{s}$
---|---|---
Matrix hydraulic conductivity $K_{3,3}$ |  $1\text{\times}{10}^{-5}$$\bm{I}$ | $\mathrm{m}\text{/}\mathrm{s}$
Fracture effective tangential hydraulic conductivity $K_{2}$ |  $1\text{\times}{10}^{-3}$$\bm{I}$ | $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$
Fracture effective normal hydraulic conductivity $\kappa_{2}$ | $20$ | $\text{/}\mathrm{s}$
Matrix porosity $\phi_{3,1},\phi_{3,2}$ | $2\text{\times}{10}^{-1}$ |
Matrix porosity $\phi_{3,3}$ | $2.5\text{\times}{10}^{-1}$ |
Fracture porosity $\phi_{2}$ | $4\text{\times}{10}^{-1}$ |
Fracture cross-sectional length $\epsilon_{2}$ | $1\text{\times}{10}^{-2}$ | $\mathrm{m}\text{/}$

##### Folder in Git repository

single

##### Results

We consider a constant time step size $\Delta
t=${10}^{7}\text{\,}\mathrm{s}\text{/}$$. The results of three different
simulations with approximately $1k$, $10k$ and $100k$ cells for the 3D domain
are collected and the comparison among the different schemes is then done on
the basis of the following results:

  1. 1.

The number of cells and degrees of freedom, according to the guidelines
described in Subsection 5.4 point 3, saved in the file results.csv;

  2. 2.

the computation of $\int_{\Omega_{3,3}}\phi_{3}c_{3}\,\mathrm{d}x$ for each
time step;

  3. 3.

the computation of
$\int_{\Omega_{2}}\varepsilon_{2}\phi_{2}c_{2}\,\mathrm{d}x$ for each time
step;

  4. 4.

the integrated flux of $c$ across the outlet boundary for each time step;

  5. 5.

the plot of the hydraulic head ${h}_{3}$ in the matrix along the line
$\left($0\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$\right)$-$\left($100\text{\,}\mathrm{m}\text{/}$,$0\text{\,}\mathrm{m}\text{/}$,$0\text{\,}\mathrm{m}\text{/}$\right)$;

  6. 6.

the plot of $c_{3}$ in the matrix, at the final simulation time, along the
line
$\left($0\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$\right)$-$\left($100\text{\,}\mathrm{m}\text{/}$,$0\text{\,}\mathrm{m}\text{/}$,$0\text{\,}\mathrm{m}\text{/}$\right)$;

  7. 7.

the plot of $c_{2}$ within the fracture, at the final simulation time, along
the line
$\left($0\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$,$80\text{\,}\mathrm{m}\text{/}$\right)$-$\left($100\text{\,}\mathrm{m}\text{/}$,$0\text{\,}\mathrm{m}\text{/}$,$20\text{\,}\mathrm{m}\text{/}$\right)$.

In files called dot_refinement_${REFINEMENT}.csv report the results ordered
column-wise in the following way: time, result from point 2, 3, and 4. The
variable ${REFINEMENT} runs from 0 to 2 for increasing refinement. Similarly,
in files called dol_refinement_${REFINEMENT}.csv report the results ordered in
the following way: arc length and result from point 5, arc length and result
from point 6, and arc length and result from point 7. See Subsection 5.4 for
more details on how to report results.

###  4.2 Case 2: regular fracture network

##### Description

The second benchmark is an extension of test case 4.1 from the benchmark study
[1] in a 3D setting. The geometry is slightly modified from [19] to allow a
full 3D solution field. The unit cube domain is
$\Omega=\left($0\text{\,}\mathrm{m}\text{/}$,$1\text{\,}\mathrm{m}\text{/}$\right)^{3}$
and 9 fractures are present, see Appendix 6.1 for the coordinates. A graphical
representation is given by Figure 2.

We define the portion of the boundary where we impose a non-homogeneous flux
condition as $\partial\Omega_{flux}$ and choose
$\partial\Omega_{h}=\\{(x,y,z)\in\partial\Omega:x,y,z>$0.875\text{\,}\mathrm{m}\text{/}$\\}$
and
$\partial\Omega_{flux}=\\{(x,y,z)\in\partial\Omega:x,y,z<$0.25\text{\,}\mathrm{m}\text{/}$\\}$,
respectively. On $\partial\Omega_{h}$ we impose
$\overline{h}=$1\text{\,}\mathrm{m}\text{/}$$, while on
$\partial\Omega_{flux}$ we impose
$\overline{u}=$-1\text{\,}\mathrm{m}\text{/}\mathrm{s}$$. On the rest of the
boundary of $\Omega$ and on the lower-dimensional objects, we impose
homogeneous flux conditions. Two test cases are considered: high-conductive
fractures (${COND} = 0) and low-conductive fractures (${COND} = 1). We define
the following matrix sub-regions

$\displaystyle\Omega_{3,0}=$ $\displaystyle\ \Omega_{3}\setminus\Omega_{3,1}$
$\displaystyle\Omega_{3,1}=$ $\displaystyle\
\\{(x,y,z)\in\Omega_{3}:x>$0.5\text{\,}\mathrm{m}\text{/}$\cap
y<$0.5\text{\,}\mathrm{m}\text{/}$\\}$
$\displaystyle\cup\\{(x,y,z)\in\Omega_{3}:x>$0.75\text{\,}\mathrm{m}\text{/}$\cap$0.5\text{\,}\mathrm{m}\text{/}$<y<$0.75\text{\,}\mathrm{m}\text{/}$\cap
z>$0.5\text{\,}\mathrm{m}\text{/}$\\}$
$\displaystyle\cup\\{(x,y,z)\in\Omega_{3}:$0.625\text{\,}\mathrm{m}\text{/}$<x<$0.75\text{\,}\mathrm{m}\text{/}$\cap$0.5\text{\,}\mathrm{m}\text{/}$<y<$0.625\text{\,}\mathrm{m}\text{/}$\cap$0.5\text{\,}\mathrm{m}\text{/}$<z<$0.75\text{\,}\mathrm{m}\text{/}$\\},$

see the right part of Figure 2, and prescribe two different matrix hydraulic
conductivities.

Figure 2: On the left, representation of the fractures and the outline of the
domain for the test case of Subsection 4.2. On the right, the part of the rock
matrix with low hydraulic conductivity $\Omega_{3,1}$.

For the transport problem, we impose unitary concentration at the inflow
boundary and time interval equal to $[0,0.25]$ $\mathrm{s}\text{/}$.

##### Parameters

| ${COND} = 0 | ${COND} = 1
---|---|---
Matrix hydraulic conductivity $K_{3}|_{\Omega_{3,0}}$ | $\bm{I}$ | $\mathrm{m}\text{/}\mathrm{s}$ | $\bm{I}$ | $\mathrm{m}\text{/}\mathrm{s}$
Matrix hydraulic conductivity $K_{3}|_{\Omega_{3,1}}$ |  $1\text{\times}{10}^{-1}$$\bm{I}$ | $\mathrm{m}\text{/}\mathrm{s}$ |  $1\text{\times}{10}^{-1}$$\bm{I}$ | $\mathrm{m}\text{/}\mathrm{s}$
Fracture effective tangential hydraulic conductivity $K_{2}$ | $\bm{I}$ | ${\mathrm{m}}^{2}\text{/}\mathrm{s}$ |  $1\text{\times}{10}^{-8}$$\bm{I}$ | ${\mathrm{m}}^{2}\text{/}\mathrm{s}$
Fracture effective normal hydraulic conductivity $\kappa_{2}$ | $2\text{\times}{10}^{8}$ | $\text{/}\mathrm{s}$ | $2$ | $\text{/}\mathrm{s}$
Intersection effective tangential hydraulic conductivity $K_{1}$ | $1\text{\times}{10}^{-4}$ | ${\mathrm{m}}^{3}\text{/}\mathrm{s}$ | $1\text{\times}{10}^{-12}$ | ${\mathrm{m}}^{3}\text{/}\mathrm{s}$
Intersection effective normal hydraulic conductivity $\kappa_{1}$ | $2\text{\times}{10}^{4}$ | $\mathrm{m}\text{/}\mathrm{s}$ | $2\text{\times}{10}^{-4}$ | $\mathrm{m}\text{/}\mathrm{s}$
Intersection effective normal hydraulic conductivity $\kappa_{0}$ | $2$ | ${\mathrm{m}}^{2}\text{/}\mathrm{s}$ | $2\text{\times}{10}^{-8}$ | ${\mathrm{m}}^{2}\text{/}\mathrm{s}$
Matrix porosity $\phi_{3}$ | $1\text{\times}{10}^{-1}$ |  | $1\text{\times}{10}^{-1}$ |
Fracture porosity $\phi_{2}$ | $9\text{\times}{10}^{-1}$ |  | $1\text{\times}{10}^{-2}$ |
Intersection porosity $\phi_{1}$ | $9\text{\times}{10}^{-1}$ |  | $1\text{\times}{10}^{-2}$ |
Fracture cross-sectional length $\epsilon_{2}$ | $1\text{\times}{10}^{-4}$ | $\mathrm{m}\text{/}$ | $1\text{\times}{10}^{-4}$ | $\mathrm{m}\text{/}$
Intersection cross-sectional area $\epsilon_{1}$ | $1\text{\times}{10}^{-8}$ | ${\mathrm{m}}^{2}\text{/}$ | $1\text{\times}{10}^{-8}$ | ${\mathrm{m}}^{2}\text{/}$
Intersection cross-sectional volume $\epsilon_{0}$ | $1\text{\times}{10}^{-12}$ | ${\mathrm{m}}^{3}\text{/}$ | $1\text{\times}{10}^{-12}$ | ${\mathrm{m}}^{3}\text{/}$

##### Folder in Git repository

regular

##### Results

Results will be collected for a sequence of 3 simulations, by taking
approximately $500$, $4k$, and $32k$ cells for the 3D domain. For the lower
dimensional objects the number of cells should be chosen accordingly. The
meshes may, if desirable, be obtained using the Gmsh mesh generator using the
file gmsh.geo in the folder regular/geometry of the repository, where the mesh
parameter h can vary in $[0.25,0.12,0.045]$ $\mathrm{m}\text{/}$. The time
interval is divided in 100 equal time steps. The results to be reported are:

  1. 1.

The number of cells and degrees of freedom, according to the guidelines
described in in Subsection 5.4 point 3, saved in two files
results_cond_${COND}.csv;

  2. 2.

for each mesh, the plot of the matrix hydraulic head ${h}_{3}$ along the line
from (0, 0, 0) to $(1,1,1)$;

  3. 3.

for the second level of mesh refinement, the average concentration on each
matrix block during time. The average concentration is defined as
$\int_{\Omega_{3,i}}c_{3}/|\Omega_{3,i}|$${\mathrm{m}}^{-3}\text{/}$, with $i$
the region identifier.

In files called dol_cond_${COND}_refinement_${REFINEMENT}.csv report the
results from point 2. The variable ${REFINEMENT} runs from 0 to 2 for
increasing refinement, while ${COND} should be 0 for the high-conductive case
and 1 for low-conductive. For point 3, the sub-domain identification indexes
$region\\_id$ are defined in Appendix 6.1, and the file color_regions.vtu
shows a graphical representation of these regions. The output is organised in
the following way: two files named dot_cond_${COND}.csv having time as the
first column, and the other columns the averaged concentration in each of the
regions.

###  4.3 Case 3: network with small features

##### Description

This test is designed to probe accuracy in the presence of small geometric
features, that may cause trouble for conforming meshing strategies. The domain
is a box
$\Omega=($0\text{\,}\mathrm{m}\text{/}$,$1\text{\,}\mathrm{m}\text{/}$)\times($0\text{\,}\mathrm{m}\text{/}$,$2.25\text{\,}\mathrm{m}\text{/}$)\times($0\text{\,}\mathrm{m}\text{/}$,$1\text{\,}\mathrm{m}\text{/}$)$,
containing eight fractures. These are defined by their vertexes, given in
Appendix 6.2, and the geometry is depicted in Figure 3.

Figure 3: Representation of the fractures and the outline of the domain for
the test case of Subsection 4.3.

A Gmsh file gmsh.geo representing this geometry is uploaded to the Git
repository in the folder small_features/geometry.

We define in and outlet boundaries as follows:

$\displaystyle\partial\Omega_{N}=\partial\Omega\setminus(\partial\Omega_{in}\cup\partial\Omega_{out})\qquad\partial\Omega_{in}=($0\text{\,}\mathrm{m}\text{/}$,$1\text{\,}\mathrm{m}\text{/}$)\times\\{$0\text{\,}\mathrm{m}\text{/}$\\}\times($1\text{/}3\text{\,}\mathrm{m}\text{/}$,$2\text{/}3\text{\,}\mathrm{m}\text{/}$)$
$\displaystyle\partial\Omega_{out}=\partial\Omega_{out,0}\cup\partial\Omega_{out,1}\quad\partial\Omega_{out,0}=($0\text{\,}\mathrm{m}\text{/}$,$1\text{\,}\mathrm{m}\text{/}$)\times\\{$2.25\text{\,}\mathrm{m}\text{/}$\\}\times($0\text{\,}\mathrm{m}\text{/}$,$1\text{/}3\text{\,}\mathrm{m}\text{/}$)$
$\displaystyle\partial\Omega_{out,1}=($0\text{\,}\mathrm{m}\text{/}$,$1\text{\,}\mathrm{m}\text{/}$)\times\\{$2.25\text{\,}\mathrm{m}\text{/}$\\}\times($2\text{/}3\text{\,}\mathrm{m}\text{/}$,$1\text{\,}\mathrm{m}\text{/}$)$

The boundary conditions for flow are: homogeneous Dirichlet conditions on
$\partial\Omega_{out}$, uniform unit inflow on $\partial\Omega_{in}$, so that
$\int_{\partial\Omega_{in}}\bm{u}_{3}\cdot\bm{n}dS=-1/3$
${\mathrm{m}}^{3}\text{/}\mathrm{s}$, and homogeneous Neumann conditions on
$\partial\Omega_{N}$. For the transport problem, we consider homogeneous
initial conditions and as boundary condition a unit concentration at
$\partial\Omega_{in}$. The overall simulation time is
$1\text{\,}\mathrm{s}\text{/}$.

##### Parameters

Matrix hydraulic conductivity $K_{3}$ | $\bm{I}$ | $\mathrm{m}\text{/}\mathrm{s}$
---|---|---
Fracture effective tangential hydraulic conductivity $K_{2}$ |  $1\text{\times}{10}^{2}$$\bm{I}$ | ${\mathrm{m}}^{2}\text{/}\mathrm{s}$
Fracture effective normal hydraulic conductivity $\kappa_{2}$ | $2\text{\times}{10}^{6}$ | $\text{/}\mathrm{s}$
Intersection effective tangential hydraulic conductivity $K_{1}$ | $1$ | ${\mathrm{m}}^{3}\text{/}\mathrm{s}$
Intersection effective normal hydraulic conductivity $\kappa_{1}$ | $2\text{\times}{10}^{4}$ | $\mathrm{m}\text{/}\mathrm{s}$
Matrix porosity $\phi_{3}$ | $2\text{\times}{10}^{-1}$ |
Fracture porosity $\phi_{2}$ | $2\text{\times}{10}^{-1}$ |
Intersection effective porosity $\phi_{1}$ | $2\text{\times}{10}^{-1}$ |
Fracture cross-sectional length $\epsilon_{2}$ | $1\text{\times}{10}^{-2}$ | $\mathrm{m}\text{/}$
Intersection cross-sectional area $\epsilon_{1}$ | $1\text{\times}{10}^{-4}$ | ${\mathrm{m}}^{2}\text{/}$

##### Folder in Git repository

small_features

##### Results

The simulations should be done on 2 grids with approximately $30k$ and $150k$
3D cells, labeled by ${REFINEMENT} 0 and 1, respectively. The time step size
is $\Delta t=$1\text{\times}{10}^{-2}\text{\,}\mathrm{s}$$. The results to be
reported are:

  1. 1.

The number of cells and degrees of freedom according to the guidelines
described in Subsection 5.4 point 3 in a file named results.csv. In addition,
the total outflow over $\partial\Omega_{out,0}$ and $\partial\Omega_{out,1}$
as columns seven and eight;

  2. 2.

the average hydraulic head ${h}_{3}$ at the inlet, $\partial\Omega_{in}$,
reported as column nine of results.csv. The interpretation of hydraulic head
at the inlet faces may differ between numerical methods; this should be seen
as part of the discretization.

  3. 3.

the hydraulic head ${h}_{3}$ in the matrix along the line
$\left($0.5\text{\,}\mathrm{m}\text{/}$,$1.1\text{\,}\mathrm{m}\text{/}$,$0\text{\,}\mathrm{m}\text{/}$\right)$-$\left($0.5\text{\,}\mathrm{m}\text{/}$,$1.1\text{\,}\mathrm{m}\text{/}$,$1\text{\,}\mathrm{m}\text{/}$\right)$;

  4. 4.

the hydraulic head ${h}_{3}$ in the matrix along
$\left($0\text{\,}\mathrm{m}\text{/}$,$2.15\text{\,}\mathrm{m}\text{/}$,$0.5\text{\,}\mathrm{m}\text{/}$\right)$-$\left($1\text{\,}\mathrm{m}\text{/}$,$2.15\text{\,}\mathrm{m}\text{/}$,$0.5\text{\,}\mathrm{m}\text{/}$\right)$;

  5. 5.

the mean tracer concentrations for each fracture and time step.

In files called dol_line_${LINE}_refinement_${REFINEMENT}.csv report the
results from point 3 and 4, respectively for ${LINE} = 0 and ${LINE} = 1. For
point 5, report in a file dot_refinement_${REFINEMENT}.csv, with time in the
first column, followed by the eight mean values ordered as in Appendix 6.2.

###  4.4 Case 4: a field case

##### Description

The fourth example has a geometry based on a post-processed outcrop from the
island of Algerøyna, outside Bergen, Norway. From the outcrop, 52 fractures
are selected, extruded in the vertical direction and then cut by a bounding
box. The resulting fracture geometry is depicted in Figure 4. The resulting
network has 106 fracture intersections, and multiple fractures intersecting
the domain boundary.

The simulation domain is the box
$\Omega=($-500\text{\,}\mathrm{m}\text{/}$,$350\text{\,}\mathrm{m}\text{/}$)\times($100\text{\,}\mathrm{m}\text{/}$,$1500\text{\,}\mathrm{m}\text{/}$)\times($-100\text{\,}\mathrm{m}\text{/}$,$500\text{\,}\mathrm{m}\text{/}$)$.
The polygons describing the fractures can be found in the file
fracture_network.csv in the folder field. A Gmsh file gmsh.geo describing the
fracture geometry is provided in the subfolder field/geometry.

Figure 4: Representation of the fractures and the outline of the domain for
the test case of Subsection 4.4.

We define in and outlet boundaries as follows:

$\displaystyle\partial\Omega_{N}=\partial\Omega\setminus(\partial\Omega_{in}\cup\partial\Omega_{out}),\qquad\partial\Omega_{in}=\partial\Omega_{in,0}\cup\partial\Omega_{in,1},\qquad\partial\Omega_{out}=\partial\Omega_{out,0}\cup\partial\Omega_{out,1},$
$\displaystyle\partial\Omega_{in,0}=($-500\text{\,}\mathrm{m}\text{/}$,$-200\text{\,}\mathrm{m}\text{/}$)\times\\{$1500\text{\,}\mathrm{m}\text{/}$\\}\times($300\text{\,}\mathrm{m}\text{/}$,$500\text{\,}\mathrm{m}\text{/}$),$
$\displaystyle\partial\Omega_{in,1}=\\{$-500\text{\,}\mathrm{m}\text{/}$\\}\times($1200\text{\,}\mathrm{m}\text{/}$,$1500\text{\,}\mathrm{m}\text{/}$)\times($300\text{\,}\mathrm{m}\text{/}$,$500\text{\,}\mathrm{m}\text{/}$),$
$\displaystyle\partial\Omega_{out,0}=\\{$-500\text{\,}\mathrm{m}\text{/}$\\}\times($100\text{\,}\mathrm{m}\text{/}$,$400\text{\,}\mathrm{m}\text{/}$)\times($-100\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$),$
$\displaystyle\partial\Omega_{out,1}=\\{$350\text{\,}\mathrm{m}\text{/}$\\}\times($100\text{\,}\mathrm{m}\text{/}$,$400\text{\,}\mathrm{m}\text{/}$)\times($-100\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$).$

The boundary conditions for flow are: homogeneous Dirichlet conditions on
$\partial\Omega_{out}$, uniform unit inflow on $\partial\Omega_{in}$, so that
$\int_{\partial\Omega_{in}}\bm{u}_{3}\cdot\bm{n}dS=$-1.2\text{\times}{10}^{5}\text{\,}{\mathrm{m}}^{3}\text{/}\mathrm{s}$$,
and homogeneous Neumann conditions on $\partial\Omega_{N}$. For the transport
problem, we consider homogeneous initial conditions, with a unit concentration
as boundary condition at $\partial\Omega_{in}$. The overall simulation time is
$5\text{\times}{10}^{3}\text{\,}\mathrm{s}\text{/}$.

##### Parameters

Matrix hydraulic conductivity $K_{3}$ | $\bm{I}$ | $\mathrm{m}\text{/}\mathrm{s}$
---|---|---
Fracture effective tangential hydraulic conductivity $K_{2}$ |  $1\text{\times}{10}^{2}$$\bm{I}$ | ${\mathrm{m}}^{2}\text{/}\mathrm{s}$
Fracture effective normal hydraulic conductivity $\kappa_{2}$ | $2\text{\times}{10}^{6}$ | $\text{/}\mathrm{s}$
Intersection effective tangential hydraulic conductivity $K_{1}$ | $1$ | ${\mathrm{m}}^{3}\text{/}\mathrm{s}$
Intersection effective normal hydraulic conductivity $\kappa_{1}$ | $2\text{\times}{10}^{4}$ | $\mathrm{m}\text{/}\mathrm{s}$
Matrix porosity $\phi_{3}$ | $2\text{\times}{10}^{-1}$ |
Fracture porosity $\phi_{2}$ | $2\text{\times}{10}^{-1}$ |
Intersection porosity $\phi_{1}$ | $2\text{\times}{10}^{-1}$ |
Fracture cross-sectional length $\epsilon_{2}$ | $1\text{\times}{10}^{-2}$ | $\mathrm{m}\text{/}$
Intersection cross-sectional area $\epsilon_{1}$ | $1\text{\times}{10}^{-4}$ | ${\mathrm{m}}^{2}\text{/}$

##### Folder in Git repository

field

##### Results

We consider a time step size of $\Delta t=$50\text{\,}\mathrm{s}\text{/}$$.
The results to be reported are:

  1. 1.

The number of cells and degrees of freedom according to the guidelines
described in Subsection 5.4 point 3 in a file named results.csv. In addition,
the total outflow over $\partial\Omega_{out,0}$ and $\partial\Omega_{out,1}$
as columns seven and eight;

  2. 2.

the average hydraulic head ${h}_{3}$ at the inlet, $\partial\Omega_{in}$,
reported as the seventh column of results.csv;

  3. 3.

the hydraulic head ${h}_{3}$ in the matrix along
$\left($350\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$,$-100\text{\,}\mathrm{m}\text{/}$\right)$-$\left($-500\text{\,}\mathrm{m}\text{/}$,$1500\text{\,}\mathrm{m}\text{/}$,$500\text{\,}\mathrm{m}\text{/}$\right)$;

  4. 4.

the hydraulic head ${h}_{3}$ in the matrix along
$\left($-500\text{\,}\mathrm{m}\text{/}$,$100\text{\,}\mathrm{m}\text{/}$,$-100\text{\,}\mathrm{m}\text{/}$\right)$-$\left($350\text{\,}\mathrm{m}\text{/}$,$1500\text{\,}\mathrm{m}\text{/}$,$500\text{\,}\mathrm{m}\text{/}$\right)$;

  5. 5.

the mean tracer concentrations throughout the simulation time in all
fractures.

In files called dol_line_${LINE}.csv report the results from point 3 and 4,
respectively for ${LINE} = 0 and ${LINE} = 1. For point 5, report in a file
dot.csv, with time in the first column, followed by the 52 mean concentration
values.

##  5 How to participate

We welcome all computational modelers to participate in the proposed benchmark
study. Three different stages of participation can be identified: 1\.
registration, 2. simulation and synchronization, 3. publication. In the
following, we will describe each stage in more detail.

###  5.1 Registration

If an interested researcher would like to participate in the benchmark study,
he/she should fill in the Google form provided at
https://goo.gl/forms/wCaPLipPI1d5Qz1j1. One form should be submitted for each
discretization method. The form asks for personal information from the
participant(s) as well as for a brief description of the discretization
method. After the application has been approved, the researcher will be asked
to sign a participation agreement that states that the participant should not
publish the results of other people before we do collectively. Once that is
signed, the researcher will be granted developer access to the Git repository
enabling him/her to download the required input data for the benchmark cases,
to view the simulation results of the other participants as well as to upload
his/her own results.

###  5.2 Simulation and synchronization

In this stage, the participant is expected to perform the computations and
convert his/her results to the requested format. In particular, he/she should
try the proposed scripts to evaluate the data as indicated in each case
description. The deadline for uploading the results to the Git repository is
31 January 2019. The participant should follow the GitLab feature branch
workflow that consists of contributing additions via a GitLab merge request of
a corresponding branch, analogous to the GitHub standard fork and pull request
workflow, see the details below in Subsection 5.4. The uploaded results will
be visible to all participants. This will initiate a synchronization phase in
which we will discuss and potentially recompute/adjust the results. In
particular, a minisymposium at the SIAM GS 19 is planned for presenting the
benchmark cases, the participants and the participating schemes as well as for
discussing the results. The synchronization phase is supposed to be finished
by uploading the final results before 15 May 2019.

###  5.3 Publication

As a final outcome of this benchmark study, we aim for a joint publication in
a peer-reviewed top-ranked journal. We would expect from each participant that
he/she helps us in writing and evaluating the manuscript before initial
submission as well as during the reviewing process. We would like to stress
that co-authorship can only be granted for at most two persons per
participating discretization method.

###  5.4 Details on reporting results

All data requested will be collected in the Git repository
{https://git.iws.uni-stuttgart.de/benchmarks/fracture-flow-3d.git}. To ensure
a unified and streamlined work flow, please observe the following structure
and formatting of the uploaded files.

The repository contains one folder for each of the four test cases. In the
results folder of these, each method will be assigned a folder with the format

$ACRONYM_INSTITUTION/$ACRONYM_NUMERICAL_SCHEME/.

No white spaces are allowed, please be consistent with upper and lower case.
As an example, results for the first test case obtained using the University
of Bergen implementation of the Mixed Virtual Element method are stored in
single/results/UiB/MVEM. Similarly, the string
$ACRONYM_INSTITUTION-$ACRONYM_NUMERICAL_SCHEME (e.g. UiB-MVEM) constitutes an
ID to be used in the handling and labeling of the results.

In the methods’ folder, three categories of files should be placed, for which
we here specify the format:

  1. 1.

Plot over line: A csv file whose first column contains the arc length and the
second column contains the corresponding values for either $c$ or ${h}$,
sampled at 2k equidistant points along the specified line. The delimiter is a
comma and there is no header for the file, an example is

        0.0,3.36948657742311
    0.5,5.38490345323433
    1.0,8.34820934803293

Generally, the file is called dol (data over line). If several refinement
levels are reported for the same test case, the files will be distinguished by
adding an appendix to the file name, see individual test case descriptions.
The plot over line can be easily done with the “Plot Over Line” filter of
ParaView, but a rearrangement of the output may be required.

  2. 2.

Plot versus time: A csv file named dot (data over time) whose first column
contains the time and the subsequent columns contain the values specified in
each test case. The same specifications apply as in the previous case.

  3. 3.

Cell numbers and other metrics: A file named results.csv where the first four
columns contain the total number of cells in $\Omega_{d},\ d=0,\dots,3$. The
next two columns contain the number of degrees of freedom dof and the number
of non-zero elements in the matrix nnz, respectively. These are followed by
other metrics in the order specified in the test case descriptions. If results
are reported for several refinement levels, the coarsest refinement level
corresponds to the first row etc. Again, the delimiter is a comma and there is
no header for the file.

Scripts for generating the plots are to be found in the scripts folders in
each of the test case folders. The scripts are named pol.py and pot.py and
contain comments on how to add new data. Once new results have been generated,
first update the Git repository. Add corresponding lines to the plot scripts
and upload the scripts and commit the corresponding data files to a newly
created feature branch. After pushing to the remote server, you can create a
corresponding merge request that will be reviewed and eventually approved by
the repository maintainers. Please verify that the scripts are working before
pushing to the remote Git repository. More details can be found at
https://docs.gitlab.com/ee/workflow/workflow.html and
https://docs.gitlab.com/ee/gitlab-basics/add-merge-request.html.

##  6 Appendix

###  6.1 Additional information for the test case of Subsection 4.2

We report the coordinates for the four corners of the fractures

$ID$ | $x_{0}$ | $y_{0}$ | $z_{0}$ | $x_{1}$ | $y_{1}$ | $z_{1}$ | $x_{2}$ | $y_{2}$ | $z_{2}$ | $x_{3}$ | $y_{3}$ | $z_{3}$
---|---|---|---|---|---|---|---|---|---|---|---|---
0 | 0.5 | 0 | 0 | 0.5 | 1 | 0 | 0.5 | 1 | 1 | 0.5 | 0 | 1
1 | 0 | 0.5 | 0 | 1 | 0.5 | 0 | 1 | 0.5 | 1 | 0 | 0.5 | 1
2 | 0 | 0 | 0.5 | 1 | 0 | 0.5 | 1 | 1 | 0.5 | 0 | 1 | 0.5
3 | 0.75 | 0.5 | 0.5 | 0.75 | 1.0 | 0.5 | 0.75 | 1.0 | 1.0 | 0.75 | 0.5 | 1.0
4 | 0.5 | 0.5 | 0.75 | 1.0 | 0.5 | 0.75 | 1.0 | 1.0 | 0.75 | 0.5 | 1.0 | 0.75
5 | 0.5 | 0.75 | 0.5 | 1.0 | 0.75 | 0.5 | 1.0 | 0.75 | 1.0 | 0.5 | 0.75 | 1.0
6 | 0.50 | 0.625 | 0.50 | 0.75 | 0.625 | 0.50 | 0.75 | 0.625 | 0.75 | 0.50 | 0.625 | 0.75
7 | 0.625 | 0.50 | 0.50 | 0.625 | 0.75 | 0.50 | 0.625 | 0.75 | 0.75 | 0.625 | 0.50 | 0.75
8 | 0.50 | 0.50 | 0.625 | 0.75 | 0.50 | 0.625 | 0.75 | 0.75 | 0.625 | 0.50 | 0.75 | 0.625

We report the $region\\_id$

$\displaystyle region\\_id(x<0.5\cap y<0.5\cap z<0.5)=0$ $\displaystyle
region\\_id(x>0.5\cap y<0.5\cap z<0.5)=1$ $\displaystyle region\\_id(x<0.5\cap
y>0.5\cap z<0.5)=2$ $\displaystyle region\\_id(x>0.5\cap y>0.5\cap z<0.5)=3$
$\displaystyle region\\_id(x<0.5\cap y<0.5\cap z>0.5)=4$ $\displaystyle
region\\_id(x>0.5\cap y<0.5\cap z>0.5)=5$ $\displaystyle region\\_id(x<0.5\cap
y>0.5\cap z>0.5)=6$ $\displaystyle region\\_id(x>0.75\cap y>0.75\cap
z>0.75)=7$ $\displaystyle region\\_id(x>0.75\cap y>0.5\cap y<0.75\cap
z>0.75)=8$ $\displaystyle region\\_id(x>0.5\cap x<0.75\cap y>0.75\cap
z>0.75)=9$ $\displaystyle region\\_id(x>0.5\cap x<0.75\cap y>0.5\cap
y<0.75\cap z>0.75)=10$ $\displaystyle region\\_id(x>0.75\cap y>0.75\cap
z>0.5\cap z<0.75)=11$ $\displaystyle region\\_id(x>0.75\cap y>0.5\cap
y<0.75\cap z>0.5\cap z<0.75)=12$ $\displaystyle region\\_id(x>0.5\cap
x<0.75\cap y>0.75\cap z>0.5\cap z<0.75)=13$ $\displaystyle
region\\_id(x>0.5\cap x<0.625\cap y>0.5\cap y<0.625\cap z>0.5\cap z<0.625)=14$
$\displaystyle region\\_id(x>0.625\cap x<0.75\cap y>0.5\cap y<0.625\cap
z>0.5\cap z<0.625)=15$ $\displaystyle region\\_id(x>0.5\cap x<0.625\cap
y>0.625\cap y<0.75\cap z>0.5\cap z<0.625)=16$ $\displaystyle
region\\_id(x>0.625\cap x<0.75\cap y>0.625\cap y<0.75\cap z>0.5\cap
z<0.625)=17$ $\displaystyle region\\_id(x>0.5\cap x<0.625\cap y>0.5\cap
y<0.625\cap z>0.625\cap z<0.75)=18$ $\displaystyle region\\_id(x>0.625\cap
x<0.75\cap y>0.5\cap y<0.625\cap z>0.625\cap z<0.75)=19$ $\displaystyle
region\\_id(x>0.5\cap x<0.625\cap y>0.625\cap y<0.75\cap z>0.625\cap
z<0.75)=20$ $\displaystyle region\\_id(x>0.625\cap x<0.75\cap y>0.625\cap
y<0.75\cap z>0.625\cap z<0.75)=21$

###  6.2 Additional information for the test case of Subsection 4.3

We report the coordinates for the four corners of the fractures

$ID$ | $x_{0}$ | $y_{0}$ | $z_{0}$ | $x_{1}$ | $y_{1}$ | $z_{1}$ | $x_{2}$ | $y_{2}$ | $z_{2}$ | $x_{3}$ | $y_{3}$ | $z_{3}$
---|---|---|---|---|---|---|---|---|---|---|---|---
0 | 0.05 | 0.25 | 0.5 | 0.95 | 0.25 | 0.5 | 0.95 | 2 | 0.5 | 0.05 | 2 | 0.5
1 | 0.5 | 0.05 | 0.05 | 0.5 | 0.05 | 0.95 | 0.5 | 0.3 | 0.95 | 0.5 | 0.3 | 0.05
2 | 0.05 | 1 | 0.5 | 0.95 | 1 | 0.5 | 0.95 | 2.2 | 0.85 | 0.05 | 2.2 | 0.85
3 | 0.05 | 1 | 0.48 | 0.95 | 1.0 | 0.48 | 0.95 | 2.2 | 0.14 | 0.05 | 2.2 | .14
4 | 0.17 | 1.9 | 0.7 | 0.17 | 1.9 | 0.3 | 0.23 | 2.2 | 0.3 | 0.23 | 2.2 | 0.7
5 | 0.23 | 1.9 | 0.7 | 0.23 | 1.9 | 0.3 | 0.17 | 2.2 | 0.3 | 0.17 | 2.2 | 0.7
6 | 0.77 | 1.9 | 0.7 | 0.77 | 1.9 | 0.3 | 0.77 | 2.2 | 0.3 | 0.77 | 2.2 | 0.7
7 | 0.83 | 1.9 | 0.7 | 0.83 | 1.9 | 0.3 | 0.83 | 2.2 | 0.3 | 0.83 | 2.2 | 0.7

## References

  * [1] Flemisch, B., Berre, I., Boon, W., Fumagalli, A., Schwenck, N., Scotti, A., Stefansson, I. & Tatomir, A.  Benchmarks for single-phase flow in fractured porous media.  _Advances in Water Resources_ 111, 239–258 (2018).  10.1016/j.advwatres.2017.10.036.
  * [2] Raviart, P.-A. & Thomas, J.-M.  A mixed finite element method for second order elliptic problems.  _Lecture Notes in Mathematics_ 606, 292–315 (1977).
  * [3] Brezzi, F. & Fortin, M.  _Mixed and Hybrid Finite Element Methods_ , vol. 15 of _Computational Mathematics_ (Springer Verlag, Berlin, 1991).
  * [4] Roberts, J. E. & Thomas, J.-M.  Mixed and hybrid methods.  In _Handbook of numerical analysis, Vol. II_ , Handb. Numer. Anal., II, 523–639 (North-Holland, Amsterdam, 1991).  10.1016/S1570-8659(05)80041-9.
  * [5] Ern, A. & Guermond, J.-L.  _Theory and Practice of Finite Elements_.  Applied mathematical sciences (Springer, 2004).  URL http://books.google.it/books?id=CCjm79FbJbcC.
  * [6] Alboin, C., Jaffré, J., Roberts, J. E., Wang, X. & Serres, C.  _Domain decomposition for some transmission problems in flow in porous media_ , vol. 552 of _Lecture Notes in Phys._ , 22–34 (Springer, Berlin, 2000).  10.1007/3-540-45467-5_2.
  * [7] Faille, I., Flauraud, E., Nataf, F., Pégaz-Fiornet, S., Schneider, F. & Willien, F.  A New Fault Model in Geological Basin Modelling. Application of Finite Volume Scheme and Domain Decomposition Methods.  In _Finite volumes for complex applications, III (Porquerolles, 2002)_ , 529–536 (Hermes Sci. Publ., Paris, 2002).
  * [8] Angot, P.  A model of fracture for elliptic problems with flux and solution jumps.  _Comptes Rendus Mathematique_ 337, 425–430 (2003).  10.1016/S1631-073X(03)00300-5.
  * [9] Martin, V., Jaffré, J. & Roberts, J. E.  Modeling Fractures and Barriers as Interfaces for Flow in Porous Media.  _SIAM J. Sci. Comput._ 26, 1667–1691 (2005).  10.1137/S1064827503429363.
  * [10] D’Angelo, C. & Scotti, A.  A mixed finite element method for Darcy flow in fractured porous media with non-matching grids.  _Mathematical Modelling and Numerical Analysis_ 46, 465–489 (2012).  10.1051/m2an/2011148.
  * [11] Fumagalli, A. & Scotti, A.  An Efficient XFEM Approximation of Darcy Flows in Arbitrarily Fractured Porous Media.  _Oil and Gas Sciences and Technologies - Revue d’IFP Energies Nouvelles_ 69, 555–564 (2014).  10.2516/ogst/2013192.
  * [12] Schwenck, N., Flemisch, B., Helmig, R. & Wohlmuth, B.  Dimensionally reduced flow models in fractured porous media: crossings and boundaries.  _Computational Geosciences_ 19, 1219–1230 (2015).  10.1007/s10596-015-9536-1.
  * [13] Flemisch, B., Fumagalli, A. & Scotti, A.  _A Review of the XFEM-Based Approximation of Flow in Fractured Porous Media_ , vol. 12 of _SEMA SIMAI Springer Series_ , chap. Advances in Discretization Methods, 47–76 (Springer International Publishing, Cham, 2016).  10.1007/978-3-319-41246-7_3.
  * [14] Boon, W. M., Nordbotten, J. M. & Yotov, I.  Robust discretization of flow in fractured porous media.  _SIAM Journal on Numerical Analysis_ 56, 2203–2233 (2018).  10.1137/17M1139102.
  * [15] Fumagalli, A. & Keilegavlen, E.  Dual virtual element methods for discrete fracture matrix models.  Tech. Rep., arXiv:1711.01818 [math.NA] (2017).  URL https://arxiv.org/abs/1711.01818.
  * [16] Nordbotten, J. M., Boon, W., Fumagalli, A. & Keilegavlen, E.  Unified approach to discretization of flow in fractured porous media.  _Computational Geosciences_ (2018).  URL https://arxiv.org/abs/1802.05961.  Accepted.
  * [17] Zielke, W., Helmig, R., Krohn, K., Shao, H. & Wollrath, J.  Discrete modelling of transport processes in fractured porous rock.  In _7th ISRM Congress_ , 57–60 (1991).  URL https://www.onepetro.org/conference-paper/ISRM-7CONGRESS-1991-013.
  * [18] Barlag, C., Hinkelmann, R., Helmig, R. & Zielke, W.  Adaptive methods for modelling transport processes in fractured subsurface systems.  In _3rd International Conference on Hydroscience and Engineering, Cottbus_ (1998).
  * [19] Geiger, S., Dentz, M. & Neuweiler, I.  A novel multi-rate dual-porosity model for improved simulation of fractured and multiporosity reservoirs.  _SPE Journal_ 18, 670–684 (2013).  10.2118/148130-PA.