[2]Alexander Litvinenko

[1]CSE Department, KAUST, Thuwal-Jeddah, Saudi Arabia

2]Department of Mathematics, RWTH Aachen, Aachen, Germany

#  Estimation of uncertainties in the density driven flow in fractured porous
media using MLMC

Dmitry Logashenko<EMAIL_ADDRESS><EMAIL_ADDRESS>Raul Tempone<EMAIL_ADDRESS>Gabriel Wittum
<EMAIL_ADDRESS>*  [

###### Abstract

We use the Multi Level Monte Carlo method to estimate uncertainties in a
Henry-like salt water intrusion problem with a fracture. The flow is induced
by the variation of the density of the fluid phase, which depends on the mass
fraction of salt. We assume that the fracture has a known fixed location but
an uncertain aperture. Other input uncertainties are the porosity and
permeability fields and the recharge. In our setting, porosity and
permeability vary spatially and recharge is time-dependent. For each
realisation of these uncertain parameters, the evolution of the mass fraction
and pressure fields is modelled by a system of non-linear and time-dependent
PDEs with a jump of the solution at the fracture. The uncertainties propagate
into the distribution of the salt concentration, which is an important
characteristic of the quality of water resources. We show that the multilevel
Monte Carlo (MLMC) method is able to reduce the overall computational cost
compared to classical Monte Carlo methods. This is achieved by balancing
discretisation and statistical errors. Multiple scenarios are evaluated at
different spatial and temporal mesh levels. The deterministic solver ug4 is
run in parallel to calculate all stochastic scenarios.

###### keywords:

uncertainty quantification, MLMC, density-driven flow, aquifer, porosity,
permeability, fracture, multigrid methods

###### Contents

  1. 1 Introduction
  2. 2 Modeling and Problem Settings
    1. 2.1 Governing equations for the flow
    2. 2.2 Model problem settings
    3. 2.3 Stochastic modeling of porosity, permeability and fracture width
  3. 3 Numerical Methods
    1. 3.1 Numerical methods for the deterministic problem
    2. 3.2 MLMC Algorithm
  4. 4 Numerical Experiments
    1. 4.1 The mean value and variance
    2. 4.2 Estimation of convergence rates $\alpha$ and $\beta$ for various QoIs
    3. 4.3 Comparison of MC and MLMC methods
  5. 5 Discussion and Conclusion

##  1 Introduction

An important and challenging problem in hydrogeology is the modeling of
densiity-driven subsurface flow and salt transport in fractured aquifers.
Essential difficulties in these settings are complicated and parameter-
sensitive flow pattern due to the fractures as well as a number of the
uncertain parameters. Uncertainties in the fracture geometry, porosity,
permeability, recharge can strongly influence the evolution of the salt
concentration in the subsurface (cf. [1, 3, 4, 33, 44]). In a real reservoir
the number of the fractures, their location and their aperture are uncertain.
One possible solution approach would be to model all these uncertainties,
build a PCE-based surrogate, perform a sensitive analysis (see [5, 45]) to
understand which uncertainties are very critical and which are not, and
estimate how they contribute to the solution. This is a very non-trivial
approach. Instead, we apply the Multi Level Monte Carlo method (MLMC) method.
Another reason for the usage of MLMC method is that the deterministic problem
is resolved via the multigrid method. So, there is a natural idea to try to
couple the MLMC and multigrid methods.

An essential problem in water production from coastal aquifers is their
salinization due to the intrusion of the seawater. Salt concentration in the
drinking water should be kept below the prescribed limits, and therefore its
forecast ois important for water resources management in coastal regions.

Fractures introduce heterogeneity in the aquifer that affects fluid flow and
makes it difficult to predict its main characteristics. A challenge is that
fractures can get sealed over time, reducing their permeability and affecting
fluid flow. In summary, adding fractures to subsurface aquifers presents
challenges related to complexity, sealing, heterogeneity and environmental
concerns. However, with careful planning and management, fractures can also
provide opportunities for increased production and improved reservoir
performance.

As a model problem, we consider the benchmark from [20] which is a
generalization of the Henry problem first introduced in [25]. The Henry
problem became a benchmark for numerical solvers for the density-driven
groundwater flow (see [50, 43, 42, 14]. In [39], the authors use the
generalised polynomial chaos expansion approximation to investigate how
incomplete knowledge of system properties affects the assessment of global
quantities. In particular, they estimated the propagation of input
uncertainties into a few dimensionless scalar parameters.

In [20, 21, 38], the authors modeled and computed the density-driven flow in
fractured porous media with deterministic settings. The fractures considered
were represented by a $(d-1)$-dimensional manifolds in a $d$-dimensional
domain. They developed a special numerical technique based on the finite
volume method and implemented it in the ug3 library. In our present work, we
apply the same model, the same discretization and solvers to compute the
multiple uncertain scenarios.

In [21], the validity of the low-dimensional representation of the fractures
has been addressed. The authors considered both $d$ and $(d-1)$-dimensional
fractures. They have shown numerically that the $(d-1)$-dimensional
representation can be a valid alternative to the full-dimensional resolution
of all the subdomains provided that the fracture apertures are in some
reasonable range. All the settings in this work were deterministic.

Other methods based on $(d-1)$-dimensional representation of fractures were
considered in [1, 3, 4, 33, 44, 41, 19].

Many techniques can be used to quantify uncertainties. One classical method is
Monte Carlo (MC) sampling. It has a well-known disadvantage — a slow
convergence $\mathcal{O}(\frac{1}{\sqrt{N}})$. Other relative recent
techniques such as surrogate models and stochastic collocation may require a
few hundred time-consuming simulations and assume a certain smoothness of the
quantity of interest (QoI).

Another class of methods is the class of perturbation methods [12]. The idea
is to decompose the QoI with respect to random parameters in a Taylor series.
The higher order terms can be neglected for small perturbations, simplifying
the analysis and numerics. These methods assume that the random perturbations
are small. For larger perturbations, these methods usually do not work.

There are a number of studies where authors have modeled uncertainties in
reservoirs (see [7, 48]). The link between stochastic methods and
hydrogeological applications was made in [6]. Further the authors made
recommendations to water suppliers in Germany about optimisation and risk
assessment. The basics of stochastic hydrogeology and an overview of
stochastic tools and uncertainty management are described in [40].

The review [46] deals with hydrogeological applications of recent advances in
uncertainty quantification, probabilistic risk assessment and decision making
under uncertainty.

This work has the following structure. Section 2 describes the model. Section
3 contains description of numerical methods. In particular, the well-known
MLMC method is reviewed in Section 3.2. Section 4 details the numerical
results, which include the numerical analysis of the model problem, the
computation of different statistics, the performance of the MLMC method, and
the performance of the parallel solver with uncertain coefficients. Finally,
we conclude this work with a discussion in Section 5.

##  2 Modeling and Problem Settings

We consider the density driven flow of the liquid phase (salt solution in
water) with variable mass fraction of the salt in a fractured aquifer. The
fractures are geological formations whose size in one of the geometrical
directions is negligible in comparison with the scales of the domain whereas
the permeability of the filling material is much higher than that of the
surrounding medium. This causes essential problems when these fractures are
considered as full-dimensional subdomains in the simulations: Specific, very
fine grids needed for the resolution of this sort of the objects as well as
the strong jump of the model parameters reduces the efficiency of the
numerical methods significantly. For this reason, several other approaches
have been developed where the fractures are represented by low-dimensional
manifolds and the fracture aperture (width) is a parameter (cf. [1, 4, 16, 13,
32, 28, 38]). Whereas it has been shown that for the density-driven flow,
these models provide a proper approximation only for a certain range of the
fracture width ([21]), their application is considered to be reasonable in
many situations. A particular advantage of such the models for the sampling
methods is the possibility to change the fracture width without remeshing of
the entire domain. In this paper, we assume that the porous matrix is immobile
and the fractures do not change their positions.

###  2.1 Governing equations for the flow

In this section, we briefly describe the PDEs and refer to [36, 38] for
details. From the modeling point of view, we consider the porous medium
$\mathscr{M}\subset\mathcal{D}$ surrounding the fracture with the porosity
$\phi_{m}:\mathscr{M}\to(0,1)$ and the permeability
${\mathbf{K}}_{m}:\mathscr{M}\to\mathbb{R}^{d\times d}$. The fluid phase in
this subdomain is characterized by salt mass fraction
$c_{m}(t,\mathbf{x}):[0,+\infty)\times\mathscr{M}\to[0,1]$ and pressure
$p_{m}(t,\mathbf{x}):[0,+\infty)\times\mathscr{M}\to\mathbb{R}$. So, the flow
in $\mathscr{M}$ obeys the standard mass conservation laws for the entire
liquid phase and the salt

$\left.\begin{array}[]{l}\partial_{t}(\phi_{m}\rho_{m})+\nabla\cdot(\rho_{m}{\mathbf{q}}_{m})=0\\\
\partial_{t}(\phi_{m}\rho_{m}c_{m})+\nabla\cdot(\rho_{m}c_{m}{\mathbf{q}}_{m}-\rho_{m}{\mathbf{D}}_{m}\nabla
c_{m})=0\end{array}\right\\}\quad x\in\mathscr{M}$ (1)

with the Darcy’s law for the velocity:

$\displaystyle{\mathbf{q}}_{m}=-\frac{{\mathbf{K}}_{m}}{\mu_{m}}(\nabla
p_{m}-\rho_{m}{\mathbf{g}}),\qquad x\in\mathscr{M},$ (2)

where $\rho_{m}=\rho(c_{m})$ and $\mu_{m}=\mu(c_{m})$ indicate the density and
the viscosity of the liquid phase,
${\mathbf{D}}_{m}(t,\mathbf{x}):[0,+\infty)\times\mathscr{M}\to\mathbb{R}^{d\times
d}$ denotes the molecular diffusion and mechanical dispersion tensor.

The fracture is assumed to be filled with a porous medium, too. In particular,
we assume the validity of the Darcy’s law inside the fracture. Without loss of
generality, we consider in this work only one fracture, see [36, 38] for the
generalization. The fracture is represented by a surface
$\mathscr{S}\subset\mathcal{D}$, $\mathscr{M}\cup\mathscr{S}=\mathcal{D}$,
$\mathscr{M}\cap\mathscr{S}=\emptyset$. However, we distinguish between the
two sides of the fracture, $\mathscr{S}^{(1)}$ and $\mathscr{S}^{(2)}$, that
geometrically coincide with $\mathscr{S}$ but are virtually interfaces of the
fracture with $\mathscr{M}$ and have opposite normals $\mathbf{n}^{(1)}$ and
$\mathbf{n}^{(2)}$.

For the fluid phase inside the fracture, let
$c_{f}:[0,+\infty)\times\mathscr{S}\to[0,1]$ and
$p_{f}:[0,+\infty)\times\mathscr{S}\to\mathbb{R}$ be the salt mass fraction
and the pressure averaged along the vertical, cf. [21]. Note that in general,
the values of $c_{f}$ and $p_{f}$ are not equal to the limits of $c_{m}$ and
$p_{m}$ at $\mathscr{S}$. Furthermore, the latter limits may be different on
$\mathscr{S}^{(1)}$ and $\mathscr{S}^{(2)}$. We denote them by $c_{m}^{(1)}$
and $c_{m}^{(2)}$, as well as $p_{m}^{(1)}$ and $p_{m}^{(2)}$ as functions on
$\mathscr{S}^{(1)}$ and $\mathscr{S}^{(2)}$, respectively. Two sets of the
equations are stated for the fracture: The laws that describe the flow along
the fracture and the ones modeling the mass exchange through the interfaces
$\mathscr{S}^{(1)}$ and $\mathscr{S}^{(2)}$.

For the lateral flow in the fracture, we impose the analogues of the mass
conservation laws (1):

$\left.\begin{array}[]{l}\partial_{t}(\phi_{f}\epsilon\rho_{f})+\nabla^{\mathscr{S}}\cdot(\epsilon\rho_{f}{\mathbf{q}}_{f})+Q_{fn}^{(1)}+Q_{fn}^{(2)}=0\\\
\partial_{t}(\phi_{f}\epsilon\rho_{f}c_{f})+\nabla^{\mathscr{S}}\cdot(\epsilon\rho_{f}c_{f}{\mathbf{q}}_{f}-\epsilon\rho_{f}D_{f}\nabla^{\mathscr{S}}c_{f})+P_{fn}^{(1)}+P_{fn}^{(2)}=0\end{array}\right\\}\quad
x\in\mathscr{S},$ (3)

