Let $\Omega$ be the computational domain in ${\rm I\kern-1.69998ptR}^{d}$, $d=1,2$, or $3$, with boundary $\partial\Omega$ and let $J=(0,T]$ be some finite interval with $T$ denoting the final time. The diffusive-advective transport of a contaminant in soil can be mathematically described by

where $c(t,{\bf x})$ denotes the concentration of the contaminant and ${\bf D}$ is the diffusion-dispersion tensor. The water content $\Theta$ and the water flux ${\bf Q}$ are determinated by solving the Richards’ equation, e.g. (4). Dirichlet, Neumann and/or Flux boundary conditions complete the model. For a more general model, including multicomponent reactive transport with equilibrium or non-equilibrium sorption we refer to [30].

The water flux ${\bf Q}(t,{\bf x})$ appearing in (1), as well as the water content $\Theta(\psi)$, are obtained by solving the mass balance equation for water, which is assumed incompressible

and the Darcy’s law

together with initial $\displaystyle\psi(t=0,{\bf x})=\Psi_{I}({\bf x})$ and boundary conditions. Above, $\psi(t,{\bf x})$ is the pressure head and $K$ the hydraulic conductivity. We neglected the influence of gravity. Combining the equations (2) and (3) one obtains the Richards equation

which is a typical mathematical model for water flow through saturated/unsaturated soil. For the coefficient functions $\displaystyle\Theta(\cdot)$ and $\displaystyle K(\cdot)$ functional dependencies of the pressure are assumed, e.g. van Genuchten-Mualem, so that the unknowns in (2) – (3) are reduced to two. Nevertheless, in this work we consider only saturated flow, e.g. $\Theta(\psi)=\Theta_{S}$ and $K(\Theta(\psi))=K_{S}$ are constants. We also restrict ourselves to two dimensional simulations ($d=2$).

The following numerical schemes are presented for homogeneous Dirichlet boundary conditions for both flow and transport problems. This is just for the sake of simplicity and does not affect in any way the generality of the schemes.

We present here a classical MHFEM scheme based on the lowest order Raviart-Thomas finite elements and two new variants of it to simulate solute transport in porous media. For higher order MHFEM schemes we refer to [10]. In what follows, $L^{2}(\Omega)$ is the space of square integrable functions on $\Omega$ and $\langle\cdot,\cdot\rangle$ for the inner product on $L^{2}(\Omega)$. The functions in $H(\mathrm{div};\Omega)$ are vector valued, having a $L^{2}$ divergence. For the discretization in time we let $N\in\mathbb{N}$ be strictly positive, and define the time step $\tau={T/N}$, as well as $t_{n}=n\tau$ $(n\in\{1,2,\dots,N\})$. Furthermore, we let ${\mathcal{T}_{h}}$ be a regular decomposition of $\Omega$ into closed $d$-simplices; $h$ stands for the mesh size. We denote by $\mathcal{S}_{h}=\mathcal{S}_{h}^{I}\cup\mathcal{S}_{h}^{D}$ the set of all faces of ${\mathcal{T}}_{h}$, where $\mathcal{S}_{h}^{I}$ are interior faces and $\mathcal{S}_{h}^{D}$ are faces on the boundary. We always denote by ${\bf n}$ the outer normal. We use the discrete subspaces $W_{h}\subset L^{2}(\Omega)$ and $V_{h}\subset H(\mathrm{div};\Omega)$ defined as

In other words, $W_{h}$ denotes the space of piecewise constant functions, while $V_{h}$ is the lowest order Raviart-Thomas space (see [11]).

The scheme for solving the water flow (2)-(3) is based on MFEM and Euler implicit and reads

To obtain a mixed formulation for the transport equation (1) we first introduce the flux variable $\bf{q}$ as additional unknown:

The fully discrete mixed approximation of (8) now reads as follows:

Unfortunately, the resulting equations (9)–(10) lead to a linear system of equations with an indefinite matrix such that standard iterative solvers cannot be applied. To overcome this difficulty, we use a hybridization technique (see [11]). Its basic idea is to relax firstly the continuity constraint of the normal components of the fluxes over inter-element faces that is implied by ${\bf v}\in H(\mathrm{div};\Omega)$. The continuity constraint is then ensured by means of an additional variational equation involving Lagrange multipliers. Precisely, the space $V_{h}$ is replaced by

