@@ -8,76 +8,77 @@ This example contains a contaminant transported by a base groundwater flow in a

![](./img/setup.png)

## Model description

Two different models are applied to simulate the system: In a first step, the groundwater velocity is evaluated under stationary conditions. Therefore the single phase model is applied.

In a second step, the contaminant gets transported based on the groundwater velocity field. It is assumed, that the dissolved contaminant does not affect density and viscosity of the groundwater and thus, it is handled as a tracer by the tracer model. The tracer model is then solved instationarily.

Two different models are applied to simulate the system: In a first step, the groundwater velocity is evaluated under stationary conditions using the single phase model.

In a second step, the contaminant is transported with the groundwater velocity field. It is assumed, that the dissolved contaminant does not affect density and viscosity of the groundwater, and thus, it is handled as a tracer by the tracer model. The tracer model is then solved instationarily.

### 1p Model

The single phase model uses Darcy's law as the equation for the momentum conservation:

```math

\textbf v = - \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho {\textbf g} \right)

\textbf v = - \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho {\textbf g} \right),

```

With the darcy velocity $` \textbf v `$, the permeability $` \textbf K`$, the dynamic viscosity $` \mu`$, the pressure $`p`$, the density $`\rho`$ and the gravity $`\textbf g`$.

with the darcy velocity $` \textbf v `$, the permeability $` \textbf K`$, the dynamic viscosity $` \mu`$, the pressure $`p`$, the density $`\rho`$ and the gravity $`\textbf g`$.

Darcy's law is inserted into the continuity equation:

Darcy's law is inserted into the mass balance equation:

```math

\phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0

```

with porosity $`\phi`$.

where $`\phi`$ is the porosity.