where $\nabla^{\mathscr{S}}$ denotes the differential operators on the
manifold $\mathscr{S}$ and $\epsilon$ the fracture width. The Darcy velocity
along the fracture is

$\displaystyle{\mathbf{q}}_{f}=-\frac{K_{f}}{\mu_{f}}(\nabla^{\mathscr{S}}p_{f}-\rho_{f}{\mathbf{g}}),\qquad
x\in\mathscr{S}.$ (4)

In (3–4), $\rho_{f}=\rho(c_{f})$ and $\mu_{f}=\mu(c_{f})$ are the density and
the viscosity of the fluid phase in the fracture. Furthermore,
$\phi_{f}:\mathscr{S}\to[0,1]$ denotes the porosity associated with the
material in the fracture, $K_{f}:\mathscr{S}\to\mathbb{R}$ its lateral
permeability and $D_{f}(t,x):[0,+\infty)\times\mathscr{S}\to\mathbb{R}$ the
lateral diffusion and dispersion tensor in the fracture.

The terms $Q_{fn}^{(k)}$ and $P_{fn}^{(k)}$, $k\in\\{1,2\\}$, are the mass
fluxes through the faces $\mathscr{S}^{(k)}$ of the fracture, i.e. the normal
fluxes. They are defined as

$\left.\begin{array}[]{rl}Q_{fn}^{(k)}&:=\rho(c_{m}^{(k)})q_{fn}^{(k)}\\\
P_{fn}^{(k)}&:=\rho(c_{m}^{(k)})c^{(k)}_{\mathrm{upwind}}q_{fn}^{(k)}-\rho(c_{m}^{(k)})D_{fn}^{(k)}\dfrac{c_{m}^{(k)}-c_{f}}{\epsilon/2}\end{array}\right\\}\quad
x\in\mathscr{S}^{(k)},$ (5)

where $c^{(k)}_{\mathrm{upwind}}=c_{m}^{(k)}$ if $q_{fn}^{(k)}<0$ but
$c^{(k)}_{\mathrm{upwind}}=c_{f}$ if $q_{fn}^{(k)}\geq 0$ with

$q_{fn}^{(k)}:=-\frac{K_{fn}^{(k)}}{\mu(c_{f}^{(k)})}\left[\frac{p_{m}^{(k)}-p_{f}}{\epsilon/2}-(\rho(c_{f}^{(k)})-\rho_{f}){\mathbf{g}}\cdot\mathbf{n}^{(k)}\right],\qquad
x\in\mathscr{S}^{(k)}.$ (6)

In (6), $K_{fn}^{(k)}:\mathscr{S}^{(k)}\to\mathbb{R}$ is the normal
permeability and
$D_{fn}^{(k)}:[0,+\infty)\times\mathscr{S}^{(k)}\to\mathbb{R}$ is the normal
diffusion and dispersion coefficient of the fracture – bulk medium interface
$\mathscr{S}^{(k)}$.

Systems (1) and (3) are coupled by the continuity of the mass fluxes of the
interfaces $\mathscr{S}^{(k)}$, $k\in\\{1,2\\}$:

$\left.\begin{array}[]{rl}\rho(c_{m}^{(k)}){\mathbf{q}}_{m}\cdot\mathbf{n}^{(k)}&=Q_{fn}^{(k)}\\\
\rho(c_{m}^{(k)})c_{m}^{(k)}{\mathbf{q}}_{m}\cdot\mathbf{n}^{(k)}+\rho(c_{m}^{(k)}){\mathbf{D}}_{m}\nabla
c_{m}\cdot\mathbf{n}^{(k)}&=P_{fn}^{(k)}\end{array}\right\\},\qquad
x\in\mathscr{S}^{(k)}$ (7)

Furthermore, for the inner edges of fractures, we impose the no-flux boundary
conditions so that no mass exchange is possible there, cf. [20, 21].

Equations (1), (3) and (7) with the velocities (2), (4) and (6) as well as
particular specifications of the parameters $\rho$, $\mu$, $\phi_{m,f}$,
${\mathbf{K}}_{m}$, ${\mathbf{D}}_{m}$, $K_{f}$, $K_{fn}$, $D_{f}$ and
$D_{fn}$ form the complete model for the density-driven flow in the porous
medium. This system must be closed by boundary conditions for $c_{m}$,
$c_{f}$, $p_{m}$ and $p_{f}$ on $\partial\mathcal{D}$ as well as initial
conditions for $c_{m}$ and $c_{f}$ at $t=0$.

###  2.2 Model problem settings

For our numerical tests, we choose a simple but very illustrative two-
dimensional problem with one fracture, proposed in [20]. It extends the
setting of the Henry problem [25, 43] with a fracture. The aquifer is
represented by a rectangular domain $\mathcal{D}=[0,2]\times[-1,0]$
$[\mathrm{m}^{2}]$ completely saturated with the liquid phase, see Fig. 1
(left). The fracture is located near the right “sea side” where the heavy
salty water intrudes into the aquifer and cuts this boundary.

Figure 1: Left: Geometry, boundary conditions and the selected points for the
modified Henry problem with a fracture. Right: The flow pattern of this
problem; color: $c_{m}$ (red for $c_{m}=1$); black lines: streamlines of the
flow field.

The model presented in Section 2.1 must be closed by specification of boundary
conditions for $c_{m,f}$ and $p_{m,f}$, as well as initial conditions for
$c_{m,f}$. In this work, we follow the settings from [20] but add the
uncertainty as described in Section 2.3 below. In particular, for the initial
conditions, we set

$\left.c_{m}\right|_{t=0}=0,\qquad\left.c_{f}\right|_{t=0}=0.$ (8)

A scheme of the boundary conditions is shown in Fig. 1 (left). On the right
boundary, the Dirichlet conditions for $c_{m,f}$ and $p_{m,f}$ model the
seawater intrusion:

$\left.c_{m,f}\right|_{x=2}=1,\qquad\left.p_{m,f}\right|_{x=2}=-\rho_{1}gy.$
(9)

On the left boundary, the inflow (recharge) of fresh water is imposed:

$\left.c_{m}\right|_{x=0}=0,\qquad\left.\rho{\mathbf{q}}_{m}\cdot\mathbf{e}_{x}\right|_{x=0}=\hat{q}_{\mathrm{in}},$
(10)

where $\mathbf{e}_{x}=(1,0)^{\top}$, and $\hat{q}_{\mathrm{in}}$ is a
prescribed function of time, see Sect. 2.3. Note that the fracture does not
cut the left side of the domain, so no boundary conditions for $c_{f}$ and
$p_{f}$ are required there. For the upper and the lower boundaries of
$\mathcal{D}$, we impose no-flux boundary conditions for all the equations.

The flow and salt transport patterns appearing in this problem are discussed
in [20]. They are presented in Fig. 1 (right): The color codes the mass
fraction $c_{m}$ with red corresponding to $c_{m}=1$ and blue to $c_{m}=0$;
the lines are the streamlines of the flow. In the left part of $\mathcal{D}$,
the flow is induced by the pure water recharge through the left boundary. As
in the Henry problem, the salt water intruding through the lower part of the
right boundary turns up and is washed out through the upper part of the right
side. In particular, the much more permeable fracture serves as a pathway for
the strong wash-out that evacuates the salty water from the domain. For that,
the part of the domain below the fracture is mostly separated from the upper
part of the domain — as much as the capacity of the fracture allows. These
phenomena depend on the width and the permeability of the fracture. In this
work, we investigate how the fracture aperture influences the concentration.

In the case of time-independent boundary conditions, the mass fraction field
converges to a steady state (cf. [20]). However, before it, there is an
initial instationary phase taking considerable time when the maximum
uncertainty in the solution w.r.t. the parameters is attained. This phase can
be recognized even for the time-dependent recharge $\hat{q}_{\mathrm{in}}$
considered below.

In this paper, as for the Henry problem (see [50, 43, 42]), we set

$\displaystyle\rho(c)=\rho_{0}+(\rho_{1}-\rho_{0})c,\qquad\qquad\mu=\text{const},$
(11)

where $\rho_{0}$ is the density of the pure water and $\rho_{1}$ is the
density of the brine (both considered as constants). Besides that, we assume
that porous medium is isotropic:

$\displaystyle{\mathbf{K}}_{m}=K_{m}\mathbf{I},$ (12)

where $K_{m}:\mathscr{M}\to\mathbb{R}$. Furthermore, in this work, we neglect
the mechanical dispersion, so that

$\displaystyle{\mathbf{D}}_{m}=\phi_{m}D_{\mathrm{mol}}\mathbf{I},\qquad{\mathbf{D}}_{f}=\phi_{f}D_{\mathrm{mol}}$
(13)

where $D_{\mathrm{mol}}\in\mathbb{R}$ is the constant molecular diffusion
coefficient of the salt in the liquid phase. For the interface between the
bulk medium and the fracture, we assume

$\displaystyle K_{fn}=K_{m},\qquad D_{fn}=\phi_{m}D_{\mathrm{mol}},\qquad
x\in\mathscr{S}.$ (14)