The discrete space for the Lagrange multiplier is defined by

The fully-discrete mixed hybrid variational formulation of the overall system (9)–(10) then reads as follows:

The hybridization increases first the complexity of the nonlinear systems. This drawback can be resolved by applying static condensation (see [11]) and eliminate now the internal degrees of freedom. Moreover, the Lagrange multipliers can be used to construct a second order accurate approximation of the primary variables [11, 8]. The last observation also stays at the basis of the development of the two new schemes presented at the end of this section. The new schemes are much more robust for convection-dominated problems.

Let now $\Theta_{T}^{n}$ and $c_{T}^{n}$ denote the water content and concentration, respectively, on the element $T$, $\{q_{TS}^{n}\}_{S\subset T}$, $\{Q_{TS}^{n}\}_{S\subset T}$ the components of the flux of contaminant and water, respectively, in the local Raviart-Thomas space basis $\{{\bf{w}}_{TS}\}_{S\subset T}({\bf{x}})=\dfrac{{\bf{x}}-{\bf{x}}_{S}}{d\,|T|}$ (cf. [11]), $\displaystyle B_{TSS^{\prime}}:=\int_{T}({\bf{\bf D}^{-1}}{\bf{w}}_{TS^{\prime}})\cdot{\bf{w}}_{TS}\,d{\bf{x}}$, and $\lambda_{S}^{n}$ be the constant Lagrange multiplier on the face $S$. We denote by $|T|:=\int_{T}d{\bf{x}}$ the area (volume) of element $T$ and by ${\bf{x}}_{S}$ the corner of $T$ which is not on $S$. We obtain from (11)–(13) the following system of nonlinear equations:

Mass conservation equation:

Equation for the flux:

Continuity of the flux over faces:

The system of equations (14)–(16) gives the classical MHFEM scheme (referred to as MHFEM 0 in Sec. 3) for solute transport in porous media.

A new MHFEM scheme is obtained by using the Lagrange multipliers, instead of the piecewise constant concentrations, for discretizing the convective term in (1). Instead of (15), this yields

The equations (14) and (16) remain unchanged. The new MHFEM scheme (referred to as MHFEM 1 in Sec. 3) is given by (14), (17) and (16). Numerical results show that the new scheme is much more robust for convection-dominated problems and has the same convergence properties as the classical scheme. A rigorous analysis of this new scheme is ongoing [33].

We also propose a full upwind MHFEM scheme (referred to as MHFEM 2 in Sec. 3) given by (14), (16) and

with

The upwind MHFEM scheme is suitable for highly convection-dominated-problems. The above schemes are, up to our knowledge, not reported so far in the literature.

The three schemes were implemented in the UG toolbox (see [6]). The linear systems are solved by a multigrid method.

For a description of the general case see [14] and for a more detailed description of Galerkin finite elements see e.g. [19].

For the space discretization we use the subspace $V_{h}^{GFEM}\subset H^{1}(\Omega)$

with $k=1,2$ where $\mathcal{P}_{k}(T)$ denotes the space of all polynomials up to the degree $k$. We restrict ourselves to Lagrange elements, i.e., all degrees of freedom are function values in fixed points. This points are called nodes. Let the nodes ${\bf a}_{i}$ ($i=1,\dots,M$) be sorted such that for $i=1,\dots,M_{1}$ ${\bf a}_{i}$ is an inner node. The basis function associated with the node ${\bf a}_{i}$ is denoted by $\varphi_{i}$.

For the time discretization the BE method or the BDF(2) method is used, i.e., the time derivative is replaced by the difference quotient

At the first time step ($n=1$) the backward Euler method is applied because there is no value $u^{-1}$.

Differentiating out in (1) and using (2) leads to

So the fully discrete formulation of the transport problem is in case of the backward Euler method:

and in case of the BDF(2) method:

All integrals without space derivatives are assembled using mass lumping. These integrals are evaluated using a nodal quadrature formula, i.e. a quadrature rule with the quadrature points ${\bf a}_{i}$ and weights

This is not possible for standard quadratic elements because some weights would be zero. Instead of that quadratic elements enriched by a bubble function (named $\overline{P}_{2}$) can be used (see [13]).

Let $c_{i}^{n}$ be the concentration value in the node ${\bf a}_{i}$ at the time $t_{n}$. From (22) and (24) respectively we obtain the following system of equations

where the difference quotient $D_{n}$ is $D^{BE}_{n}$ in case of (22) and $D^{BDF2}_{n}$ in case of (24). The integrals are evaluated with an appropriate quadrature rule.

For the convection dominated case a finite volume stabilization for linear elements is implemented. Let $\Omega_{k}$ denote the control volume (Donald-Diagram) associated to the node ${\bf a}_{k}$. The integral stemming from the convective term is approximated by

Another point of view is that we carry out a finite volume discretization for the whole partial differential equation. Let us consider the finite volume scheme treated in [19, chapter 6.2]. In this FV scheme, the discretization of all terms except the advective coincide with that one of the linear finite element method using mass lumping. Hence we only have to replace the advective term in the linear finite element discretization to get a finite volume scheme. The boundary integral on the right hand side of (27) can be discretized using standard finite volume techniques. We have implemented full upwinding, exponential upwinding and partial upwinding (see e.g. [19, chapter 6.2]).

The Galerkin finite elements are used only for the transport problem (1). The flow problem (2)-(3) is solved using mixed hybrid finite elements (see sec. 2.2).

The model was implemented using a software kernel for parallel computations in the field of PDEs called M++ (see [42]).

In this section we briefly present the classical finite volume method and a generalization of the method to higher order trial spaces. In order to keep the presentation simple, we restrict ourselves to two dimensions only. However, generalizations to higher dimensions are straightforward. The classic FVM can be found in e.g. [5], a second order FVM in [22] and recently an arbitrary order FVM has been described in [39] and [40].

The FVM uses a reformulation of the equations (1) and (2). For an arbitrary control volume $B\subset\Omega$ the two equations

express the conservation property of the system. In order to discretize these equations we will first choose an appropriate trial space for the unknown functions are $c(t,{\bf x})$ and $\psi(t,{\bf x})$ and then describe our choice of control volumes $B$.

By $\mathcal{P}_{k}$ we denote again the space of polynomials of degree at most $k$. $C^{0}(\Omega)$ is the space of continuous functions on $\Omega$ and we will use a trial function space for triangles defined by

An appropriate choice of basis functions are given by Lagrange interpolants for interpolation nodes on the triangles. Therefore, for each triangle $T\in\mathcal{T}_{h}$ denote its corners by $\left\{\mathbf{a}_{0},\mathbf{a}_{1},\mathbf{a}_{2}\right\}$ ($\mathbf{a}_{i}\in\mathbb{R}^{2}$) and use barycentric coordinates $\left(\lambda_{0},\lambda_{1},\lambda_{2}\right)$ in order to describe any point ${\bf x}\in T$, i.e., one has a representation

On each triangle we have a set of $\binom{k+2}{2}$ nodes for a given degree $k$ defined to be the points $n_{i_{1},i_{2}}$ of $T$ with barycentric coordinates

We introduce a second decomposition $\mathcal{B}_{h}(\Omega):=\left\{B_{0},B_{1},B_{2},\dots\right\}$ of the domain $\Omega$ into (open) control volumes satisfying

The shape of the control volume $B_{i}$ will be restricted to polygons, though this is not necessary in general. Furthermore, we will construct the control volumes the way that one and only one degree of freedom of the trial space will be contained in one control volume.

In the case of a linear trial space the degrees of freedom are only located in the corners of the triangles and the control volumes can be constructed by connecting the barycenter of every triangle $T\in\mathcal{T}_{h}$ by straight lines with each edge midpoint. Hence, every triangle is partitioned into three parts (subcontrol volumes) and every part is associated with one node of the triangle. Now the control volume $B_{i}$ containing the node $i$ is defined by the union of all subcontrol volumes being associated with this node. A visualization is given in Fig. 1a).