The equation is discretized using a cell-centered finite volume scheme as spatial discretization for the pressure as primary variable. For details on the discretization scheme, have a look at the dumux [handbook](https://dumux.org/handbook).

### Tracer Model

The transport of the contaminant component $`\kappa`$ is based on the previously evaluated velocity field $`\textbf v`$ with the help the following mass balance equation:

The transport of the contaminant component $`\kappa`$ is based on the previously evaluated velocity field $`\textbf v`$ with the help of the following mass balance equation:

With the porosity $`\phi`$, the mass fraction of the contaminant component $`\kappa`$: $`X^\kappa`$, the porous medium diffusivity $` D^\kappa_\text{pm} `$ and the density of the fluid phase $`\varrho`$.

where $`X^\kappa`$ is the mass fraction of the contaminant component $`\kappa`$ and $` D^\kappa_\text{pm} `$ is the effective diffusivity.

The porous medium diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$, the porosity $`\phi`$ and the porous medium tortuosity $`\tau`$ by the following equation:

The effective diffusivity is a function of the diffusion coefficient of the component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous medium:

```math

D^\kappa_\text{pm}= \phi \tau D^\kappa

D^\kappa_\text{pm}= \phi \tau D^\kappa.

```

The primary variable of this model is the mass fraction $`X^\kappa`$. We apply the same spatial discretization as in the single phase model and use the implicit Euler method for time discretization. For more information, have a look at the dumux handbook.

The primary variable of this model is the mass fraction $`X^\kappa`$. We apply the same spatial discretization as in the single phase model and use the implicit Euler method for time discretization. For more information, have a look at the dumux [handbook](https://dumux.org/handbook).

In the following, we take a close look at the files containing the set-up: At first, boundary conditions and spatially distributed parameters are set in `problem_1p.hh` and `spatialparams_1p.hh`, respectively, for the single phase model and subsequently in `problem_tracer.hh` and `spatialparams_tracer.hh` for the tracer model. Afterwards, we show the different steps for solving the model in the source file `main.cc`. At the end, we show some simulation results.

In the following, we take a close look at the files containing the set-up: The boundary conditions and spatially distributed parameters for the single phase model are set in `problem_1p.hh` and `spatialparams_1p.hh`.

For the tracer model, this is done in the files `problem_tracer.hh` and `spatialparams_tracer.hh`, respectively. Afterwards, we show the different steps for solving the model in the source file `main.cc`. Finally, some simulation results are shown.

## The file `spatialparams_1p.hh`

In this file, we generate the random permeability field in the constructor of the `OnePTestSpatialParams` class. Thereafter, spatial properties of the porous medium such as the permeability and the porosity are defined in various functions for the 1p problem.

We want to generate a random permeability field. For this, we use a random number generation of the C++ standard library.

In this file, we generate a random permeability field in the constructor of the `OnePTestSpatialParams` class.

For this, we use the random number generation facilities provided by the C++ standard library.

```cpp

#include<random>

```

In the file `properties.hh` all properties are declared.

We use the properties for porous medium flow models, declared in the file `properties.hh`.

```cpp

#include<dumux/porousmediumflow/properties.hh>

```

We include the spatial parameters for single-phase, finite volumes from which we will inherit.

We include the spatial parameters class for single-phase models discretized by finite volume schemes.

The spatial parameters defined for this example will inherit from those.

```cpp

#include<dumux/material/spatialparams/fv1p.hh>

namespaceDumux{

```

In the `OnePTestSpatialParams` class, we define all functions needed to describe the porous matrix, e.g. porosity and permeability for the 1p_problem.

In the `OnePTestSpatialParams` class, we define all functions needed to describe the porous medium, e.g. porosity and permeability for the 1p_problem.

```cpp

template<classGridGeometry,classScalar>

classOnePTestSpatialParams

:publicFVSpatialParamsOneP<GridGeometry,Scalar,

OnePTestSpatialParams<GridGeometry,Scalar>>

{

```

We introduce `using` declarations that are derived from the property system, which we need in this class.

We declare aliases for types that we are going to need in this class.

We generate random fields for the permeability using a lognormal distribution, with the `permeability_` as mean value and 10 % of it as standard deviation. A seperate distribution is used for the lens using `permeabilityLens_`.

We generate random fields for the permeability using lognormal distributions, with `permeability_` as mean value and 10 % of it as standard deviation.

A separate distribution is used for the lens using `permeabilityLens_`. A permeability value is created for each element of the grid and is stored in the vector `K_`.

@@ -120,7 +126,9 @@ We generate random fields for the permeability using a lognormal distribution, w

}

```

### Properties of the porous matrix

We define the (intrinsic) permeability $`[m^2]`$ using the generated random permeability field. In this test, we use element-wise distributed permeabilities.

This function returns the permeability $`[m^2]`$ to be used within a sub-control volume (`scv`) inside the element `element`.

One can define the permeability as function of the primary variables on the element, which are given in the provided `ElementSolution`.

Here, we use element-wise distributed permeabilities that were randomly generated in the constructor (see above).

The top and bottom of our domain are Neumann boundaries:

```cpp

else

values.setAllNeumann();

returnvalues;

...

...

@@ -363,7 +372,7 @@ we retreive again the global position

constauto&pos=scvf.ipGlobal();

PrimaryVariablesvalues(0);

```

we assign pressure values in [Pa] according to a pressure gradient to 1e5 Pa at the top and 1.1e5 Pa at the bottom.

and assign pressure values in [Pa] according to a pressure gradient to 1e5 Pa at the top and 1.1e5 Pa at the bottom.

```cpp

values[0]=1.0e+5*(1.1-pos[dimWorld-1]*0.1);

returnvalues;

...

...

@@ -391,12 +400,13 @@ We leave the namespace Dumux.

## The file `spatialparams_tracer.hh`

In this file, we define spatial properties of the porous medium such as the permeability and the porosity in various functions for the tracer problem. Further, spatial dependent properties of the tracer fluid system are defined and in the end two functions handel the calculated volume fluxes from the solution of the 1p problem.

In the file `properties.hh`, all properties are declared.

In this file, we define spatial properties of the porous medium such as the permeability and the porosity in various functions for the tracer problem.

Furthermore, spatial dependent properties of the tracer fluid system are defined and in the end two functions handle the calculated volume fluxes from the solution of the 1p problem.

We use the properties for porous medium flow models, declared in the file `properties.hh`.

```cpp

#include<dumux/porousmediumflow/properties.hh>

```

As in the 1p spatialparams, we inherit from the spatial parameters for single-phase, finite volumes, which we include here.

As in the 1p spatialparams, we inherit from the spatial parameters for single-phase models using finite volumes, which we include here.

```cpp

#include<dumux/material/spatialparams/fv1p.hh>

```

...

...

@@ -406,7 +416,6 @@ namespace Dumux {

```

In the `TracerTestSpatialParams` class, we define all functions needed to describe spatially dependent parameters for the `tracer_problem`.

```cpp

template<classGridGeometry,classScalar>

classTracerTestSpatialParams

:publicFVSpatialParamsOneP<GridGeometry,Scalar,

...

...

@@ -443,24 +452,30 @@ We do not consider dispersivity for the tracer transport. So we set the dispersi

{return0;}

```

### Properties of the fluid system

In the following, we define fluid properties that are spatial parameters in the tracer model. They can possible vary with space but are usually constants. Further spatially constant values of the fluid system are defined in the `TracerFluidSystem` class in `problem.hh`.

In the following, we define fluid properties that are spatial parameters in the tracer model.

They can possible vary in space but are usually constants.

Furthermore, spatially constant values of the fluid system are defined in the `TracerFluidSystem` class in `problem.hh`.

We define the fluid density to a constant value of 1000 $`\frac{kg}{m^3}`$.

```cpp

ScalarfluidDensity(constElement&element,

constSubControlVolume&scv)const

{return1000;}

```

We define the fluid molar mass.

This interface defines the fluid molar mass within the sub-control volume `scv`.

We define a function which returns the field of volume fluxes. This is e.g. used to calculate the transport of the tracer.

We define a function which returns the volume flux across the given sub-control volume face `scvf`.

This flux is obtained from the vector `volumeFlux_` that contains the fluxes across al sub-control volume faces of the discretization.

This vector can be set using the `setVolumeFlux` function.

```cpp

template<classElementVolumeVariables>

ScalarvolumeFlux(constElement&element,

...

...

@@ -471,7 +486,8 @@ We define a function which returns the field of volume fluxes. This is e.g. used

returnvolumeFlux_[scvf.index()];

}

```

We define a function to set the volume flux. This is used in the main function to set the volume flux to the calculated value based on the solution of the 1p problem.

We define a function that allows setting the volume fluxes for all sub-control volume faces of the discretization.

This is used in the main function after these fluxes have been based on the pressure solution obtained with the single-phase model.

```cpp

voidsetVolumeFlux(conststd::vector<Scalar>&f)

{volumeFlux_=f;}

...

...

@@ -491,7 +507,7 @@ private:

Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties.

### Include files

Again, we have to include the dune grid interface:

Again, we use YaspGrid, the implementation of the dune grid interface for structured grids:

```cpp

#include<dune/grid/yaspgrid.hh>

```

...

...

@@ -507,7 +523,7 @@ We include again the porous medium problem class that this class is derived from

```cpp

#include<dumux/porousmediumflow/problem.hh>

```

and the base fluidsystem

and the base fluidsystem. We will define a custom fluid system that inherits from that class.

```cpp

#include<dumux/material/fluidsystems/base.hh>

```

...

...

@@ -561,30 +577,34 @@ We define the spatial parameters for our tracer simulation:

@@ -798,18 +823,22 @@ and another one for in- and output.

```cpp

#include<iostream>

```

Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and linear solvers. So we need some includes from that.

Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and linear solvers.

Here, we include classes related to parallel computations, time measurements and file I/O.

In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file `properties.hh`.

In Dumux, the property system is used to specify classes and compile-time options to be used by the model.

For this, different properties are defined containing type definitions, values and methods.

All properties are declared in the file `properties.hh`.

```cpp

#include<dumux/common/properties.hh>

```

The following file contains the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.

The following file contains the parameter class, which manages the definition and retrieval of input

parameters by a default value, the inputfile or the command line.

```cpp

#include<dumux/common/parameters.hh>

```

...

...

@@ -843,7 +872,7 @@ int main(int argc, char** argv) try

{

usingnamespaceDumux;

```

We define the type tags for the two problems. They are created in the individual problem files.

Convenience aliases for the type tags of the two problems, which are defined in the individual problem files.

@@ -862,7 +891,11 @@ We parse the command line arguments.

Parameters::init(argc,argv);

```

### Create the grid

A gridmanager tries to create the grid either from a grid file or the input file. We solve both problems on the same grid. Hence, the grid is only created once using the type tag of the 1p problem.

The `GridManager` class creates the grid from information given in the input file.

This can either be a grid file, or in the case of structured grids, by specifying the coordinates

of the corners of the grid and the number of cells to be used to discretize each spatial direction.

Here, we solve both the single-phase and the tracer problem on the same grid.

Hence, the grid is only created once using the grid type defined by the type tag of the 1p problem.

In the following section, we set up and solve the 1p problem. As the result of this problem, we obtain the pressure distribution in the problem domain.

In the following section, we set up and solve the 1p problem. As the result of this problem, we obtain the pressure distribution in the domain.

#### Set-up

We create and initialize the finite volume grid geometry, the problem, the linear system, including the jacobian matrix, the residual and the solution vector and the gridvariables.

We need the finite volume geometry to build up the subcontrolvolumes (scv) and subcontrolvolume faces (scvf) for each element of the grid partition.

...

...

@@ -895,14 +928,14 @@ The jacobian matrix (`A`), the solution vector (`p`) and the residual (`r`) are

autoA=std::make_shared<JacobianMatrix>();

autor=std::make_shared<SolutionVector>();

```

The grid variables store variables on scv and scvf (volume and flux variables).

The grid variables store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables).

We calculate the volume flux for every subcontrolvolume face, which is not on a Neumann boundary (is not on the boundary or is on a Dirichlet boundary).

We calculate the volume fluxes for all sub-controlvolume faces except for Neumann boundary faces