Deterministic values of these parameters that are not varied in the scenarios
are presented in Table 1.

We investigate the mass fraction $c_{m}$ at 6 locations where the largest
variance is observed. This corresponds to the situation when the water quality
of wells reaching these points is controlled. The list of selected points
follows (see Fig. 1):

$\displaystyle\\{\mathbf{x}_{i}:=(x_{i},y_{i})_{i=1,\ldots,6}\\}=\\{$
$\displaystyle(1.1,-0.8),(1.2,-0.8),(1.3,-0.8),$
$\displaystyle(1.1,-0.9),(1.2,-0.9),(1.3,-0.9).\\}$ (15)

Whereas the other points can be considered, the reason to select these 6
points (and 6 small subdomains) is that not all other points are
“interesting”, i.e., not all points have significant variation in $c_{m}$.
MLMC reduces the variance, but if the initial variance is small, there is no
need to use MLMC.

Symbol Quantity Value Unit $D_{\mathrm{mol}}$ Deterministic diffusion
coefficient in the medium $18.8571\cdot 10^{-6}$ [$\mathrm{m}^{2}\textrm{
}\mathrm{s}^{-1}$] $\mathbf{g}$ Gravity $9.8$ [$\mathrm{m}\textrm{
}\mathrm{s}^{-2}$] $K_{m}$ Deterministic permeability of the medium
$1.019368\cdot 10^{-9}$ [$\mathrm{m}^{2}$] $K_{f}$ Deterministic permeability
of the fracture $1.019368\cdot 10^{-6}$ [$\mathrm{m}^{2}$] $\phi_{m}$
Deterministic porosity of the medium $0.35$ - $\phi_{f}$ Deterministic
porosity of the fracture $0.7$ - $\mu$ Viscosity $10^{-3}$
[$\mathrm{kg}\textrm{ }\mathrm{m}^{-1}\textrm{ }\mathrm{s}^{-1}$] $\rho^{pW}$
Density of water $1\cdot 10^{3}$ [$\mathrm{kg}\textrm{ }\mathrm{m}^{-3}$]
$\rho^{pB}$ Density of brine $1.025\cdot 10^{3}$ [$\mathrm{kg}\textrm{
}\mathrm{m}^{-3}$] $\epsilon$ Fracture width
$\mathbb{E}\left[\epsilon\right]=5.05\cdot 10^{-3}$ [$\mathrm{m}$]

Table 1: Simulation parameters for the model problem from Section 2.2

###  2.3 Stochastic modeling of porosity, permeability and fracture width

Stochastic modeling is used to simulate the behavior of a fractured reservoir
under uncertainty. The primary sources of uncertainty are the hydrogeological
properties of the porous medium — porosity ($\phi_{m,f}$) and permeability
($K_{m,f}$) fields of the porous matrix and the fracture aperture ($\epsilon$)
— as well as the freshwater recharge intensity $\hat{q}_{x}$. In this work we
do not consider the uncertainty in $\phi_{f}$ and $K_{f}$. The QoIs are
related to the mass fraction $c_{m}$, a function of $\phi_{m}$, $K_{m}$, and
$\hat{q}_{x}$. We model the uncertain $\phi_{m}$ using a random field and
assume $K_{m}$ to dependent on $\phi_{m}$:

$\displaystyle K_{m}=K_{m}(\phi_{m})\in\mathbb{R}.$ (16)

The distribution of $\phi(\mathbf{x},\bm{\xi})$, $\mathbf{x}\in\mathcal{D}$,
is determined by a set of stochastic parameters
$\bm{\xi}=(\xi_{1},\ldots,\xi_{M},...)$. Each component $\xi_{i}$ is a random
variable.

The dependence in (16) is specific for every material. We refer to [34, 35,
11] for a detailed discussion. In this work, we use a Kozeny-Carman-like law

$\displaystyle
K_{m}(\phi_{m})=\kappa_{KC}\cdot\dfrac{\phi_{m}^{3}}{1-\phi_{m}^{2}},$ (17)

where $\kappa_{KC}$ is a constant scaling factor. The recharge inflow flux is
kept constant across the left boundary but depends on the stochastic variable
$\hat{q}_{x}$. We also assume that it is independent of $\phi$ and
${\mathbf{K}}$.

The uncertain width of the fracture, the recharge, and the porosity are
modeled as follows, $\xi_{1},\xi_{2},\xi_{3}\in U[-1,1]$,