For the higher order trial spaces we can reduce the construction of the control volumes to the procedure of the linear case. Indeed, we subdivide the triangles into a set of smaller subtriangles using the isolines

On all resulting subtriangles, nodes are located in the corners only so that the structure of the subtriangles is equivalent to the structure of triangles with linear trial functions. Now we use the construction described above on these finer triangles. An example for a quadratic triangle is given in Fig. 1b) and for a set of triangles the control volumes are shown in Fig. 1c).

Now we can set up the FVM schemes using BE time stepping:

The BDF(2) time stepping scheme is obtained by replacing the time discretization $D_{n}^{BE}$ with $D_{n}^{BDF2}$ as described in the previous section. The FVM-schemes of arbitrary order are implemented in the UG toolbox (see [6]). In Sec. 3 we denote by FV $i$ the FVM scheme of order $i$.

The GRW algorithm is a generalization of the particle tracking (PT) method which increases the speed of the computations and considerably improves the accuracy of the numerical simulations [38]. For a two-dimensional problem and constant porosity, the solution of a parabolic equation of form (1) is described using $N$ particles which move in a grid, undergoing advective displacements and diffusive jumps according to the random walk law. The concentration field at a given time $t=k\delta t$ and a point $(x_{1},x_{2})=(i_{1}\delta x_{1},i_{2}\delta x_{2})$ is given by

where $\Delta_{l}=2s_{l}\delta x_{l}$, $l=1,2$, are the lengths of the symmetrical intervals centered at $x_{l}$ and $n(i_{1},i_{2},k)$ is the number of particles which at time step $k$ lie at the grid point $(i_{1},i_{2})$.

The one-dimensional GRW algorithm describes the scattering of the $n(i,k)$ particles from $(x_{i},t_{k})$ by

where $v_{j}=Q_{j}\delta t/\delta x$ are discrete displacements produced by the velocity field and $d$ describes the diffusive jumps. The quantities $\delta n$ are Bernoulli random variables and describe respectively, the number of particles which remain at the same grid site after an advective displacement and the number of particles jumping to the left and to the right of the advected position $j+v_{j}$. The distribution of the particles at the next time $(k+1)\delta t$ is given by

The average number of particles undergoing diffusive jumps and the average number of particles remaining at the same node after the displacement $v_{j}$ are given by the relations

where $0\leq r\leq 1$. The diffusion coefficient $D$ is related to the grid steps by the relation

For two and three-dimensional cases, the same procedure is repeated for all space directions.

Because the total number of particles $N$ contained in the grid is conserved, the GRW algorithm is stable. The condition $r\leq 1$, ensures that there is no numerical diffusion. In [38] it was shown that for Gaussian diffusion the numerical solution converges as $\mathcal{O}(\delta x^{2})$ $+\mathcal{O}(N^{-1/2})$, i.e. for large numbers of particles the convergence order is $\mathcal{O}(\delta x^{2})$, the same as for the finite differences scheme.

The “reduced fluctuations” GRW algorithm is defined by

where $n=n(j,k)-\delta n(j,j+v_{j},k)$, $[n/2]$ is the integer part of $n/2$ and $\theta$ is a variable taking the values $0$ and $1$ with probability $1/2$. This algorithm is appropriate for large scale problems, for two reasons. Firstly, the diffusion front does not extend beyond the limit concentration defined by one particle at a grid point, keeping a physical significant shape (unlike in finite differences where a pure diffusion front has a cubic shape of side $\sim(2Dt)^{1/2}$). Secondly, the “reduced fluctuations” algorithm requires only a minimum number of calls of the random number generator.

Let us compare GRW with the classical PT scheme. PT approximates the solution of (1) by (38), after having generated (sequentially) $N$ particles’ trajectories and estimating (through a post-processing) the number of occurrences $n$ of the $N$ particles in sampling volumes $\Delta_{l}$. PT uses $N$ particles and $T$ random numbers to generate a single particle trajectory. To do the same job, GRW needs a single call of the uniformly distributed random numbers generator to move all $n$ particles from at the most as many points as the grid size $L$, over $T$ time steps. Hence, (GRW-CPU time) / (PT-CPU time)=$\mathcal{O}(L/N)$. For instance, if $N=10^{24}$ (Avogadro’s number) and $L=10^{6}$ grid points, a huge speed-up of computations by a factor of $10^{18}$ is achieved. A comparison with a PT code (diffusion over ten time steps of $N$ particle starting at the center of a cubic grid) shows that while for the GRW algorithm there were practically no limitations concerning the total number of particles and the computation time was of about one second, PT simulations for $N=10^{9}$ particles already required a computing time of about one hour and 256 processors on a CRAY T3E parallel machine [38].

To compute the variance of particle displacements, $s_{ll}^{2}$, $l=1,2$, a more accurate result is obtained if instead of the concentration (38) one uses the point density of the number of particles $n(i_{1},i_{2},k)$:

Using (39), the effective coefficients are computed as

Let us consider $N_{{}_{X_{0}}}$ points uniformly distributed inside the initial plume, $N/N_{{}_{X_{0}}}$ particles at each initial point and let $n(i_{1},i_{2},k;i_{01},i_{02})$ be the distribution of particles at the time step $k$ given by the GRW procedure for a diffusion process starting at $(i_{01}\delta x_{1},i_{02}\delta x_{2})$. Writing the distribution for the extended plume as

the averages from (39) can be rewritten in the form

where $\alpha$ stands for $i_{l}$ and $i_{l}^{2}$ respectively. It follows from (40) that the variance (39) is an average over the trajectories of the diffusion process starting at given initial positions and over the distribution of the initial positions.

In this section we consider a deterministic problem with a known analytical solution. The problem represents the idealized case of solute transport without sorption or reaction in a heterogeneous soil. We quantify the numerical diffusion for the various schemes. Additionally, we study the influence of the Péclet number on numerical diffusion. Consequently, we solve the equation (1) with $\Theta=1$ and a constant water flux ${\bf Q}=\left(\begin{array}[]{c}Q_{x}\\ Q_{y}\end{array}\right)$. The computational domain is $[0,10]\times[0,10]$. The diffusion-dispersion coefficient is given by ${\bf D}=D\left(\begin{array}[]{cc}1&0\\ 0&1\end{array}\right)$, with $D$ constant. We consider adequate initial and Dirichlet boundary conditions such that equation (1) admits the analytical solution

where $\epsilon,\mu_{x},\mu_{y}$ are constants. We vary the value of $Q_{x}$ to investigate the dependence of the numerical diffusion on the Péclet number defined by Pe = $\dfrac{L\sqrt{Q_{x}^{2}+Q_{y}^{2}}}{D}$, with $L$ denoting the length of the computational domain. For this we compute at each time step the first moments

and the second moments

As results from above, the half-slopes of the centered second moments give the diffusion coefficients:

The difference between the numerically computed half-slopes and $D$ furnishes the numerical diffusion. The integrals are always aproximated like (this holds also for (T2))

with $f_{T}$ denoting the value of $f$ in the center of mass of $T$ and $|T|$ stands for the area of $T$.

We performed simulations with $D=0.01$, $\mu_{x}=\mu_{y}=5$, $Q_{x}=0,0.1,1,2$, $Q_{y}=0$ and $\epsilon=0.1$. Accordingly, we have global Péclet numbers of 0, 100, 1000 and 2000. To see the effect of the local Péclet number (Pe_l = $h\dfrac{\sqrt{Q_{x}^{2}+Q_{y}^{2}}}{D}$) we also performed simulations on two meshes: a coarse one, with a mesh diameter $h=0.1$ and a fine one with the mesh diameter $h=0.05$. The time step is taken $\tau=0.05$ and $\tau=0.025$ for the computations on the coarse and fine mesh, respectively. The local Péclet number is for the first case 0, 1, 10, and 20, whereas for the second one 0, 0.5, 5, and 10. The results are presented in Fig. 2 – 7. The difference of $D_{x}$ and $D_{y}$ to $D=0.01$ divided by $D=0.01$ gives the relative numerical diffusion. We did computations with the classical, new and upwind MHFEM schemes (see Sec. 2.2), linear, quadratic and upwind GFEM schemes (see Sec. 2.3) and first, second and third order FVM schemes (see Sec. 2.4). The numerical diffusion is very small for low Péclet numbers and obviously increases for high Péclet numbers. It can also be seen that the local Péclet number is more important than the global one, the results on a finer grid showing a clear improvement of the numerical diffusion. Nevertheless, the local Péclet number alone can not be used as a criterion for how diffusive the scheme is: the results with $Q_{x}=1$ on the coarse grid are much worse than the results for $Q_{x}=2$ on the fine grid, although the local Péclet numbers are identical.