$\epsilon(\xi_{1})=0.01\cdot((1-0.01)\cdot\xi_{1}+(1+0.01))/2,$ (18)
$\hat{q}_{x}(t,\xi_{3})=3.3\cdot 10^{-6}\cdot(1+0.1\xi_{3})(1+0.1\sin(\pi
t/40)),$ (19) $\phi_{m}(x,y,\xi_{2})=0.35\cdot(1+0.02\cdot(\xi_{2}\cos(\pi
x/2)+\xi_{2}\sin(2\pi y)).$ (20)

The corresponding permeability is defined by (17).

Unknown porosity can be also modeled by a random field. One of the ways to
compute this random field is to make an assumption about the covariance matrix
or to estimate it from the available measurement data. Then an auxiliary
eigenvalue problem is solved and the truncated Karhunen-Loéve expansion (KLE)
is constructed [29]. Due to the high complexity of this approach, we implement
a simpler one. We use $L_{2}$ orthogonal functions ($\sin()$ and $\cos()$),
which mimic the $L_{2}$-orthogonal functions used in the KLE. These functions
are then multiplied by the uniform random variables $\bm{\xi}$ as in (20).
Typically one takes one of three options for $\bm{\xi}$: uniform, Gaussian or
log-normal.

##  3 Numerical Methods

###  3.1 Numerical methods for the deterministic problem

The system (1)-(3) is numerically solved in the domain
$\mathcal{D}\times[0,T]$. $\mathcal{D}$ is covered with an unstructured grid
$\mathcal{D}_{h}$ of triangles and quadrilaterals. We denote the
characteristic mesh size by $h$. We apply a special technique for resolving
the jump of the solution at the fracture, cf. [36] for details. The vertex-
centered finite-volume scheme is used for the discretization of the system
(1–4) in space (cf. [36, 38]). The number of degrees of freedom associated
with $\mathcal{D}_{h}$ is denoted by $n$. There are two degrees of freedom per
grid vertex in $\mathcal{D}_{h}$: one for the mass fraction $c_{m,f}$ and
another for the pressure $p$. At the fracture, several grid vertices share the
same geometrical position, cf. [36, 38]. We use the regular refinement to
create the grid hierarchy so that $h=\mathcal{O}(n^{-1/d})$. The implicit
Euler method with a time step $\tau$ is used for the time discretization of
(1–4). The number of the computed time steps is $r=T/\tau$. In this paper, we
focus on the mass fraction $c_{m}$ of the salt in the bulk porous medium. We
denote the part of the discrete solution of the model approximating $c_{m}$ by
$c_{m,h,\tau}$.

We use the full upwind for the convective terms (cf. [15]). Therefore, the
discretization error is of the first order w.r.t. the spatial mesh size $h$.
Furthermore, the Euler method provides the first-order discretization error
w.r.t. $\tau$. Thus, as $d=2$,

$\|c_{m}-c_{m,h,\tau}\|_{2}=\mathcal{O}(h+\tau)=\mathcal{O}(n^{-1/2}+r^{-1}),$
(21)

which is consistent with our numerical tests.

The implicit time-stepping scheme is unconditionally stable but requires
solution of the extensive nonlinear discretized algebraic system with $n$
unknowns in every time step. The Newton’s method is used for this. Linear
systems in the Newton iteration are solved using the BiCGStab method (cf. [2])
preconditioned with the geometric multigrid method (V-cycle, cf. [22]). In the
multigrid cycle, the ILUβ-smothers [23] and Gaussian elimination are used as
the coarse grid solver.

###  3.2 MLMC Algorithm

To reduce the total computing cost, we apply the MLMC method, which is a
natural idea because the deterministic solver uses a multigrid method (see
Section 3.1). The MLMC method gains efficiency by combining samples computed
on different grids in the hierarchy. A more in-depth description of these
techniques is found in [9, 10, 17, 18, 24, 47, 30].

Let $\bm{\xi}$ be a vector of random variables, and $g(\bm{\xi})$ the quantity
of interest. In this work, we tried different QoI $g$: solution $c_{m}$ at a
point $(t,\mathbf{x})$, integrals over a small subdomain $\Delta_{i}$ and over
the whole computational domain $\mathcal{D}$, see Sect. 4. The MLMC method
aims to approximate the expected value $\mathbb{E}\left[g\right]$ with an
optimal computational cost. MLMC constructs a telescoping sum, defined over a
sequence of spatial and temporal meshes, $\ell=0,\ldots,L$, as described next,
to achieve this goal. The QoI $g$, numerically evaluated on level $\ell$, is
denoted by $g_{h_{\ell},\tau_{\ell},\ell}$ or, for simplicity, by just
$g_{\ell}$, where $h_{\ell}$ and $\tau_{\ell}$ are the discretization steps in
space and time on level $\ell$.

Let $s_{0}$ denote the maximum computational cost of the evaluation of one
realization of $g_{0}$. Similarly, $s_{\ell}$ denotes the computing cost of
evaluating $g_{\ell}-g_{\ell-1}$. For simplicity, we assume that $s_{\ell}$
for $g_{\ell}-g_{\ell-1}$ is almost the same as $s_{\ell}$ for $g_{\ell}$.
Since the number of required iterations is variable, the cost of computing a
sample of $g_{\ell}-g_{\ell-1}$ may fluctuate for various realizations.

The MLMC method calculates
$\mathbb{E}\left[g_{L}\right]\approx\mathbb{E}\left[g\right]$ using the
following telescopic sum:

$\displaystyle\mathbb{E}\left[g_{L}\right]$
$\displaystyle=\mathbb{E}\left[g_{0}\right]+\sum_{\ell=1}^{L}\mathbb{E}\left[g_{\ell}-g_{\ell-1}\right]$
(22) $\displaystyle\approx
m_{0}^{-1}\sum_{i=1}^{m_{0}}g_{0}^{(0,i)}+\sum_{\ell=1}^{L}\left(m_{\ell}^{-1}\sum_{i=1}^{m_{\ell}}(g_{\ell}^{(\ell,i)}-g_{\ell-1}^{(\ell,i)})\right).\
$ (23)

In the above equation, level $\ell$ in the superscript $(\ell,i)$ indicates
that independent samples are used at each correction level.

As $\ell$ increases, the variance of $g_{\ell}-g_{\ell-1}$ decreases. Thus,
the total computational cost can be reduced by taking fewer samples on finer
meshes.

In the following numerical experiments
$n_{\ell}=4n_{\ell-1}=\ldots=4^{\ell}n_{0}=2^{d\ell}n_{0}$, $d=2$, and
$r_{\ell}=2r_{\ell-1}=\ldots=2^{\ell}r_{0}$. In the case of uniform,
equidistant mesh, we could also write similar formulas for step sizes:
$h_{\ell}=h_{\ell-1}\cdot 2^{-1}=h_{\ell-2}\cdot 2^{-2}=\ldots=h_{0}\cdot
2^{-\ell}$ and $\tau_{\ell}=\tau_{0}\cdot 2^{-\ell}$.

The average cost $s_{\ell}$ of generating one sample of $g_{\ell}$ (the cost
of one deterministic simulation for one random realization) is

$s_{\ell}=\mathcal{O}(n_{\ell}r_{\ell})=\mathcal{O}(4^{\ell}n_{0}\cdot
2^{\ell}r_{0})=\mathcal{O}(2^{2\ell}n_{0}\cdot
2^{\ell}r_{0})=\mathcal{O}(2^{\hat{d}\ell\gamma}n_{0}r_{0}),$ (24)

where $\hat{d}=d+1=3$ and $\gamma=1$.

######  Definition 1.

Let
$Y_{\ell}:=m_{\ell}^{-1}\sum_{i=1}^{m_{\ell}}(g_{\ell}^{(\ell,i)}-g_{\ell-1}^{(\ell,i)})$,
where $g_{-1}\equiv 0$, so that

$\displaystyle\mathbb{E}\left[Y_{\ell}\right]:=\begin{cases}\mathbb{E}\left[g_{0}\right],&\ell=0\\\
\mathbb{E}\left[g_{\ell}-g_{\ell-1}\right],&\ell>0\\\ \end{cases}.$ (25)

Denote by $Y:=\sum_{\ell=0}^{L}Y_{\ell}$ the multilevel estimator of
$\mathbb{E}\left[g\right]$ based on $L+1$ levels and $m_{\ell}$ independent
samples on level $\ell$, where $\ell=0,\ldots,L$.

Furthermore, we denote $V_{0}=\mathbb{V}\left[g_{0}\right]$ and for $\ell\geq
1$, let $V_{\ell}$ be the variance of $g_{\ell}-g_{\ell-1}$:
$V_{\ell}:=\mathbb{V}\left[g_{\ell}-g_{\ell-1}\right]$.

The standard theory states the following facts for the mean and the variance:

$\mathbb{E}\left[Y\right]=\mathbb{E}\left[g_{L}\right],\qquad\mathbb{V}\left[Y\right]={\textstyle\sum_{\ell=0}^{L}m_{\ell}^{-1}V_{\ell}}.$
(26)

The cost of the multilevel estimator $Y$ is

$S:={\textstyle\sum_{\ell=0}^{L}m_{\ell}s_{\ell}}.$ (27)

The mean squared error (MSE) is used to measure the quality of the multilevel
estimator:

$\text{MSE}:=\mathbb{E}\left[(Y-\mathbb{E}\left[g\right])^{2}\right]=\mathbb{V}\left[Y\right]+\left(\mathbb{E}\left[Y\right]-\mathbb{E}\left[g\right]\right)^{2},$
(28)

where $Y$ is what we computed via MLMC, and $\mathbb{E}\left[g\right]$ what
actually should be computed. To achieve

$\text{MSE}\leq\varepsilon^{2}$

for some prescribed tolerance $\varepsilon$, we ensure that

$\left(\mathbb{E}\left[Y\right]-\mathbb{E}\left[g\right]\right)^{2}=(\mathbb{E}\left[g_{L}-g\right])^{2}\leq\tfrac{1}{2}\varepsilon^{2}$
(29)

and

$\mathbb{V}\left[Y\right]\leq\tfrac{1}{2}\varepsilon^{2}.$ (30)

The bias error $\mathbb{E}\left[g_{L}-g\right]$ corresponds to the
discretization error discussed in Sec. 3.1. Later we will see that
$\mathbb{E}\left[Y\right]-\mathbb{E}\left[g\right]=O(2^{-\alpha L})$ with
$\alpha\approx 1$. The bias error can be made smaller than $\varepsilon^{2}/2$
by choosing a sufficiently large $L$. Then, for this $L$, we can compute
optimal $m_{0},\dots,m_{L}$ by formula in (34) to provide (30).

In the following, we repeat the well-known [18] results on the computation of
the sequence $m_{0},\dots,m_{L}$. For a fixed variance
$\mathbb{V}\left[Y\right]=:\varepsilon^{2}/2$, the cost $S$ is minimized by
choosing as $m_{\ell}$ the solution of the optimization problem

$\mbox{min}_{m_{0},\ldots,m_{L}}F(m_{0},\ldots,m_{L}),$ (31)

where
$F(m_{0},\ldots,m_{L}):=\sum_{\ell=0}^{L}\left(m_{\ell}s_{\ell}+\mu^{2}\frac{V_{\ell}}{m_{\ell}}\right)$,
$\mu^{2}$ is a Lagrange multiplier. Thus, the derivatives of $F$ w.r.t.
$m_{\ell}$ are equal to zero:

$\frac{\partial F(m_{0},\ldots,m_{L})}{\partial
m_{\ell}}:=s_{\ell}-\mu^{2}\frac{V_{\ell}}{m_{\ell}^{2}}=0.$ (32)

Solving the system (32), we obtain

$m_{\ell}^{2}=\mu^{2}\frac{V_{\ell}}{s_{\ell}},\quad\text{i.e.}\quad
m_{\ell}=\mu\sqrt{\frac{V_{\ell}}{s_{\ell}}}.$ (33)

Taking into account that the variation $\mathbb{V}\left[Y\right]$ is fixed and
substituting (33) into (26), i.e.
$\sum_{\ell=0}^{L}V_{\ell}/m_{\ell}=\varepsilon^{2}/2$, we obtain an equation
for $\mu$:

$\sum_{\ell=0}^{L}\frac{V_{\ell}}{\mu\sqrt{\frac{V_{\ell}}{s_{\ell}}}}=\tfrac{1}{2}\varepsilon^{2}.$

From this equation, we get
$\mu=2\varepsilon^{-2}\sum_{\ell=0}^{L}\sqrt{V_{\ell}s_{\ell}}$, and therefore

$m_{\ell}=2\varepsilon^{-2}\cdot\sqrt{\frac{V_{\ell}}{s_{\ell}}}\cdot\sum_{i=0}^{L}\sqrt{V_{i}s_{i}}.$
(34)

For this set of $m_{\ell}$, the total computational cost of $Y$ is

$S=2\varepsilon^{-2}\left(\sum_{\ell=0}^{L}\sqrt{V_{\ell}s_{\ell}}\right)^{2}.$
(35)

For further analysis of this sum, see [18], p.4.

The cost increases exponentially with $\ell$ while the weak error
$\mathbb{E}\left[g_{L}-g\right]$ and multilevel correction variance $V_{\ell}$
decrease exponentially leads to the following theorem (cf. Theorem 1, p. 6 in
[18]):

######  Theorem 2.

Consider a fixed $t=t^{*}$. Suppose positive constants $\alpha,\beta,\gamma>0$
exist such that $\alpha\geq\frac{1}{2}\mbox{min}(\beta,\gamma\hat{d})$, and

$\displaystyle|\mathbb{E}\left[g_{\ell}-g\right]|$ $\displaystyle\leq
c_{1}2^{-\alpha\ell}$ (36a) $\displaystyle V_{\ell}$ $\displaystyle\leq
C_{2}2^{-\beta\ell}$ (36b) $\displaystyle s_{\ell}$ $\displaystyle\leq
c_{3}2^{\hat{d}\gamma\ell}.$ (36c)

Then, for any accuracy $\varepsilon<e^{-1}$, a constant $c_{4}>0$ and a
sequence of realizations $\\{m_{\ell}\\}_{\ell=0}^{L}$ exist, such that
$\text{MSE}<\varepsilon^{2}$, where MSE is defined in (28), and the
computational cost is

$\displaystyle S=\begin{cases}c_{4}{\varepsilon^{-2}},&\beta>\hat{d}\gamma\\\
c_{4}{\varepsilon^{-2}\left(\log(\varepsilon)\right)^{2}},&\beta=\hat{d}\gamma\\\
c_{4}{\varepsilon^{-\left(2+\frac{\hat{d}\gamma-\beta}{\alpha}\right)}},&\beta<\hat{d}\gamma.\\\
\end{cases}$ (37)

This theorem (see also [27, 26, 8, 9, 17]) indicates that, even in the worst-
case scenario, the MLMC algorithm has a lower computational cost than that of
the traditional (single-level) MC method, which scales as
$\mathcal{O}(\varepsilon^{-2-\hat{d}\gamma/\alpha})$.

######  Remark 1.

In Theorem 2, the factors $C_{1}$, $C_{2}$, $C_{3}$, $C_{4}$ as well as the
exponents $\alpha$, $\beta$ and $\gamma$ depend on the time point $t$. This
makes $L$ and $m_{\ell}$ time-dependent, too.

######  Remark 2.

One possible choice for a scaling factor in (29) and (30) is
$E_{0}:=|\mathbb{E}\left[g_{0}(t^{*},\mathbf{x}^{*})\right]|$, where
$g_{0}(t^{*},\mathbf{x}^{*})$ is the solution computed on level $\ell=0$ at
point $(t^{*},\mathbf{x}^{*})$.

######  Remark 3.

We consider the error relatively to $g$. For this, in (29) and (30), we
replace $\varepsilon$ by $\varepsilon\cdot E_{0}$. Equivalently, we can divide
(rescale) $\mathbb{E}\left[g_{L}-g\right]$ by $E_{0}$ and $V_{i}$ by
$E_{0}^{2}$. Therefore, we get:

${|\mathbb{E}\left[g_{L}-g\right]|}\leq C_{1}2^{-\alpha L}.$ (38)

Now, to satisfy (29), we want

${|\mathbb{E}\left[g_{L}-g\right]|}\leq\frac{1}{\sqrt{2}}\varepsilon E_{0}.$
(39)

From this inequality, we can estimate $L$:

$C_{1}2^{-\alpha L}=\frac{\varepsilon E_{0}}{\sqrt{2}}$ (40)
$L=-\tfrac{1}{\alpha}\log_{2}\frac{\varepsilon E_{0}}{\sqrt{2}C_{1}}.$ (41)

Equations (34) and (35) attain the form

$m_{\ell}=\frac{2\varepsilon^{-2}}{{E_{0}}^{2}}\cdot\sqrt{\frac{V_{\ell}}{s_{\ell}}}\cdot\sum_{i=0}^{L}\sqrt{V_{i}s_{i}},\qquad
S=\frac{2\varepsilon^{-2}}{E_{0}^{2}}\left(\sum_{\ell=0}^{L}\sqrt{V_{\ell}s_{\ell}}\right)^{2}.$
(42)

Using preliminary numerical tests (see Fig. 5) and following preprocessing
steps as in Algorithm 1, we can estimate the convergence rates $\alpha$ for
the mean (the so-called weak convergence) and $\beta$ for the variance (the
so-called strong convergence), as well as the constants $C_{1}$ and $C_{2}$.
In addition, $\alpha$ is strongly connected to the order of the discretization
error (see Section 3.1), which equals $1$. Note that precise estimates of
parameters $\alpha$ and $\beta$ are crucial to distribute the computational
effort optimally.

Algorithm 1 Preprocessing for MLMC algorithm

Input: $L=3$, $m_{\ell}=10$, $\ell=0,\ldots,L$

Compute $g_{\ell}(\omega_{i})$, $g_{\ell}(\omega_{i})-g_{\ell-1}(\omega_{i})$
for $i=1\ldots m_{\ell}$ and $\ell=1\ldots L$

Estimate $\mathbb{E}\left[g_{\ell}(\omega_{i})-g_{\ell-1}(\omega_{i})\right]$,
$\mathbb{V}\left[g_{\ell}(\omega_{i})-g_{\ell-1}(\omega_{i})\right]$

Estimate convergence rates $\alpha$ and $\beta$, complexity rate $\gamma$,
constants $C_{1},C_{2},C_{3}$

Algorithm 2 MLMC algorithm

Input: MSE Error $\varepsilon^{2}$

Estimate $L$

Compute $m_{\ell}$, $\ell=0,\ldots,L$

Compute $\mathbb{E}\left[g_{\ell}(\omega_{i})-g_{\ell-1}(\omega_{i})\right]$,
$i=1\ldots m_{\ell}$

Build the telescopic sum to estimate $\mathbb{E}\left[g\right]$

##  4 Numerical Experiments

We performed numerical experiments with four uncertain parameters: thickness,
recharge, permeability and porosity. Permeability and porosity are mutually
dependent. To check that uncertainties in all four parameters contribute
significantly to the solution, we performed auxiliary experiments with only
1-2 random parameters. First, only one of these four parameters was uncertain,
then thickness and recharge were uncertain. In this way we checked that there
were no parameters whose contribution could be neglected.

We consider the following QoIs: 1) the solution $c_{m}$ at a point
$(t,\mathbf{x})$, i.e. $g=c_{m}(t,\mathbf{x})$, and 2) an integral over a
small sub-domain

$I_{i}(t,\omega):=\int_{\mathbf{x}\in\Delta_{i}}c_{m}(t,\mathbf{x},\omega)\rho(c_{m}(t,\mathbf{x},\omega))d\mathbf{x},\quad\text{where}$
(43) $\Delta_{i}:=[x_{i}-0.1,x_{i}+0.1]\times[y_{i}-0.1,y_{i}+0.1],\quad
i=1,\ldots,6.$ (44)

The list of all points $\mathbf{x}_{i}=(x_{i},y_{i})$ is defined in (2.2). The
value of $I_{i}$ is the mass of the salt in a subdomain $\Delta_{i}$. The size
of each $\Delta_{i}$ is small ($0.2^{2}=0.04$), compared to $\mathcal{D}$.

The numerical scheme provides only the first order of the accuracy to compute
$c_{m}(t,\mathbf{x}_{i},\omega)$ and $I_{i}(t,\omega)$, i.e., the convergence
rate $\alpha$ (weak convergence in (36a)) should be $\approx 1$. The numerical
results below show that, indeed, $\alpha\approx 1$.

Table 2 contains computing times needed to compute the solution $c_{m}$ at
each level $\ell$. The fifth column contains the average computing time, and
the sixth and seventh columns contain the shortest and longest computing
times. The computing time for each simulation varies depending on the number
of iterations, which depends on the porosity and permeability. We observed
that, after $\approx 6016$ sec., the solution is almost unchanging; thus, we
perform the experiment only for $t\in[0,T]$, where $T=6016$ sec. For example,
if the number of time steps is $r_{\ell}=188$ (Level 0 in Table 2), then the
time step $\tau=\frac{T}{r_{\ell}}=\frac{6016}{188}=32$ sec.

Level $\ell$ | $n_{\ell}$, ($\frac{n_{\ell}}{n_{\ell-1}}$) | $r_{\ell}$, ($\frac{r_{\ell}}{r_{\ell-1}}$) | $\tau_{\ell}=6016/r_{\ell}$ | Computing times ($s_{\ell}$), ($\frac{s_{\ell}}{s_{\ell-1}}$)
---|---|---|---|---
average | min. | max.
0 | 608 | 188 | 32 | 3 | 2.4 | 3.4
1 | 2368 (3.9) | 376 (2) | 16 | 22 (7.3) | 15.5 | 27.8
2 | 9344 (3.9) | 752 (2) | 8 | 189 (8.6) | 115 | 237
3 | 37120 (4) | 1504 (2) | 4 | 1831 (10) | 882 | 2363
4 | 147968 (4) | 3008 (2) | 2 | 18580 (10) | 7865 | 25418
Table 2: Number of degrees of freedom $n_{\ell}$, number of time steps
$r_{\ell}$, step size in time $\tau_{\ell}$, average, minimal, and maximal
computing times on each level $\ell$.

In Figure 2 (top left) we visualise the coefficient of variation $CV(g)$,
which is the ratio of the standard deviation to the mean. The variance
$V_{\ell}:=\mathbb{V}\left[g_{\ell}\right]$ is shown on the top right, the
mean $\mathbb{E}\left[g_{\ell}-g_{\ell-1}\right]$ on the bottom left and the
variance $\mathbb{V}\left[g_{\ell}-g_{\ell-1}\right]$ on the bottom right. The
QoI is $g:=c_{m}(t,\mathbf{x}_{1})$ and the time $t=0,\ldots,47\tau$ is along
the $x$ axis.

Figure 2: (top left) The coefficient of variance $CV_{\ell}:=CV(g_{\ell})$,
(top right) the variance $\mathbb{V}\left[g_{\ell}\right]$, (bottom left) the
mean $\mathbb{E}\left[g_{\ell}-g_{\ell-1}\right]$, (bottom right) the variance
$\mathbb{V}\left[g_{\ell}-g_{\ell-1}\right]$. The QoI is
$g=c_{m}(t,\mathbf{x}_{1})$. The small oscillations in the two lower pictures
are due to the dependence of the recharge $\hat{q}_{x}$ on the time, cf. (19).

###  4.1 The mean value and variance

Fig. 3 shows the mean value $\mathbb{E}\left[c_{m}(t,\mathbf{x})\right]$ at
time points $\\{7\tau,19\tau,40\tau,94\tau\\}$. In all cases
$\mathbb{E}\left[c_{m}(t,\mathbf{x})\right]\in[0,1]$, the dark red color
corresponds to the value 1 and dark blue to the value 0.

Similarly, Fig. 4 presents the variance of the salt mass fraction in the
computational domain for $t=\\{7\tau,19\tau,40\tau,94\tau\\}$. The maximal
values of $\mathbb{V}\left[c_{m}\right]$ (denoted by the dark red colour) are
$1.9\cdot 10^{-3}$, $3.4\cdot 10^{-3}$, $2.9\cdot 10^{-3}$, $2.4\cdot 10^{-3}$
respectively.

The general flow and transport patters in this setting were discussed in Sect.
2.2. The aperture $\epsilon$ of the fracture has the most essential influence
on the $c_{m}$. If $\epsilon$ is small (i.e. for $\xi_{1}=-1$, cf. (18)), the
fracture does not increase the permeability of the medium significantly so the
flow and transport are very similar to the Henry problem. For thick fractures
$(\xi_{1}=1)$, the separation of the upper and the lower part of the domain
(see Sect. 2.2) is very strong. This is clearly visible in the mean value of
$c_{m}(t,\mathbf{x})$ in Fig. 3. Note that $\epsilon$ changes the flow under
the left end of the fracture where fluid enters it. This area is characterized
by the high value of the variance, see Fig. 4. The influence of the other
uncertainties is weaker but nevertheless important, too.

Figure 3: The mean value $\mathbb{E}\left[c_{m}(t,\mathbf{x})\right]$ at
$t=\\{7\tau,19\tau,40\tau,94\tau\\}$. In all cases,
$\mathbb{E}\left[c_{m}\right](t,\mathbf{x})\in[0,1]$.

Figure 4: The variance $\mathbb{V}\left[c_{m}(t,\mathbf{x})\right]$ at
$t=\\{7\tau,19\tau,40\tau,94\tau\\}$. Maximal values (dark red colour) of
$\mathbb{V}\left[c_{m}\right]$ are $1.9\cdot 10^{-3}$, $3.4\cdot 10^{-3}$,
$2.9\cdot 10^{-3}$, $2.4\cdot 10^{-3}$ respectively. The dark blue colour
corresponds to a zero value.

###  4.2 Estimation of convergence rates $\alpha$ and $\beta$ for various QoIs

Figure 5 (left) shows the decay of
$\mathbb{E}\left[g_{\ell}-g_{\ell-1}\right]$ as a function of
$g_{\ell}-g_{\ell-1}$. This dependence is fitted by a
$2^{-\alpha\ell+\zeta_{1}}$ function. The parameter $\alpha$ indicates the
weak convergence rate. Figure 5 (right) shows the decay of
$\mathbb{V}\left[g_{\ell}-g_{\ell-1}\right]$ with respect to
$g_{\ell}-g_{\ell-1}$. This dependence is fitted by a
$2^{-\beta\ell+\zeta_{2}}$ curve, where the parameter $\beta$ indicates the
strong convergence rate. The constant are $C_{1}=0.47$, $C_{2}=3.7\cdot
10^{-3}$. The QoI is $g:=c_{m}(t_{15},\mathbf{x}_{1})$. The calculated rates
$\alpha=1.07$ and $\beta=1.97$ are very close to the theoretical values, which
are $1$ and $2$ respectively.

Figure 5: The weak and the strong convergences, QoI is
$g:=c_{m}(t_{15},\mathbf{x}_{1})$, $\alpha=1.07$, $\zeta_{1}=-1.1$,
$\beta=1.97$, $\zeta_{2}=-8$.

In Fig. 6(left) the decays of $\mathbb{E}\left[g_{\ell}-g_{\ell-1}\right]$ and
$\mathbb{E}\left[g_{\ell}\right]$ are plotted against the level $\ell$. The
QoI is $c_{m}(t_{18},\mathbf{x}_{1})$, $\mathbf{x}_{1}=(1.1,-0.8)$. Similarly,
Figure 6 (right) shows that $\mathbb{V}\left[g_{\ell}-g_{\ell-1}\right]$
decays much faster than $\mathbb{V}\left[g_{\ell}\right]$.