The upwind schemes (MHFEM 2, partial or full upwind GFEM) should be applied only for high local Péclet numbers. These methods are stabilizing by introducing artificial diffusion [19] and this should be, if possible, avoided. We mention here that even for $Q_{x}=0.1$ the full upwind schemes show artificial diffusion (like $10-20\%$ but because of the scale of the plot this can not be seen in the pictures Fig. 2 or Fig. 4). The partial upwind GFEM scheme is constructed in such a way that no artificial diffusion is introduced if the local Péclet number is smaller or equal to one. Nevertheless, without applying any stabilization technique the solution becomes (uncontrolled) unphysical, at extreme the solution even blows up. The results for the coarse grid and high local Péclet number are clearly showing anomalous high numerical diffusion for all discretization schemes under consideration. This comportment is due to an oscillating solution, i.e. the solution shows extremely high and/or negative values. In this case the upwind schemes furnish much better results. The results obtained with FVM are better than MHFEM or GFEM. We also mention that the upwind MHFEM scheme (MHFEM 2) is less diffusive even than the partial upwind GFEM scheme. Further, the new MHFEM scheme (MHFEM 1) is better than the classical scheme (MHFEM 0). On the finer grid, for $Q_{x}=1$ the schemes without upwind are again better, i.e. the mesh is now fine enough. But for $Q_{x}=2$ we still have an oscillating solution even on the fine mesh and the upwind schemes are better. If we would refine once more we will again have the schemes without upwind in advantage! Consequently, when dealing with a convection-dominated-problem we have two possibilities: we decrease the size of the discretization parameters (time step and mesh size) till the solution shows no oscilliations or apply some stabilization techniques. The first method can be very time costly, the second one may have too much numerical diffusion. The art is normally to find a compromise of the two!

The higher order schemes in space (quadratic GFEM or second or third order FVM) are showing a bit less numerical diffusion than the lower order ones on the same grid (but the number of unknowns is much higher, hence higher computational time). This is more obvious in time: BDF(2) is much less numerically diffusive than BE (see Fig. 8 and Fig. 9). We also point out that the differences are much higher in the flow direction (for $D_{xx}$) than in the transverse one (for $D_{yy}$). But the superiority of the higher order schemes is clearly seen only for ’smooth’ problems like the one in (T1), for real application problems as considered in (T2), Sec. 3.2 the differences between the higher and lower order schemes are negligible (but not the differences in the computational time!).

In this section we study how diffusive the numerical schemes under consideration are for a realistic simulation of solute transport in heterogeneous soil. We consider here a saturated soil, i. e. $\Theta(\psi)=\Theta_{S}$ and $K(\Theta)=K_{S}$ are constants. We set $\Theta_{S}=1$ and we generate a statistically normally distributed log-hydraulic conductivity as explained below. The water flux is obtained by solving (2)–(3). Let $\Omega=[0,210]\times[0,85]$ be the computational domain and $T=200$ the final time. The diffusion-dispersion coefficient is given again by ${\bf D}=D\left(\begin{array}[]{cc}1&0\\ 0&1\end{array}\right)$, with $D=0.01$. The initial and the boundary conditions for both pressure and concentration are given in Fig. 10. The total solute mass is one. The physical units in this section are milligrams, meters and days and are not written every time explicitly.