Figure 6: Decay comparison of (left)
$\mathbb{E}\left[g_{\ell}-g_{\ell-1}\right]$ and
$\mathbb{E}\left[g_{\ell}\right]$ vs. $\ell$; (right)
$\mathbb{V}\left[g_{\ell}-g_{\ell-1}\right]$ and
$\mathbb{V}\left[g_{\ell}\right]$. QoI $g=c_{m}(t_{18},\mathbf{x}_{1})$,
$\mathbf{x}_{1}=(1.1,-0.8)$.

In Table 3 we provide estimations for the weak and strong convergence rates
$\alpha$ and $\beta$, as well as constants $C_{1}$ and $C_{2}$ for different
QoIs. Each QoI is the solution at a point $\mathbf{x}_{i}$ or an integral
$I_{i}$ as in (43). We see that the value $\alpha$ is very close to one. This
is fully corresponding to the theory. The estimated value of $\beta$ is also
very close to the theoretical value, which is equal to $2\alpha=2$. The values
of $\alpha$ and $C_{1}$ are used for estimation of the number of needed levels
$L$, see [31].

QoI | $\alpha$ | $c_{1}$ | $\beta$ | $C_{2}$
---|---|---|---|---
$c_{m}(\mathbf{x}_{1},15\tau)$ | -1.07 | 0.47 | -1.97 | $4\cdot 10^{-3}$
$c_{m}(\mathbf{x}_{2},15\tau)$ | -0.96 | 0.25 | -1.6 | $3.1\cdot 10^{-4}$
$c_{m}(\mathbf{x}_{3},15\tau)$ | -0.9 | 0.12 | -1.9 | $4.5\cdot 10^{-5}$
$c_{m}(\mathbf{x}_{4},15\tau)$ | -1.07 | 0.47 | -1.94 | $3.6\cdot 10^{-3}$
$c_{m}(\mathbf{x}_{5},15\tau)$ | -0.95 | 0.25 | -1.5 | $1.9\cdot 10^{-4}$
$c_{m}(\mathbf{x}_{6},15\tau)$ | -0.89 | 0.1 | -1.82 | $9.2\cdot 10^{-5}$
$I_{2}(c_{m})$ | -1.5 | 8.3 | -2.6 | 0.4
Table 3: Estimated the weak and strong convergence rates $\alpha$ and $\beta$,
and constants $C_{1}$ and $C_{2}$ for different QoIs. QoI are solution in a
point $\mathbf{x}_{i}$ and an integral $I_{i}$ as in (43).

In Table 4 we show the numbers $m_{\ell}$ needed to get the MSE error below
$\varepsilon^{2}$, and $\varepsilon$ should be understood as a relative error,
as it is defined in (39). All numbers are obtained for the QoI
$g=c_{m}(t_{15},\mathbf{x}_{1})$. The numbers $m_{ell}$ can be very different
for other QoIs.

$\varepsilon$ | $m_{0}$ | $m_{1}$ | $m_{2}$ | $m_{3}$ | $m_{4}$
---|---|---|---|---|---
0.2 | 28 | 0 | 0 | 0 | 0
0.1 | 10 | 1 | 0 | 0 | 0
0.05 | 57 | 4 | 1 | 0 | 0
0.025 | 342 | 22 | 4 | 1 | 0
0.01 | 3278 | 205 | 35 | 6 | 1
Table 4: The number of samples $m_{\ell}$ vs. $\varepsilon$, QoI
$g=c_{m}(t_{15},\mathbf{x}_{1})$.

###  4.3 Comparison of MC and MLMC methods

In Figure 7 we compare the estimated costs of the MC (denoted by a blue dashed
line) and MLMC (solid red line) methods, see the two bottom graphics. We also
plot the $\varepsilon^{2}$ graphic (yellow dashed line) and the theoretical MC
and MLMC costs, purple dashed line and green dashed line respectively. The
theoretical MC cost is
$\mathcal{O}(\varepsilon^{-2-\frac{\hat{d}\gamma}{\alpha}})$, with
$\hat{d}=2+1=3$, $\gamma=1$, $\alpha=1.08$. One can see that from a certain
epsilon $\varepsilon\approx 8\cdot 10^{-3}$ the MLMC outperforms the MC
method. It can also be seen that the theoretical MC and MLMC lines (two top
graphics) have the same slope as the calculated MC and MLMC lines (two bottom
graphics). The theoretical MC and MLMC graphics are scaled by a constant for
better visibility.

Figure 7: Complexity comparison of ML and MLMC methods against $\varepsilon$
(horizontal axis) in log-log scale.

##  5 Discussion and Conclusion

In this work we have considered a setting that mimics the Henry problem [42,
43], modeling seawater intrusion into a 2D coastal aquifer. The pure water
recharge from the “land side” resists the salinisation of the aquifer due to
the influx of saline water through the “sea side”, thereby achieving some
equilibrium in the salt concentration. In our setting, following [20], we
consider a fracture on the sea side that significantly increases the
permeability of the porous medium.

The flow and transport essentially depend on the geological parameters of the
porous medium, including the fracture. We investigated the effects of various
uncertainties on saltwater intrusion. We assumed uncertainties in the fracture
width, the porosity of the bulk medium, its permeability and the pure water
recharge from the land side. The porosity and permeability were modeled by
random fields, the recharge by a random but periodic intensity and the
thickness by a random variable. We calculated the mean and variance of the
salt mass fraction, which is also uncertain.

The main question we investigated in this work was how well the MLMC method
can be used to compute statistics of different QoIs. We found that the answer
depends on the choice of the QoI. First, not every QoI requires a hierarchy of
meshes and MLMC. Second, MLMC requires stable convergence rates for
$\mathbb{E}\left[g_{\ell}-g_{\ell-1}\right]$ and
$\mathbb{V}\left[g_{\ell}-g_{\ell-1}\right]$. These rates should be
independent of $\ell$. If these convergence rates vary for different $\ell$,
then it will be hard to estimate $L$ and $m_{\ell}$, and MLMC will either not
work or be suboptimal. We were not able to get stable convergence rates for
all levels $\ell=1,\ldots,5$ when the QoI was an integral as in (43). We found
that for $\ell=1,\ldots 4$ and $\ell=5$ the rate $\alpha$ was different.
Further investigation is needed to find the reason for this. Another
difficulty is the dependence on time, i.e. the number of levels $L$ and the
number of sums $m_{\ell}$ depend on $t$. At the beginning the variability is
small, then it increases, and after the process of mixing salt and fresh water
has stopped, the variance decreases again.

The number of random samples required at each level was estimated by
calculating the decay of the variances and the computational cost for each
level. These estimates depend on the minimisation function in the MLMC
algorithm.

To achieve the efficiency of the MLMC approach presented in this work, it is
essential that the complexity of the numerical solution of each random
realisation is proportional to the number of grid vertices on the grid levels.
This holds only under special choice the numerical solvers, and we applied the
geometric multigrid method implemented in the ug4 software toolkit, cf. [37,
49]. Our numerical tests confirm the expected efficiency. The random
realizations do not depend on each other and have been computed concurrently
on different nodes of the supercomputer Shaheen II at the King Abdullah
University of Science and Technology.

In the numerical experiments, we computed the expected value and variance of
the mass fraction in the whole domain and the solution at a few preselected
points $(t,\mathbf{x})$. It turns out that the grid hierarchy of only 2 to 3
levels is sufficient for some QoIs. Note that a different type of porosity in
(20) may lead to a different conclusion.

Our numerical results confirm that the MLMC is principally more numerically
efficient than the MC method. Thus, sampling at different mesh levels helps to
reduce the overall computational cost.
Constraints. 1\. The time dependence is challenging. The optimal number of
samples depends on the time moment $t$, and the spatial location $\mathbf{x}$,
and may be small for some points and large for others. 2. Each new QoI
requires new estimates of all parameters $\alpha$, $\beta$, $\gamma$, and the
resulting MLMC performance graphics can be different.
Future work. A more realistic modeling of the fracture, porosity and
permeability with more random variables can be considered.

Acknowledgments For computing time, this research used Shaheen II, which is
managed by the Supercomputing Core Laboratory at the King Abdullah University
of Science and Technology (KAUST) in Thuwal, Saudi Arabia. We thank the KAUST
HPC support team for their assistance with Shaheen II. This work was supported
by the Alexander von Humboldt foundation.

## References

  *   * Angot et al [2009] Angot P, Boyer F, Hubert F (2009) Asymptotic and numerical modelling of flow in fractured porous media. ESAIM: M2AN Mathematical Modelling and Numerical Analysis 43:239–275
  * Barrett et al [1994] Barrett R, Berry M, Chan TF, et al (1994) Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Society for Industrial and Applied Mathematics, 10.1137/1.9781611971538
  * Bastian et al [2000a] Bastian P, Birken K, Johannsen K, et al (2000a) Parallel solution of partial differential equations with adaptive multigrid methods on unstructured grids. In: Jäger W, Krause E (eds) High performance computing in science and engineering. Springer, Berlin, pp 506–519
  * Bastian et al [2000b] Bastian P, Chen Z, Ewing RE, et al (2000b) Numerical simulation of multiphase flow in fractured porous media. In: Chen Z, et al. (eds) Numerical treatment of multiphase flows in porous media. Proceedings of the international workshop, Beijing, China, August 2–6, 1999., Lect. Notes Phys., vol 552. Springer, Berlin, pp 50–68
  * Blatman and Sudret [2010] Blatman G, Sudret B (2010) An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis. Probabilistic Engineering Mechanics 25(2):183 – 197. http://dx.doi.org/10.1016/j.probengmech.2009.10.003
  * Bode et al [2018] Bode F, Ferré T, Zigelli N, et al (2018) Reconnecting stochastic methods with hydrogeological applications: A utilitarian uncertainty analysis and risk assessment approach for the design of optimal monitoring networks. Water Resources Research 54(3):2270–2287. https://doi.org/10.1002/2017WR020919
  * Carrera [1993] Carrera J (1993) An overview of uncertainties in modelling groundwater solute transport. Journal of Contaminant Hydrology 13(1):23 – 48. https://doi.org/10.1016/0169-7722(93)90049-X, chemistry and Migration of Actinides and Fission Products
  * Charrier et al [2013] Charrier J, Scheichl R, Teckentrup AL (2013) Finite element error analysis of elliptic pdes with random coefficients and its application to multilevel monte carlo methods. SIAM Journal on Numerical Analysis 51(1):322–352. 10.1137/110853054
  * Cliffe et al [2011] Cliffe K, Giles M, Scheichl R, et al (2011) Multilevel monte carlo methods and applications to elliptic pdes with random coefficients. Computing and Visualization in Science 14(3):3–15. 10.1007/s00791-011-0160-x
  * Collier et al [2015] Collier N, Haji-Ali AL, Nobile F, et al (2015) A continuation multilevel monte carlo algorithm. BIT Numerical Mathematics 55(2):399–432
  * Costa [2006] Costa A (2006) Permeability-porosity relationship: A reexamination of the kozeny-carman equation based on a fractal pore-space geometry assumption. Geophysical Research Letters 33(2). 10.1029/2005GL025134
  * Cremer et al [2015] Cremer C, , Graf T (2015) Generation of dense plume fingers in saturated–unsaturated homogeneous porous media. Journal of Contaminant Hydrology 173:69 – 82. https://doi.org/10.1016/j.jconhyd.2014.11.008
  * D’Angelo and Scotti [2012] D’Angelo C, Scotti A (2012) A mixed finite element method for darcy flow in fractured porous media with non-matching grids. ESAIM: Mathematical Modelling and Numerical Analysis 46(02):465–489. 10.1051/m2an/2011148
  * Dhal and Swain [2022] Dhal L, Swain S (2022) Understanding and modeling the process of seawater intrusion: a review, pp 269–290. 10.1016/B978-0-12-823830-1.00009-2
  * Frolkovič and De Schepper [2001] Frolkovič P, De Schepper H (2001) Numerical modelling of convection dominated transport coupled with density driven flow in porous media. Advances in Water Resources 24(1):63–72. 10.1016/S0309-1708(00)00025-7
  * Fumagalli and Scotti [2013] Fumagalli A, Scotti A (2013) A numerical method for two-phase flow in fractured porous media with non-matching grids. Advances in Water Resources 62:454–464. 10.1016/j.advwatres.2013.04.001, computational Methods in Geologic CO2 Sequestration
  * Giles [2008] Giles MB (2008) Multilevel Monte Carlo path simulation. Operations Research 56(3):607–617
  * Giles [2015] Giles MB (2015) Multilevel Monte Carlo methods. Acta Numerica 24:259–328
  * Graf and Therrien [2007] Graf T, Therrien R (2007) Variable-density Groundwater Flow and Solute Transport in Irregular 2D Fracture Networks. Adv Water Resour 30:455–468
  * Grillo et al [2010] Grillo A, Logashenko D, Stichel S, et al (2010) Simulation of density-driven flow in fractured porous media. Advances in Water Resources 33(12):1494–1507. 10.1016/j.advwatres.2010.08.004
  * Grillo et al [2012] Grillo A, Lampe M, Logashenko D, et al (2012) Simulation of salinity- and thermohaline-driven flow in fractured porous media. Journal of Porous Media 15(5):439–458. 10.1615/JPorMedia.v15.i5.40
  * Hackbusch [1985] Hackbusch W (1985) Multi-Grid Methods and Applications. Springer, Berlin
  * Hackbusch [1994] Hackbusch W (1994) Iterative Solution of Large Sparse Systems of Equations. Springer, New-York
  * Haji-Ali et al [2016] Haji-Ali AL, Nobile F, von Schwerin E, et al (2016) Optimization of mesh hierarchies in multilevel Monte Carlo samplers. Stoch Partial Differ Equ Anal Comput 4(1):76–112. 10.1007/s40072-015-0049-7
  * Henry [1964] Henry HR (1964) Effects of dispersion on salt encroachment in coastal aquifers, in ’seawater in coastal aquifers’. US Geological Survey, Water Supply Paper 1613:C70–C80
  * Hoel et al [2012] Hoel H, Von Schwerin E, Szepessy A, et al (2012) Adaptive multilevel Monte Carlo simulation. In: Numerical Analysis of Multiscale Computations. Springer, p 217–234
  * Hoel et al [2014] Hoel H, von Schwerin E, Szepessy A, et al (2014) Monte Carlo Methods and Applications 20(1):1–41. doi:10.1515/mcma-2013-0014
  * Hoteit and Firoozabadi [2008] Hoteit H, Firoozabadi A (2008) An efficient numerical model for incompressible two-phase flow in fractured media. Advances in Water Resources 31(6):891–905. 10.1016/j.advwatres.2008.02.004
  * Khoromskij et al [2009] Khoromskij BN, Litvinenko A, Matthies H (2009) Application of hierarchical matrices for computing the Karhunen-Loève expansion. Computing 84(1-2):49–67. 10.1007/s00607-008-0018-3
  * Litvinenko et al [2019] Litvinenko A, Yucel AC, Bagci H, et al (2019) Computation of electromagnetic fields scattered from objects with uncertain shapes using multilevel monte carlo method. IEEE Journal on Multiscale and Multiphysics Computational Techniques 4:37–50. 10.1109/JMMCT.2019.2897490
  * Logashenko et al [2024] Logashenko D, Litvinenko A, Tempone R, et al (2024) Uncertainty quantification in the henry problem using the multilevel monte carlo method. Journal of Computational Physics 503:112854. https://doi.org/10.1016/j.jcp.2024.112854
  * Martin et al [2005] Martin V, Jaffré J, Roberts JE (2005) Modeling fractures and barriers as interfaces for flow in porous media modeling fractures and barriers as interfaces for flow in porous media. SIAM J Sci Comput 26(5):1667—1691. 10.1137/S1064827503429363
  * Martinez-Landa and Carrera [2006] Martinez-Landa L, Carrera J (2006) A methodology to interpret cross-hole tests in a granite block. Journal of Hydrogeology 325(1-4):222–240
  * Panda and Lake [1994] Panda M, Lake W (1994) Estimation of single-phase permeability from parameters of particle-size distribution. AAPG Bull 78:1028–1039
  * Pape et al [1999] Pape H, Clauser C, Iffland J (1999) Permeability prediction based on fractal pore-space geometry. Geophysics 64(5):1447–1460. 10.1190/1.1444649
  * Reiter et al [2012] Reiter S, Logashenko D, Grillo A, et al (2012) Preparation of grids for simulations of groundwater flow in fractured porous media. Computing and Visualization in Science 15:209–225. 10.1007/s00791-013-0210-7
  * Reiter et al [2013] Reiter S, Vogel A, Heppner I, et al (2013) A massively parallel geometric multigrid solver on hierarchically distributed grids. Computing and Visualization in Science 16(4):151–164. 10.1007/s00791-014-0231-x
  * Reiter et al [2014] Reiter S, Logashenko D, Stichel S, et al (2014) Models and simulations of variable-density flow in fractured porous media. International Journal of Computational Science and Engineering 9(5-6):416–432. 10.1504/IJCSE.2014.064527
  * Riva et al [2015] Riva M, Guadagnini A, Dell’Oca A (2015) Probabilistic assessment of seawater intrusion under multiple sources of uncertainty. Advances in Water Resources 75(C):93–104. 10.1016/j.advwatres.2014.11.002
  * Rubin [2003] Rubin Y (2003) Applied Stochastic Hydrogeology. Oxford University Press, 10.1093/oso/9780195138047.001.0001
  * Shikaze et al [1998] Shikaze SG, Sudicky EA, Schwartz FW (1998) Density-dependent solute transport in discretely-fractured geologic media: is prediction possible? Journal of Contaminant Hydrogeology 34:273–291
  * Simpson and Clement [2003] Simpson MJ, Clement T (2003) Theoretical Analysis of the worthiness of Henry and Elder problems as benchmarks of density-dependent groundwater flow models. Adv Water Resour 26:17–31
  * Simpson and Clement [2004] Simpson MJ, Clement TP (2004) Improving the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models. Water Resources Research 40(1):W01504. 10.1029/2003WR002199
  * Sorek et al [2001] Sorek S, Borisov V, Yakirevich A (2001) A two-dimensional Areal Model for Density Dependent Flow Regime. Transport in Porous Media 43:87–105
  * Sudret [2008] Sudret B (2008) Global sensitivity analysis using polynomial chaos expansions. Reliability engineering & system safety 93(7):964–979
  * Tartakovsky [2013] Tartakovsky D (2013) Assessment and management of risk in subsurface hydrology: A review and perspective. Advances in Water Resources 51:247 – 260. https://doi.org/10.1016/j.advwatres.2012.04.007, 35th Year Anniversary Issue
  * Teckentrup et al [2013] Teckentrup A, Scheichl R, Giles M, et al (2013) Further analysis of multilevel monte carlo methods for elliptic pdes with random coefficients. Numerische Mathematik 125(3):569–600
  * Vereecken et al [2016] Vereecken H, Schnepf A, Hopmans J, et al (2016) Modeling soil processes: Review, key challenges, and new perspectives. Vadose Zone Journal 15(5):vzj2015.09.0131. https://doi.org/10.2136/vzj2015.09.0131
  * Vogel et al [2013] Vogel A, Reiter S, Rupp M, et al (2013) Ug 4: A novel flexible software system for simulating pde based models on high performance computers. Computing and Visualization in Science 16(4):165–179. 10.1007/s00791-014-0232-9
  * Voss and Souza [1987] Voss C, Souza W (1987) Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone. Water Resources Research 23(10):1851–1866. 10.1029/WR023i010p01851


# Simple Natural Language Processing Tools
for Danish

Leon Derczynski
Department of Computer Science
ITU Copenhagen
Denmark DK2300
<EMAIL_ADDRESS>

###### Abstract

This technical note describes a set of baseline tools for automatic processing
of Danish text. The tools are machine-learning based, using natural language
processing models trained over previously annotated documents. They are
maintained at ITU Copenhagen and will always be freely available.

_K_ eywords natural language processing $\cdot$ Danish $\cdot$ open access

##  1 Introduction

Danish, a language with few resources for automatic processing, is spoken by
six million people, largely concentrated in the Scandinavian country of
Denmark. At the NLP (natural language processing) group at ITU Copenhagen, one
of our foci is Nordic and Danish languages. Thus, we aim to include local
languages in our NLP research whenever possible, and to provide tools for
others working on these languages. Despite a modest number of researcher in
the country, there has been a gap in usable Danish NLP tools that are
available to the wider community. As in every time that new technology is not
effectively mediated to those who can use it, this restriction stymies the
impact and development NLP in Denmark. To bridge this artificial lacuna, a
basic set of easy-to-run utilities that build on existing datasets and open
tools has been developed at ITU, which this note introduces.

The tools are available at https://nlp.itu.dk/resources/.

##  2 daner: Named Entity Recognition

Named entity recognition means finding words in text that refer to things by a
specific name. This could be e.g. using “Nanna" to refer to a specific
individual, instead of “them" or “she" or “the person"; or referring to
“Moscow" instead of “the nearest city". The important thing is that one
mentions a name specifically. Named entity recognition (NER) finds those
names. In the case of daner, names of Locations, Organisation and Persons are
tagged.

The daner tool wraps the CoreNLP [1] named entity recognition component [2],
using the DKIE data [3] developed as part of the GATE tool [4], which is
derived from the UD part of the Copenhagen Dependency Treebank [5], itself
including data from the Danish Dependency Treebank [6]. Further ad-hoc data is
added as required, from newswire and other sources. The data overall comprises
20K tokens from the CDT (largely 1990s newswire), 3K tokens from non-newswire
Danish, and 10K tokens from post-2015 newswire. The mode performs best over
Danish newswire text. The data is annotated for three classes, person (PER),
location (LOC), and organisation (ORG). Tags are in BIO format.

The tool recognises names of people, of organisations, and of locations. It
performs automatic tokenization, and outputs in slashed-tag format. For
example:

> En/O stor/O reform/O skal/O derfor/O blandt/O andet/O styrke/O tilliden/O
> til/O politikere/O og/O medier/O ,/O genopbygge/O tilliden/O til/O Skat/O
> og/O mindske/O de/O økonomiske/O forskelle/O i/O Danmark/B-LOC ./O

Here, “Danmark" is tagged as a location. The B indicates that “Danmark" is the
first token in the location name (the beginning ).

Additional data is added to an internal repository that is also used to train
daner; some of this data has copyright restrictions, and so cannot be
distributed, though the aggregated statistical representation in daner’s model
is free from that restriction, and so can be downloaded by anyone.

The NER model is configured to use contextual features, including n-grams up
to 5-grams, and to use type sequences and word shape features.

daner uses Brown clusters induced over large Danish-language corpus. These are
created using the Generalised Brown formulation [7], which offers more
flexibility in use than the classical version – specifically, one clustering
run can be used to generate clusterings of any size. For daner, 2500
clusters111A size that works well for NER in English [8] are induced from a
clustering with window size of 5000, over a set of 134M Danish tokens. These
134M tokens are taken from Danish Wikipedia and Danish content from the Common
Crawl corpus filtered through the FastText language ID [9], both using June
2019 data. The clusters are used inside the model, but also available
separately at http://itu.dk/~leod/dansk-brown.tar.bz2.

##  3 dapipe: Tokenization, Part-of-speech tagging, and Dependency Parsing

Another utility in the toolkit, dapipe, provides three functions: breaking
text files into words and sentences, finding what class each word is through
PoS tagging, and determining one part of the syntax that links a sentence’s
words together, through dependency parsing.

Defining where words and sentences start and stop is very important if we are
going to do NLP. These processes are both called tokenization. dapipe uses the
Universal dependencies method for tokenizing, breaking text into sentences and
words.

The next step is to give a “part of speech" to each word. This would be
something like a noun (substantiv), verb (verbum), or preposition
(præposition). This tells us the function each particular word has. dapipe
works generally well, though for social media text, I recommend
structbilty,222See https://github.com/bplank/bilstm-aux which is more
resilient [10, 11] to the unusual words and noise found in this colloquial
setting [12].

Once the words and their classes are known, one may continue by parsing the
words in the sentence. The dapipe tool uses the Universal Dependencies
schema333See https://universaldependencies.org/introduction.html for this,
finding the main grammatical relations between each word to construct a
dependency tree that links every word in the sentence together.

All these tasks are achieved by using UDPipe [13] based on the universal
dependencies [14] data for Danish. Sample output follows.

    leon@blade:~/dapipe$ ./dapipe dktest.txt

    ...

    # text = Derfor presser EU-siden hrdt p, for at amerikanerne ikke saboterer aftalen, der
    blandt andet lader udenlandske observatrer inspicere iranske atomanlg.
    1ΨDerforΨderforΨADVΨ_Ψ_Ψ2ΨadvmodΨ_Ψ_
    2ΨpresserΨpresseΨVERBΨ_ΨMood=Ind|Tense=Pres|VerbForm=Fin|Voice=ActΨ0ΨrootΨ_Ψ_
    3ΨEU-sidenΨEU-sidΨNOUNΨ_ΨDefinite=Def|Gender=Com|Number=SingΨ2ΨnsubjΨ_Ψ_
    4ΨhrdtΨhrdtΨADVΨ_ΨDegree=PosΨ2ΨadvmodΨ_Ψ_
    5ΨpΨpΨADPΨ_ΨAdpType=PrepΨ4ΨcaseΨ_ΨSpaceAfter=No
    6Ψ,Ψ,ΨPUNCTΨ_Ψ_Ψ2ΨpunctΨ_Ψ_
    7ΨforΨforΨADPΨ_ΨAdpType=PrepΨ11ΨmarkΨ_Ψ_
    8ΨatΨatΨSCONJΨ_Ψ_Ψ11ΨmarkΨ_Ψ_
    9ΨamerikanerneΨamerikanerΨNOUNΨ_ΨDefinite=Def|Gender=Com|Number=PlurΨ11ΨnsubjΨ_Ψ_
    10ΨikkeΨikkeΨADVΨ_Ψ_Ψ11ΨadvmodΨ_Ψ_
    11ΨsabotererΨsabotereΨVERBΨ_ΨMood=Ind|Tense=Pres|VerbForm=Fin|Voice=ActΨ2ΨadvclΨ_Ψ_
    12ΨaftalenΨaftaleΨNOUNΨ_ΨDefinite=Def|Gender=Com|Number=SingΨ11ΨobjΨ_ΨSpaceAfter=No
    13Ψ,Ψ,ΨPUNCTΨ_Ψ_Ψ12ΨpunctΨ_Ψ_
    14ΨderΨderΨPRONΨ_ΨPartType=InfΨ17ΨnsubjΨ_Ψ_
    15ΨblandtΨblandtΨADPΨ_ΨAdpType=PrepΨ17ΨadvmodΨ_Ψ_
    16ΨandetΨandenΨPRONΨ_ΨGender=Neut|Number=Sing|PronType=IndΨ15ΨfixedΨ_Ψ_
    17ΨladerΨladeΨVERBΨ_ΨMood=Ind|Tense=Pres|VerbForm=Fin|Voice=ActΨ12Ψacl:relclΨ_Ψ_
    18ΨudenlandskeΨudenlandskΨADJΨ_ΨDegree=Pos|Number=PlurΨ19ΨamodΨ_Ψ_
    19ΨobservatrerΨobservatrΨNOUNΨ_ΨDefinite=Ind|Gender=Com|Number=PlurΨ17ΨobjΨ_Ψ_
    20ΨinspicereΨinspicereΨADVΨ_ΨDegree=CmpΨ21ΨadvmodΨ_Ψ_
    21ΨiranskeΨiranskΨADJΨ_ΨDegree=Pos|Number=PlurΨ22ΨamodΨ_Ψ_
    22ΨatomanlgΨatomanlgΨNOUNΨ_ΨDefinite=Ind|Gender=Neut|Number=PlurΨ17ΨobjΨ_ΨSpaceAfter=No
    23Ψ.Ψ.ΨPUNCTΨ_Ψ_Ψ2ΨpunctΨ_ΨSpacesAfter=\n

##  4 Summary

This note describes a free toolkit for basic natural language processing. The
two tools daner and dapipe are built on freely-available resources and
maintaned at ITU Copenhagen.

Feedback, requests, and so on can be addressed to the ITU NLP research group:
http://nlp.itu.dk/

## References

  * [1] Christopher Manning, Mihai Surdeanu, John Bauer, Jenny Finkel, Steven Bethard, and David McClosky.  The Stanford CoreNLP natural language processing toolkit.  In Proceedings of 52nd annual meeting of the Association for Computational Linguistics, pages 55–60. Association for Computational Linguistics, 2014.
  * [2] Jenny Rose Finkel, Trond Grenager, and Christopher Manning.  Incorporating non-local information into information extraction systems by Gibbs sampling.  In Proceedings of the 43rd annual meeting on association for computational linguistics, pages 363–370. Association for Computational Linguistics, 2005.
  * [3] Leon Derczynski and Kenneth S Bøgh.  DKIE: Open source information extraction for Danish.  In Proceedings of the Demonstrations at the 14th Conference of the European Chapter of the Association for Computational Linguistics, pages 61–64. Association for Computational Linguistics, 2014.
  * [4] Hamish Cunningham, Diana Maynard, Kalina Bontcheva, Valentin Tablan, Niraj Aswani, Ian Roberts, Genevieve Gorrell, Adam Funk, Angus Roberts, Danica Damljanovic, Leon Derczynski, et al.  Developing language processing components with GATE version 8.0.  University of Sheffield Department of Computer Science, 2012.
  * [5] Anders Johannsen, Héctor Martínez Alonso, and Barbara Plank.  Universal dependencies for danish.  In Proc. International Workshop on Treebanks and Linguistic Theories, page 157, 2015.
  * [6] Matthias T Kromann, Line Mikkelsen, and Stine Kern Lynge.  Danish dependency treebank.  In Proc. International Workshop on Treebanks and Linguistic Theories, pages 217–220, 2003.
  * [7] Leon Derczynski and Sean Chester.  Generalised brown clustering and roll-up feature generation.  In Proc. AAAI Conference on Artificial Intelligence, 2016.
  * [8] Leon Derczynski, Sean Chester, and Kenneth S Bøgh.  Tune your brown clustering, please.  In Proc. International Conference Recent Advances in Natural Language Processing, RANLP, volume 2015, pages 110–117. Association for Computational Linguistics, 2015.
  * [9] Armand Joulin, Edouard Grave, Piotr Bojanowski, and Tomas Mikolov.  Bag of tricks for efficient text classification.  arXiv preprint arXiv:1607.01759, 2016.
  * [10] Barbara Plank and Željko Agić.  Distant supervision from disparate sources for low-resource part-of-speech tagging.  In Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing, pages 614–620. Association for Computational Linguistics, 2018.
  * [11] Barbara Plank, Anders Søgaard, and Yoav Goldberg.  Multilingual Part-of-Speech Tagging with Bidirectional Long Short-Term Memory Models and Auxiliary Loss.  In ACL 2016, arXiv preprint arXiv:1604.05529, 2016.
  * [12] Leon Derczynski, Diana Maynard, Giuseppe Rizzo, Marieke Van Erp, Genevieve Gorrell, Raphaël Troncy, Johann Petrak, and Kalina Bontcheva.  Analysis of named entity recognition and linking for tweets.  Information Processing & Management, 51(2):32–49, 2015.
  * [13] Milan Straka, Jan Hajic, and Jana Straková.  UDPipe: Trainable Pipeline for Processing CoNLL-U Files Performing Tokenization, Morphological Analysis, POS Tagging and Parsing.  In Proc. LREC, 2016.
  * [14] Joakim Nivre, Marie-Catherine De Marneffe, Filip Ginter, Yoav Goldberg, Jan Hajic, Christopher D Manning, Ryan McDonald, Slav Petrov, Sampo Pyysalo, Natalia Silveira, et al.  Universal dependencies v1: A multilingual treebank collection.  In Proceedings of the Tenth International Conference on Language Resources and Evaluation (LREC 2016), pages 1659–1666, 2016.