A robust upwind mixed hybrid finite element method for transport in variably saturated porous media

The mixed finite element (MFE) method is well adapted for the simulation of fluid flow in heterogeneous porous media. However, when employed for the transport equation, it can generate solutions with strong unphysical oscillations because of the hyperbolic nature of advection. In this work, a robust upwind MFE scheme is proposed to avoid such unphysical oscillations. The new scheme is a combination of the upwind edge/face centered finite volume method with the hybrid formulation of the MFE method. The scheme ensures continuity of both advective and dispersive fluxes between adjacent elements and allows to maintain the time derivative continuous, which permits employment of high-order time integration methods via the method of lines (MOL).

Numerical simulations are performed in both saturated and unsaturated porous media to investigate the robustness of the new upwind MFE scheme. Results show that, contrarily to the standard scheme, the upwind MFE method generates stable solutions without under and overshoots. The simulation of contaminant transport into a variably saturated porous medium highlights the robustness of the proposed upwind scheme when combined with the MOL for solving nonlinear problems.

1 Introduction

The mixed finite element (MFE) method (Raviart and Thomas, 1977; Brezzi et al., 1985; Chavent and Jaffré, 1986; Brezzi and Fortin, 1991; Younes et al., 2010) is known to be a robust numerical scheme for solving elliptic diffusion problems such as the fluid flow in heterogeneous porous media. The method combines advantages of the finite volumes, by ensuring local mass conservation and continuity of fluxes between adjacent cells, and advantages of finite elements by easily handling heterogeneous domains with discontinuous parameter distributions and unstructured meshes. As a consequence, the MFE method has been largely used for flow in porous media (see, for instance, the review of Younes et al. (2010) and references therein). The hybridization technique has been largely used with the MFE method to improve its efficiency (Chavent and Roberts, 1991; Traverso et al., 2013). This technique allows to reduce the total number of unknowns and produces a final system with a symmetric positive definite matrix. The unknowns with the hybrid MFE method are the Lagrange multipliers which correspond to the traces of the scalar variable at edges/faces (Chavent and Jaffré, 1986).

When applied to transient diffusion equations with small time steps, the hybrid MFE method can produce solutions with small unphysical over- and undershoots (Hoteit et al., 2002a, b; Mazzia, 2008). A lumped formulation of the hybrid MFE method was developed by Younes et al. (2006) to improve its monotonicity and reduce nonphysical oscillations. The lumped formulation ensures that the maximum principle is respected for parabolic diffusion equations on acute triangulations (Younes et al., 2006). For more general 2D and 3D element shapes, the lumping procedure allows to significantly improve the monotonous character of the hybrid MFE solution (Younes et al., 2006; Koohbor et al., 2020). As an illustration, the lumped formulation was shown to be more efficient and more robust than the standard hybrid formulation for the simulation of the challenging nonlinear problem of water infiltration into an initially dry soil (Belfort et al., 2009). The lumped formulation has recently been used for flow discretization in the case of density-driven flow in saturated–unsaturated porous media (Younes et al., 2022a).

However, the MFE method remains little used for the discretization of the full transport equation. When employed to the advection–dispersion equation, the MFE method can generate solutions with strong numerical instabilities in the case of advection-dominated transport because of the hyperbolic nature of the advection operator. To avoid these instabilities, one of the most popular and easiest ways is to use an upwind scheme. Indeed, although upwind schemes introduce some numerical diffusion leading to an artificial smearing of the numerical solution, they avoid unphysical oscillations and remain useful, especially for large domains and regional field simulations. In the literature, some upwind mixed finite element schemes have been employed to improve the robustness of the MFE method for advection-dominated problems (Dawson, 1998; Dawson and Aizinger, 1999; Radu et al., 2011; Vohralík, 2007; Brunner et al., 2014).

The main idea of an upwind scheme for an element E is to calculate the mass flux exchanged with its adjacent element E^′ using the concentration from E in the case of an outflow and the concentration from E^′ in the case of an inflow. However, this idea cannot be applied as such with the hybrid MFE method since the hybridization procedure requires to express the flux at the element interface as only a function of variables at the element E (and not E^′). To overcome this difficulty, Radu et al. (2011) and Brunner et al. (2014) proposed an upwind MFE method where, in the case of an inflow, the concentration at the adjacent element E^′ is replaced by an approximation using the concentration at E and the trace of concentration at the interface ${\partial }_{E{E}^{\prime }}$ by assuming that the edge concentration is the mean of the concentrations in E and E^′. However, this assumption cannot be verified for a general configuration. Furthermore, with such an assumption, each of the advective and dispersive fluxes is discontinuous at the element interfaces, and continuity is only fulfilled for the total flux.

In this work, a new upwind MFE method is proposed for solving the full transport equation without requiring any approximation of the upwind concentration. The new scheme is a combination of the upwind edge/face centered finite volume (FV) scheme with the lumped formulation of the MFE method. It guarantees continuity of both advective and dispersive fluxes at element interfaces. Further, the new upwind MFE scheme maintains the time derivative continuous, and thus allows to employ high-order time integration methods via the method of lines (MOL), which was shown to be very efficient for solving nonlinear problems (see, for instance, Fahs et al., 2009 and Younes et al., 2009).

This article is structured as follows. In Sect. 2, we recall the hybrid MFE method for the discretization of the transport equation. In Sect. 3, we introduce the new upwind MFE method based on the combination of the upwind edge/face FV scheme with the lumped formulation of the MFE method. In Sect. 4, numerical experiments are performed for transport in saturated and unsaturated porous media to investigate the robustness of the new developed upwind MFE scheme. Some conclusions are given in the last section of the article.

2 The hybrid MFE method for the advection–dispersion equation

The mass conservation of the contaminant in variably saturated porous media is

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \left(\mathit{\theta }C\right)}{\partial t}}+\mathrm{\nabla }\cdot \left({\stackrel{\mathrm{̃}}{q}}_{\mathrm{a}}+{\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}\right)=\mathrm{0},\end{array}$$

where C is the normalized concentration [–], θ is water content [L³ L⁻³], t is time [T], ${\stackrel{\mathrm{̃}}{q}}_{\mathrm{a}}=qC$ is the advective flux with q the Darcy velocity [L T⁻¹], and ${\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}$ the dispersive flux given by

$$\begin{array}{}\text{(2)}& {\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}=-D\mathrm{\nabla }C,\end{array}$$

with D, the dispersion tensor, expressed by

$$\begin{array}{}\text{(3)}& D=D_{\mathrm{m}}I+\left({\mathit{\alpha }}_{\mathrm{L}}-{\mathit{\alpha }}_{\mathrm{T}}\right)q\otimes q/\left|q\right|+{\mathit{\alpha }}_{\mathrm{T}}\left|q\right|I,\end{array}$$

in which α_L and α_T are the longitudinal and transverse dispersivities [L], D_m is the pore water diffusion coefficient [L² T⁻¹], and I is the unit tensor.

Figure 1 Vectorial basis functions for the MFE method.

The water content θ and the Darcy velocity q are linked by the fluid mass conservation equation in variably saturated porous media:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}+\mathrm{\nabla }\cdot q=\mathrm{0}.\end{array}$$

Substituting Eq. (4) into Eq. (1) yields the following advection–dispersion equation:

$$\begin{array}{}\text{(5)}& \mathit{\theta }{\displaystyle \frac{\partial C}{\partial t}}+\mathrm{\nabla }\cdot \left({\stackrel{\mathrm{̃}}{q}}_{\mathrm{a}}+{\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}\right)-C\mathrm{\nabla }\cdot q=\mathrm{0}.\end{array}$$

In this work, we consider that the velocity q is obtained by solving Richards' equation using the hybrid MFE method. For a two-dimensional domain with a triangular mesh, q is approximated inside each triangle E using the lowest-order Raviart–Thomas (RT0) vectorial basis functions $w_j^E$:

$$\begin{array}{}\text{(6)}& q=\sum _{j=\mathrm{1}}^{\mathrm{3}}{Q}_j^E{w}_j^E,\end{array}$$

where $Q_j^E$ is the water flux across the edge E_j of E (see Fig. 1) and $w_j^E=\frac{\mathrm{1}}{\mathrm{2}\left|E\right|}\left(\begin{array}{c}x-x_j^E\\ y-y_j^E\end{array}\right)$ is the typical RT0 basis functions (Younes et al., 1999) with $\left(x_j^E,y_j^E\right)$ the coordinates of the node j opposite to the edge E_j of E, and $\left|E\right|$ the area of E.

To apply the hybrid MFE method to the transport Eq. (5), we approximate the dispersive flux ${\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}$ with RT0 vectorial basis functions as

$$\begin{array}{}\text{(7)}& {\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}=\sum _{j=\mathrm{1}}^{\mathrm{3}}{\stackrel{\mathrm{̃}}{Q}}_j^{\mathrm{d},E}{w}_j^E,\end{array}$$

where ${\stackrel{\mathrm{̃}}{Q}}_j^{\mathrm{d},E}=\underset{E_j}{\int }{\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}\cdot {\mathit{\eta }}_j^E$ is the dispersive flux across the edge E_j of the element E and ${\mathit{\eta }}_j^E$ is the outward unit normal vector to the edge E_j.

The variational formulation of Eq. (2) using the test function $w_i^E$ yields

$$\begin{array}{}\text{(8)}& \underset{E}{\int }{D}^{-\mathrm{1}}{\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}{w}_i^E=\underset{E}{\int }C\mathrm{\nabla }\cdot w_i^E-\sum _j\underset{E_j}{\int }C{w}_i^E\cdot {\mathit{\eta }}_j^E.\end{array}$$

Substituting Eq. (7) into Eq. (8) and using properties of the basis functions $w_j^E$ give

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\sum _j{\stackrel{\mathrm{̃}}{Q}}_j^{\mathrm{d},E}\underset{E}{\int }\left(D_E^{-\mathrm{1}}{w}_j^E\right)& \cdot w_i^E={\displaystyle \frac{\mathrm{1}}{\left|E\right|}}\underset{E}{\int }C-{\displaystyle \frac{\mathrm{1}}{\left|E_i\right|}}\underset{E_i}{\int }C\\ & =C_E-{\mathrm{TC}}_i^E,\end{array}\end{array}$$

in which D_E is the local dispersion tensor at the element E, C_E is the mean concentration at E, and TC${}_i^E$ is the edge (trace) concentration (Lagrange multiplier) at the edge E_i.

Denoting the local matrix ${\stackrel{\mathrm{̃}}{\mathbf{B}}}_{i,j}^E=\underset{E}{\int }\left(D_E^{-\mathrm{1}}{w}_j^E\right)\cdot w_i^E$, the inversion of the system of Eq. (9) gives the expression for the dispersive flux ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}$:

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}=\sum _j{\stackrel{\mathrm{̃}}{\mathbf{B}}}_{i,j}^{E,-\mathrm{1}}\left(C_E-{\mathrm{TC}}_j^E\right).\end{array}$$

Besides, the integration of the mass conservation Eq. (6) over the element E writes

$$\begin{array}{}\text{(11)}& \underset{E}{\int }\mathit{\theta }{\displaystyle \frac{\partial C}{\partial t}}+\underset{E}{\int }\mathrm{\nabla }\cdot {\stackrel{\mathrm{̃}}{q}}_{\mathrm{a}}+\underset{E}{\int }\mathrm{\nabla }\cdot {\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}-\underset{E}{\int }C\mathrm{\nabla }\cdot q=\mathrm{0},\end{array}$$

which becomes, using Green's formula,

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathit{\theta }}_E\left|E\right|{\displaystyle \frac{\partial C_E}{\partial t}}& +\sum _i\underset{E_i}{\int }Cq\cdot {\mathit{\eta }}_i^E+\sum _i\underset{E_i}{\int }{\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}\cdot {\mathit{\eta }}_i^E\\ & -\underset{E}{\int }C\mathrm{\nabla }\cdot q=\mathrm{0},\end{array}\end{array}$$

where θ_E is the water content of the element E.

Substituting Eq. (2) into Eq. (12) yields

$$\begin{array}{}\text{(13)}& {\mathit{\theta }}_E\left|E\right|{\displaystyle \frac{\partial C_E}{\partial t}}+\sum _i\underset{{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{t},E}}{\underbrace{\left({\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}+{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}\right)}}-C_E\sum _i{Q}_i^E=\mathrm{0},\end{array}$$

in which ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{t},E}={\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}+{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}$ is the total flux at the edge E_i with ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}$ the advective flux given by ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}=Q_i^E{\mathrm{TC}}_i^E$ and ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}$ the dispersive flux given by Eq. (10).

The hybridization of the MFE method is performed in the following two steps.

The flux Eq. (10) is substituted into the mass conservation Eq. (13), which is then discretized in time using the first-order implicit Euler scheme,

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\mathit{\theta }}_E{\displaystyle \frac{\left|E\right|}{\mathrm{\Delta }t}}\left(C_E^{n+\mathrm{1}}-C_E^n\right)& +\sum _i{Q}_i^E{\mathrm{TC}}_i^{E,n+\mathrm{1}}-C_E^{n+\mathrm{1}}\sum _i{Q}_i^E\\ & +{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}^E{C}_E^{n+\mathrm{1}}-\sum _i{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E{\mathrm{TC}}_i^{E,n+\mathrm{1}}=\mathrm{0},\end{array}\end{array}$$

in which ${\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E=\sum _j{\stackrel{\mathrm{̃}}{B}}_{i,j}^{E,-\mathrm{1}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\alpha }}}^E=\sum _i{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E$.

Hence, the mean concentration at the new time level $C_E^{n+\mathrm{1}}$ can be expressed as a function of ${\mathrm{TC}}_i^{E,n+\mathrm{1}}$, the concentration at the edges of E, as follows:

$$\begin{array}{}\text{(15)}& C_E^{n+\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\beta }}_E}}\sum _i\left({\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E-Q_i^E\right){\mathrm{TC}}_i^{E,n+\mathrm{1}}+{\displaystyle \frac{{\mathit{\lambda }}_E}{{\mathit{\beta }}_E}}{C}_E^n,\end{array}$$

in which ${\mathit{\lambda }}_E={\mathit{\theta }}_E\frac{\left|E\right|}{\mathrm{\Delta }t}$ and ${\mathit{\beta }}_E=\left({\mathit{\lambda }}_E+{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}^E-\sum _i{Q}_i^E\right)$.

The mean concentration given by Eq. (15) is then substituted into the flux Eq. (10), which allows expressing the dispersive flux ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E,n+\mathrm{1}}$ (the subscript n+1 will be omitted to alleviate the notations) as only a function of the traces of concentration at edges ${\mathrm{TC}}_i^{E,n+\mathrm{1}}$:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}& =\sum _j\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E}{{\mathit{\beta }}_E}}\left({\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_j^E-Q_j^E\right)-{\stackrel{\mathrm{̃}}{B}}_{i,j}^{E,-\mathrm{1}}\right){\mathrm{TC}}_j^{E,n+\mathrm{1}}\\ & +{\displaystyle \frac{{\mathit{\lambda }}_E}{{\mathit{\beta }}_E}}{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E{C}_E^n.\end{array}\end{array}$$

The system to be solved is obtained by imposing the continuity of the total flux $\left({\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{t},E}+{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{t},E^{\prime }}=\mathrm{0}\right)$ as well as the continuity of the trace of concentration $\left({\mathrm{TC}}_i^{E,n+\mathrm{1}}={\mathrm{TC}}_i^{E^{\prime },n+\mathrm{1}}\right)$ at the edge i between the two elements E and E^′ (Fig. 2).

Figure 2 Continuity of concentration and total flux between adjacent elements with the hybrid MFE method.

Note that the advective flux ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}$ is continuous between E and E^′ because of the continuity of the water flux and the continuity of the trace of concentration at the interface. Thus, for the continuity of the total flux $\left({\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{t},E}+{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{t},E^{\prime }}=\mathrm{0}\right)$, it is required that the dispersive flux is continuous:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{t},E}& +{\stackrel{\mathrm{̃}}{Q}}_i^{t,E^{\prime }}=\left(Q_i^E+Q_i^{E^{\prime }}\right){\mathrm{TC}}_i^{E,n+\mathrm{1}}+{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}\\ & +{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E^{\prime }}={\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}+{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E^{\prime }}=\mathrm{0}.\end{array}\end{array}$$

Using Eq. (16), we obtain

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\sum _j\left({\stackrel{\mathrm{̃}}{B}}_{i,j}^{E,-\mathrm{1}}\right.& \left.-{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E}{{\mathit{\beta }}_E}}\left({\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_j^E-Q_j^E\right)\right){\mathrm{TC}}_j^{E,n+\mathrm{1}}\\ & +\sum _j\left({\stackrel{\mathrm{̃}}{B}}_{i,j}^{E^{\prime },-\mathrm{1}}-{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^{E^{\prime }}}{{\mathit{\beta }}_{E^{\prime }}}}\left({\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_j^{E^{\prime }}-Q_j^{E^{\prime }}\right)\right)\\ & {\mathrm{TC}}_j^{E^{\prime },n+\mathrm{1}}={\displaystyle \frac{{\mathit{\lambda }}_E}{{\mathit{\beta }}_E}}{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E{C}_E^n+{\displaystyle \frac{{\mathit{\lambda }}_{E^{\prime }}}{{\mathit{\beta }}_{E^{\prime }}}}{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^{E^{\prime }}{C}_{E^{\prime }}^n.\end{array}\end{array}$$

The continuity Eq. (18) is written for all mesh edges, and the resulting equations form the final system to be solved for the traces of concentration at edges ${\mathrm{TC}}_i^{E,n+\mathrm{1}}$ as unknowns.

Note that the hybrid MFE Eq. (18), obtained by approximating the dispersive flux with RT0 basis functions, is equivalent to the new MFE method proposed in Radu et al. (2011).

3 The upwind and lumped MFE approaches

In this section, we recall the main principles of two existing approaches, developed to improve the stability of the MFE solution of the transport equation. The first approach is the upwind hybrid MFE scheme of Radu et al. (2011), developed for advection dominated transport. The second approach is the lumped hybrid MFE method of Younes et al. (2006), developed for dispersive transport.

3.1 The upwind hybrid MFE of Radu et al. (2011)

In the case of advection-dominated transport, solving the hybrid MFE Eq. (18) can yield solutions with strong instabilities. A common way to avoid such instabilities is to use an upwind scheme for the advective flux. Thus, for an element E, the advective flux ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}=Q_i^E{\mathrm{TC}}_i^E$ at the edge i (common with the element E^′) has to be calculated using either the concentration from E (if $Q_i^E>\mathrm{0}$) or the concentration from E^′ (if $Q_i^E<\mathrm{0}$). Radu et al. (2011) suggested replacing the advective flux ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}=Q_i^E{\mathrm{TC}}_i^E$ at the interface by

$$\begin{array}{}\text{(19)}& {\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}=\left\{\begin{array}{lc}{Q}_i^E{C}^E& \mathrm{if}{Q}_i^E>\mathrm{0}\\ Q_i^E{C}^{E^{\prime }}& \mathrm{if}{Q}_i^E<\mathrm{0}\end{array}\right.\end{array}$$

The advective term is now calculated using the upwind mean concentration, which can be that of the element E or of its adjacent element E^′.

The advective flux of Eq. (19) is rewritten in the following condensed form:

$$\begin{array}{}\text{(20)}& {\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}=Q_i^E\left({\mathit{\tau }}_i^E{C}^E+\left(\mathrm{1}-{\mathit{\tau }}_i^E\right)C^{E^{\prime }}\right),\end{array}$$

with ${\mathit{\tau }}_i^E=\mathrm{1}$ for an outflow $\left(Q_i^E>\mathrm{0}\right)$ and ${\mathit{\tau }}_i^E=\mathrm{0}$ for an inflow $\left(Q_i^E<\mathrm{0}\right)$.

However, this expression is incompatible with the hybridization procedure. Indeed, if we replace, in the Eq. (14), the advective term $Q_i^E{\mathrm{TC}}_i^E$ by Eq. (20), the latter will contain both C^E and $C^{E^{\prime }}$. Thus, the first step of the hybridization procedure cannot allow expressing $C_E^{n+\mathrm{1}}$ as only a function of ${\mathrm{TC}}_i^{E,n+\mathrm{1}}$ as in the Eq. (15).

To avoid this difficulty, Radu et al. (2011) suggested replacing $C^{E^{\prime }}$ by the following expression:

$$\begin{array}{}\text{(21)}& C^{E^{\prime }}\simeq \mathrm{2}{\mathrm{TC}}_i^E-C^E.\end{array}$$

This approximation is based on the assumption that ${\mathrm{TC}}_i^E\simeq \left(C^E+C^{E^{\prime }}\right)/\mathrm{2}$.

Plugging Eq. (21) into Eq. (20), the advective flux ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}$ depends only on the variables of the element E (mean concentration C^E and edge concentration ${\mathrm{TC}}_i^E$):

$$\begin{array}{}\text{(22)}& {\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}=Q_i^E\left({\mathit{\tau }}_i^E{C}^E-\left(\mathrm{1}-{\mathit{\tau }}_i^E\right)C^E+\mathrm{2}\left(\mathrm{1}-{\mathit{\tau }}_i^E\right){\mathrm{TC}}_i^E\right).\end{array}$$

Equation (22) can then be used to replace the advective term $Q_i^E{\mathrm{TC}}_i^{E,n+\mathrm{1}}$ in Eq. (14), and thus the hybridization procedure allows to express $C_E^{n+\mathrm{1}}$ as a function of ${\mathrm{TC}}_i^{E,n+\mathrm{1}}$ as in the Eq. (15). Then, the expression of $C_E^{n+\mathrm{1}}$ is substituted into the dispersive flux Eq. (10), and the final system is obtained by prescribing continuity of the total flux $\left({\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{t},E}+{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{t},E^{\prime }}=\mathrm{0}\right)$ at the interface between E and E^′. This scheme was shown to be more efficient (by using a sparser system matrix with fewer unknowns) than the non-hybrid upwind mixed method of Dawson (1998). The two methods yielded optimal first-order convergence in time and space (Brunner et al., 2014).

The assumption given by Eq. (21) can be a rough approximation, especially in the case of a heterogeneous domain where dispersion can vary with several orders of magnitudes from element to element. For such a situation, the edge concentration can be significantly different from the average of the mean concentrations of adjacent elements. Furthermore, the advective flux is not uniquely defined at the interface and can be different for the two adjacent elements E and E^′. For instance, in the case of $Q_i^E=Q>\mathrm{0}$, the advective flux leaving the element E is ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E}=Q{C}^E$, whereas the flux entering the element E^′ is ${\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{a},E^{\prime }}=Q\left(\mathrm{2}{\mathrm{TC}}_i^E-C^{E^{\prime }}\right)$ which could be different as ${\mathrm{TC}}_i^E$ is not necessarily the mean of C^E and $C^{E^{\prime }}$. In this situation, because of the discontinuity of the advective flux, the dispersive flux will not be continuous at the interface since the continuity is prescribed only for the total flux.

3.2 The lumped hybrid MFE scheme for dispersion transport

In this section, we recall the main principles of the lumped hybrid MFE method of Younes et al. (2006), developed to improve the stability of the MFE solution in the case of dispersive transport.

Considering only dispersion, Eq. (5) simplifies to

$$\begin{array}{}\text{(23)}& \mathit{\theta }{\displaystyle \frac{\partial C}{\partial t}}+\mathrm{\nabla }\cdot {\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}=\mathrm{0}.\end{array}$$

As detailed above, the hybrid MFE method for Eq. (23) is based on two stages:

Stage 1. Discretization of the transient mass conservation equation over the element E: the integration of the mass conservation Eq. (23) over the element E gives (see Eq. 13) the following:

$$\begin{array}{}\text{(24)}& {\mathit{\theta }}_E\left|E\right|{\displaystyle \frac{\partial C_E}{\partial t}}+\sum _i{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}=\mathrm{0}.\end{array}$$

Stage 2. Imposing the continuity of the flux across the edge i sharing the two elements E and E^′:

$$\begin{array}{}\text{(25)}& {\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E}+{\stackrel{\mathrm{̃}}{Q}}_i^{\mathrm{d},E^{\prime }}=\mathrm{0}.\end{array}$$

Note that the continuity Eq. (25) can be interpreted as a steady-state mass conservation equation at the edge level. Hence, the hybrid MFE discretization uses the transient mass conservation equation at the element level, given by Eq. (24), and the steady-state mass conservation at the edge level, given by Eq. (25). With the lumped hybrid MFE method of Younes et al. (2006), the transient term is taken into account at the edge level. Hence, the lumped formulation uses a steady-state mass conservation equation at the element level and a transient mass conservation equation at the edge level. The two stages of the lumped hybrid MFE are as follows:

Stage 1. Discretization of the steady-state mass conservation equation over E: the steady-state transport over the element E writes

$$\begin{array}{}\text{(26)}& \sum _i{\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_i^{\mathrm{d},E}=\mathrm{0},\end{array}$$

where ${\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_i^{\mathrm{d},E}$ is the steady-state dispersive flux across the edge E_i.

Therefore, the mean concentration of Eq. (15) becomes

$$\begin{array}{}\text{(27)}& C_E=\sum _i{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E}{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}^E}}{\mathrm{TC}}_i^E,\end{array}$$

and using Eq. (16), the steady-state dispersive flux writes

$$\begin{array}{}\text{(28)}& {\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_i^{\mathrm{d},E}=\sum _j\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_j^E}{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}^E}}-{\stackrel{\mathrm{̃}}{B}}_{i,j}^{E,-\mathrm{1}}\right){\mathrm{TC}}_j^E.\end{array}$$

Stage 2. Discretization of the transient mass conservation equation over the lumping region R_i: the edge centered finite volume discretization of the transient transport Eq. (23) over the lumping region R_i (hatched area in Fig. 3), associated with the edge i, writes

$$\begin{array}{}\text{(29)}& \underset{R_i}{\int }\mathit{\theta }{\displaystyle \frac{\partial C}{\partial t}}+\underset{R_i}{\int }\mathrm{\nabla }\cdot {\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}=\mathrm{0},\end{array}$$

where the lumping regions R_i is formed by the two simplex regions $S_i^E$ and $S_i^{E^{\prime }}$, for an inner edge i sharing the two elements E and E^′, and by the sole simplex region $S_i^E$ for a boundary edge. The simplex region $S_i^E$ is defined by joining the center of E with the nodes j and k forming the edge i.

Figure 3 The lumping region R_i associated with the edge i, sharing the elements E and E^′, and formed by the two simplex regions $S_i^E$ and $S_i^{E^{\prime }}$.

Associating the edge concentration ${\mathrm{TC}}_i^E$ to R_i (see Fig. 3 for notations), Eq. 29) gives

$$\begin{array}{}\text{(30)}& \left\{{\displaystyle \frac{\left|E\right|}{\mathrm{3}}}{\mathit{\theta }}_E{\displaystyle \frac{\partial {\mathrm{TC}}_i^E}{\partial t}}+{\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_{ij}^{\mathrm{d},E}+{\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_{ik}^{\mathrm{d},E}\right\}+{}^{\prime }=\mathrm{0},\end{array}$$

in which ${\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_{ij}^{\mathrm{d},E}$ and ${\mathrm{TC}}_{ij}^E$ are respectively the dispersive flux and the concentration at the interior interface (ij)^E between the simplex regions $S_i^E$ and $S_j^E$. The shortcut { }^′ designates the same contribution as { }, but of the adjacent element E^′, in the case of Eq. (30), it corresponds to $\frac{\left|E^{\prime }\right|}{\mathrm{3}}{\mathit{\theta }}_{E^{\prime }}\frac{\partial {\mathrm{TC}}_i^E}{\partial t}+{\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_{ij}^{\mathrm{d},E^{\prime }}+{\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_{ik}^{\mathrm{d},E^{\prime }}$.

Besides, applying the steady-state dispersive transport Eq. (26) on the simplex region $S_i^E$ yields

$$\begin{array}{}\text{(31)}& {\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_{ij}^{\mathrm{d},E}+{\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_{ik}^{\mathrm{d},E}+{\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_i^{\mathrm{d},E}=\mathrm{0}.\end{array}$$

Finally, substituting Eqs. (28 ) and (31) into the transport Eq. (30) gives the final system to solve with the lumped hybrid MFE scheme:

$$\begin{array}{}\text{(32)}& \left\{{\displaystyle \frac{\left|E\right|}{\mathrm{3}}}{\mathit{\theta }}_E{\displaystyle \frac{\partial {\mathrm{TC}}_i^E}{\partial t}}+\sum _j\left({\stackrel{\mathrm{̃}}{B}}_{i,j}^{E,-\mathrm{1}}-{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_j^E}{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}^E}}\right){\mathrm{TC}}_j^E\right\}+{}^{\prime }=\mathrm{0}..\end{array}$$

Note that the lumped hybrid formulation Eq. (32) and the standard hybrid formulation (Eqs. 24 and 25) are exactly the same in the case of steady-state diffusion transport.

In the lumped formulation Eq. (32), the term of mass (with time derivative) has a contribution only on the diagonal term of the final system matrix. This improves the monotonous character of the solution (see Younes et al., 2006). For instance, in the case of an acute triangulation, the maximum principle is respected by the lumped formulation Eq. (32) whatever the heterogeneity of the porous medium (Younes et al., 2006).

Contrarily to the standard hybrid MFE scheme, where the discretization of the temporal derivative performed in Eq. (14) was necessary to obtain the final system given by Eq. (18), the lumped scheme given by Eq. (32) keeps the time derivative continuous which allows the use of efficient high-order temporal discretization methods via the MOL.

In the case of 2D triangular elements, the lumped formulation Eq. (32) is algebraically equivalent to the nonconforming Crouzeix–Raviart (Crouzeix and Raviart, 1973) finite element method (see Younes et al., 2008). The nonconforming Crouzeix–Raviart method uses the chapeau functions as basis functions to approximate the concentration, like the standard finite element method, but seed nodes are the midpoints of the edges.

4 The new upwind hybrid MFE scheme for advection–dispersion transport

To avoid the rough approximation (Eq. 21), we develop hereafter a new upwind MFE scheme where the advection term is calculated using upwind edge concentration instead of upwind mean concentration of the element E. The idea of the scheme is to extend the lumped hybrid MFE procedure to transport by both advection and dispersion and to use an upwind edge centered FV scheme to avoid unphysical oscillations caused by the hyperbolic nature of advection.

The integration of the whole mass conservation Eq. (5) over the lumping region R_i writes

$$\begin{array}{}\text{(33)}& \underset{R_i}{\int }\mathit{\theta }{\displaystyle \frac{\partial C}{\partial t}}+\underset{R_i}{\int }\mathrm{\nabla }\cdot \left(qC\right)+\underset{R_i}{\int }\mathrm{\nabla }\cdot {\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}-\underset{R_i}{\int }C\mathrm{\nabla }\cdot q=\mathrm{0}.\end{array}$$

Using notations of Fig. 3, we obtain

$$\begin{array}{}\text{(34)}& \begin{array}{rl}\left\{{\displaystyle \frac{\left|E\right|}{\mathrm{3}}}{\mathit{\theta }}_E{\displaystyle \frac{\partial {\mathrm{TC}}_i^E}{\partial t}}\right.& +Q_{ij}^E{\mathrm{TC}}_{ij}^E+Q_{ik}^E{\mathrm{TC}}_{ik}^E+{\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_{ij}^{\mathrm{d},E}+{\underset{\mathrm{¯}}{\stackrel{\mathrm{̃}}{Q}}}_{ik}^{\mathrm{d},E}\\ & \left.-{\mathrm{TC}}_i^E\left(Q_{ij}^E+Q_{ik}^E\right)\right\}+\mathit{\{}{\mathit{\}}}^{\prime }=\mathrm{0},\end{array}\end{array}$$

in which $Q_{ij}^E$ is the water flux at the interior interface (ij)^E, evaluated using the RT0 approximation of the velocity given by Eq. (6), which yields

$$\begin{array}{}\text{(35)}& Q_{ij}^E={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\left(Q_j^E-Q_i^E\right).\end{array}$$

Using Eqs. (28) and (31) and denoting ${\mathit{\lambda }}_E={\mathit{\theta }}_E\frac{\left|E\right|}{\mathrm{3}}$, Eq. (34) becomes

$$\begin{array}{}\text{(36)}& \begin{array}{rl}\left\{{\mathit{\lambda }}_E{\displaystyle \frac{\partial {\mathrm{TC}}_i^E}{\partial t}}\right.& +\sum _j\left({\stackrel{\mathrm{̃}}{B}}_{i,j}^{E,-\mathrm{1}}-{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_j^E}{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}^E}}\right){\mathrm{TC}}_j^E+Q_{ij}^E{\mathrm{TC}}_{ij}^E\\ & \left.+Q_{ik}^E{\mathrm{TC}}_{ik}^E-\left(Q_{ij}^E+Q_{ik}^E\right){\mathrm{TC}}_i^E\right\}+\mathit{\{}{\mathit{\}}}^{\prime }=\mathrm{0}.\end{array}\end{array}$$

The interior concentration ${\mathrm{TC}}_{ij}^E$ at the interface between the simplex regions $S_i^E$ and $S_j^E$ is calculated using an upwind scheme (see Fig. 3) defined by

$$\begin{array}{}\text{(37)}& {\mathrm{TC}}_{ij}^E={\mathit{\tau }}_{ij}^E{\mathrm{TC}}_i^E+\left(\mathrm{1}-{\mathit{\tau }}_{ij}^E\right){\mathrm{TC}}_j^E,\end{array}$$

with ${\mathit{\tau }}_{ij}^E=\mathrm{1}$ if $\left(Q_{ij}^E\ge \mathrm{0}\right)$, else ${\mathit{\tau }}_{ij}^E=\mathrm{0}$.

Thus, the final system to solve becomes

$$\begin{array}{}\text{(38)}& \begin{array}{rl}\left\{{\mathit{\lambda }}_E{\displaystyle \frac{\partial {\mathrm{TC}}_i^E}{\partial t}}\right.& +\sum _j\left({\stackrel{\mathrm{̃}}{B}}_{i,j}^{E,-\mathrm{1}}-{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_j^E}{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}^E}}\right){\mathrm{TC}}_j^E\\ & +Q_{ij}^E\left(\mathrm{1}-{\mathit{\tau }}_{ij}^E\right)\left({\mathrm{TC}}_j^E-{\mathrm{TC}}_i^E\right)\\ & \left.+Q_{ik}^E\left(\mathrm{1}-{\mathit{\tau }}_{ik}^E\right)\left({\mathrm{TC}}_k^E-{\mathrm{TC}}_i^E\right)\right\}+\mathit{\{}{\mathit{\}}}^{\prime }=\mathrm{0}.\end{array}\end{array}$$

Figure 4 Description of the problem of the contamination of a 2D saturated porous medium.

In the case of a first-order Euler implicit time discretization, Eq. (37) becomes

$$\begin{array}{}\text{(39)}& \begin{array}{rl}& \left\{\begin{array}{l}\sum _j\left({\stackrel{\mathrm{̃}}{B}}_{i,j}^{E,-\mathrm{1}}-\frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_i^E{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_j^E}{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}^E}\right){\mathrm{TC}}_j^{E,n+\mathrm{1}}+{\mathit{\lambda }}_E{\mathrm{TC}}_i^{E,n+\mathrm{1}}\\ +Q_{ij}^E\left(\mathrm{1}-{\mathit{\tau }}_{ij}^E\right)\left({\mathrm{TC}}_j^{E,n+\mathrm{1}}-{\mathrm{TC}}_i^{E,n+\mathrm{1}}\right)\\ +Q_{ik}^E\left(\mathrm{1}-{\mathit{\tau }}_{ik}^E\right)\left({\mathrm{TC}}_k^{E,n+\mathrm{1}}-{\mathrm{TC}}_i^{E,n+\mathrm{1}}\right)-{\mathit{\lambda }}_E{\mathrm{TC}}_i^{E,n}\end{array}\right\}\\ & +\mathit{\{}{\mathit{\}}}^{\prime }=\mathrm{0},\end{array}\end{array}$$

where ${\mathit{\lambda }}_E={\mathit{\theta }}_E\frac{\left|E\right|}{\mathrm{3}\mathrm{\Delta }t}$.

It is easy to see that, due to upwinding, the system matrix corresponding to Eq. (28) is always an M-matrix (a non-singular matrix with m_ii>0, m_ij≤0) in the case of transport by advection. The M-matrix property insures the stability of the scheme since it guaranties the respect of the discrete maximum principle i.e., local maxima or minima will not appear in the C solution in a domain without local sources or sinks.

Further, Eq. (37) expresses the total exchange between E and E^′ and therefore reflects the continuity of the total advection–dispersion flux between them. Both advective and dispersive fluxes are continuous between the adjacent elements E and E^′. The advective flux, calculated using the upwind edge concentration, is uniquely defined at the interface of the lumping region and is therefore continuous. As a consequence, the dispersive flux is also continuous between E and E^′ since the total flux is continuous at the interface between them.

5 Numerical experiments

In this section, a first test case dealing with transport in saturated porous media is simulated with the standard hybrid MFE and the new upwind MFE schemes. The results are compared against an analytical solution in order to validate the new developed scheme and to show its robustness for solving advection-dominated transport problems compared to the standard one. The second test case deals with transport in the unsaturated zone and aims to investigate the robustness of the new scheme when combined with the MOL for solving highly nonlinear problems.

5.1 Transport in saturated porous media: comparison against a 2D analytical solution

The hybrid and upwind MFE formulations are compared against the analytical solution developed by Leij and Dane (1990) for a simplified 2D transport problem (Fig. 4). The test case has been employed by Putti et al. (1990) and Siegel et al. (1997) for the verification of transport codes. It deals with the contamination from the left boundary of a 2D rectangular domain of dimension (0, 100) × (0, 40).

The boundary conditions for the transport are of Dirichlet type at the inflow (left vertical boundary), with

$$\begin{array}{}\text{(40)}& C=\left\{\begin{array}{cl}\mathrm{0}& \mathrm{for}x=\mathrm{0}\mathrm{and}\mathrm{0}\le y<\mathrm{12}\\ \mathrm{1}& \mathrm{for}x=\mathrm{0}\mathrm{and}\mathrm{12}\le y\le \mathrm{28}\\ \mathrm{0}& \mathrm{for}x=\mathrm{0}\mathrm{and}\mathrm{28}<y\le \mathrm{40}\end{array}\right..\end{array}$$

A zero diffusive flux is imposed at the right vertical outflow boundary. The top and bottom are impermeable boundaries. A uniform horizontal flow occurs from left to right with a constant flux q_x=0.5 m d⁻¹ prescribed at the left vertical boundary and a fixed head H=100 m at the right vertical boundary. The longitudinal and transverse dispersivities are α_L=0.2 m and α_T=0.05 m, respectively. The domain is discretized with a fine unstructured triangular mesh formed by 33 216 elements, and the simulation is performed for a final simulation time T=30 d using the Euler implicit time discretization with a fixed time step of 0.1 d. The linear systems are solved in each time step with a direct solver using an unsymmetric-pattern multi-frontal method and a direct sparse LU factorization (UMFPACK).

The analytical solution of this test case for an infinite domain is given by Leij and Dane (1990),

$$\begin{array}{}\text{(41)}& \begin{array}{rl}{C}_{\mathrm{analy}}(x,y,t)& ={\displaystyle \frac{x}{{\left(\mathrm{16}\mathit{\pi }{\mathit{\alpha }}_{\mathrm{L}}\right)}^{\mathrm{1}/\mathrm{2}}}}\underset{\mathrm{0}}{\overset{T}{\int }}{\mathit{\tau }}^{-\mathrm{3}/\mathrm{2}}\left\{\mathrm{erf}\left[{\displaystyle \frac{y-\mathrm{12}}{{\left(\mathrm{4}{\mathit{\alpha }}_{\mathrm{T}}\mathit{\tau }\right)}^{\mathrm{1}/\mathrm{2}}}}\right]\right.\\ & \left.+\mathrm{erf}\left[{\displaystyle \frac{\mathrm{28}-y}{{\left(\mathrm{4}{\mathit{\alpha }}_{\mathrm{T}}\mathit{\tau }\right)}^{\mathrm{1}/\mathrm{2}}}}\right]\right\}\exp \left[-{\displaystyle \frac{(x-\mathit{\tau }{)}^{\mathrm{2}}}{\mathrm{4}{\mathit{\alpha }}_{\mathrm{L}}\mathit{\tau }}}\right]\mathrm{d}\mathit{\tau },\end{array}\end{array}$$

with $\mathrm{erf}\left(x\right)=\frac{\mathrm{2}}{\sqrt{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{x}{\int }}\exp \left(-{\mathit{\tau }}^{\mathrm{2}}\right)\mathrm{d}\mathit{\tau }$.

The final distributions of the concentration with both hybrid MFE and upwind MFE schemes are depicted in Fig. 5. Although we have used an unstructured mesh, the two schemes yield almost symmetrical results. The hybrid MFE scheme (Fig. 5a) yields a solution with unphysical oscillations. Indeed, around 1.2 % of the contaminated region (i.e., the region with $\left|C\right|\ge {\mathrm{10}}^{-\mathrm{5}}$) exhibits unphysical oscillations with 0.4 % of the contaminated region with $C\le -{\mathrm{10}}^{-\mathrm{3}}$ and 0.8 % of the contaminated region with C≥1.001. These unphysical oscillations, although they seem moderate, can be dramatic, for instance, when dealing with reactive transport where some reactions occur only if the concentration exceeds a certain threshold. The solution obtained with the new upwind formulation (Fig. 5b) is monotone (all concentrations are between 0 and 1) which is in agreement with the physics. However, these results come at the expense of some numerical diffusion added to the solution. To appreciate the quality of both solutions and validate the upwind MFE method, we compare the concentration profile of the two methods to the analytical solution of Leij and Dane (1990) for a horizontal section located at y=20 m and a vertical section located at x=20 m.

Figure 5 Concentration distribution with the hybrid MFE and the upwind MFE methods for the 2D saturated transport problem (only the region x≤70 m is depicted).

The results of Fig. 6 show that the solution of both hybrid MFE and upwind MFE methods are in very good agreement with the analytical solution, which validates the new upwind MFE numerical model. Note, however, that a small numerical diffusion is observed with the upwind MFE solution, which is especially visible in Fig. 6b. Indeed, for the simulated problem, the transverse dispersivity is much smaller than the longitudinal one, and, as a consequence, the concentration front is sharper in the vertical section than in the horizontal one. This explains why the numerical diffusion generated by the upwind MFE method is more pronounced in Fig. 6b than in Fig. 6a.

Figure 6 Concentration profiles at y=20 m (a) and x=20 m (b) with the analytical, hybrid MFE, and upwind MFE solutions.

The test problem is then simulated using different mesh refinements to investigate the order of convergence of the new method. We start with a uniform mesh formed by 1000 triangles and a time step Δt=0.1 s. In each level of refinement, each triangle is subdivided into four similar triangles, by joining the three mid-edges and the time step Δt is halved. The following error is computed (Brunner et al., 2014):

$$\begin{array}{}\text{(42)}& \begin{array}{rl}\mathrm{Er}& =\left\{{∥C_{\mathrm{analyt}}\left(t^N\right)-C\left(t^N\right)∥}_{\mathrm{0}}^{\mathrm{2}}\right.\\ & {\left.+\mathrm{\Delta }t\sum _{n=\mathrm{1}}^N{∥{\stackrel{\mathrm{̃}}{q}}_{\mathrm{analyt}}^t\left(t^n\right)-{\stackrel{\mathrm{̃}}{q}}^t\left(t^n\right)∥}_{\mathrm{0}}^{\mathrm{2}}\right\}}^{\mathrm{1}/\mathrm{2}},\end{array}\end{array}$$

where ${\stackrel{\mathrm{̃}}{q}}^t={\stackrel{\mathrm{̃}}{q}}_{\mathrm{a}}+{\stackrel{\mathrm{̃}}{q}}_{\mathrm{d}}$ is the total advection–dispersion flux and N the total number of time steps.

The runs are performed on a single computer with an Intel Xeon E-2246G processor and 32 GB memory. The results of the computations, summarized in Table 1, clearly show optimal first-order convergence in space and time for the developed upwind hybrid MFE method.

Table 1 Numerical results for the new upwind hybrid MFE method.

5.2 Transport in a variably-saturated porous medium

In this test case, the developed upwind MFE method is combined with the MOL for solving contaminant transport in a variably-saturated porous medium. The advection–dispersion equation is transformed to an ordinary differential equation (ODE) using the new upwind MFE formulation for the spatial discretization, whereas the time derivative is maintained continuously. Therefore, high-order time integration methods included in efficient ODE solvers can be employed. With these solvers, both the time step size and the order of the time integration can vary during the simulation to deliver accurate results in an acceptable computational time.

To investigate the robustness and efficiency of the combination of the developed upwind MFE method with the MOL, we simulate in this section the problem of contaminant infiltration into a variably-saturated porous medium.

The domain (Fig. 7) is a rectangular box of 3 m × 2 m, filled with sand, with an initial water table at 0.65 m and hydrostatic pressure distribution. An infiltration of a tracer contaminant is applied over the left-most 0.1 m of the surface with a constant flux of 10⁻⁶ m s⁻¹. The right vertical side has a fixed head H=0.65 m below the water table and an impermeable boundary above it. The left vertical side and the upper (except the infiltration zone) and bottom boundaries are impermeable boundaries.

Figure 7 Description of the problem of contaminant infiltration into a 2D variably-saturated porous medium.

In this problem, the flow and transport are coupled by the velocity, which is obtained by solving the following pressure-head form of the nonlinear Richards' equation:

$$\begin{array}{}\text{(43)}& {\displaystyle }& {\displaystyle }\left(c\left(h\right)+S_{\mathrm{S}}{\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{S}}}}\right){\displaystyle \frac{\partial H}{\partial t}}+\mathrm{\nabla }\cdot q=\mathrm{0}\text{(44)}& {\displaystyle }& {\displaystyle }q=-k_{\mathrm{r}}K\mathrm{\nabla }H,\end{array}$$

with S_S the specific mass storativity related to head changes [L⁻¹], $H=h+y$ the equivalent head [L], $h=\frac{P}{\mathit{\rho }g}$ the pressure head, P the pressure [Pa], ρ the fluid density [M L⁻³], g the gravity acceleration [L T⁻²], y the upward vertical coordinate [L], c(h) the specific moisture capacity [L⁻¹], θ_S the saturated water content [L³ L⁻³], q the Darcy velocity [L T⁻¹], $K=\frac{\mathit{\rho }g}{\mathit{\mu }}k$ the hydraulic conductivity [L T⁻¹], k the permeability [L²], μ the fluid dynamic viscosity [M L⁻¹ T⁻¹], and k_r the relative conductivity [–].

We use the standard van Genuchten (1980) model for the relationship between water content and pressure head,

$$\begin{array}{}\text{(45)}& S_{\mathrm{e}}={\displaystyle \frac{\mathit{\theta }\left(h\right)-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}=\left\{\begin{array}{ll}\frac{\mathrm{1}}{{\left(\mathrm{1}+|\mathit{\alpha }h{|}^n\right)}^m}& h<\mathrm{0}\\ \mathrm{1}& h\ge \mathrm{0}\\ ,& \end{array}\right.\end{array}$$

where α [L⁻¹] and n [–] are the van Genuchten parameters, $m=\mathrm{1}-\mathrm{1}/n$, S_e [–] is the effective saturation, and θ_r [–] is the residual water content. The conductivity–saturation relationship is derived from the Mualem (1976) model,

$$\begin{array}{}\text{(46)}& k_{\mathrm{r}}=S_{\mathrm{e}}^{\mathrm{1}/\mathrm{2}}{\left[\mathrm{1}-{\left(\mathrm{1}-S_{\mathrm{e}}^{\mathrm{1}/m}\right)}^m\right]}^{\mathrm{2}}.\end{array}$$

The material properties of the test problem are given in Table 2.

Table 2 Parameters for the problem of infiltration into a 2D variably-saturated porous medium.

Figure 8 Concentration distribution, with the hybrid MFE (a) and the upwind MFE (b) schemes for the transport problem with high dispersion in a variably-saturated porous medium.

The simulation is performed for 80 h using a triangular mesh formed by 4273 triangular elements. Two test cases are investigated. In the first test case, the longitudinal and transverse dispersivities are α_L=0.03 m and α_T=0.003 m, respectively. The second test case is less diffusive with α_L=0.01 m and α_T=0.001 m.

The coupled nonlinear flow–transport system is solved using the MOL, which allows the use of efficient high-order time integration methods, for both the hybrid MFE and the upwind MFE schemes. To this aim, a hybrid MFE formulation with continuous time derivative was developed by extending the lumping procedure, developed in Younes et al. (2006) for the flow equation, to the advection–dispersion transport Eq. (5).

The time integration is performed with the DASPK time solver which uses an efficient automatic time-stepping scheme based on the fixed leading coefficient backward difference formulas (FLCBDF). The linear systems arising at each time step are solved with the preconditioned Krylov iterative method. The nonlinear problem is linearized using the Newton method with a numerical approximation of the Jacobian matrix.

The results of the hybrid MFE and the upwind MFE methods are depicted in Fig. 8 for the first test case involving high dispersion. Good agreement can be observed between the results of the hybrid MFE (Fig. 8a) and upwind MFE (Fig. 8b) schemes when combined with the MOL. In these figures, the contaminant progresses essentially vertically through the unsaturated zone of the soil. When the saturated zone is reached, the contaminant progresses horizontally and remains close to the water table. Note that the results of both schemes are stable and free from unphysical oscillations (Fig. 8a and b).

For the second test case with lower dispersion (α_L=0.01 m, α_T=0.001 m), the hybrid MFE method yields unstable results containing unphysical oscillations (red color in Fig. 9a). These oscillations hamper the convergence of the numerical model, and severe convergence issues can be encountered if we further decrease the dispersivity values. The results of the upwind MFE scheme are monotone and do not contain any unphysical oscillation (Fig. 9b). These results point out the robustness of the new upwind MFE method for transport in saturated and unsaturated porous media. The developed transport scheme has recently been successfully combined with the MFE method for fluid flow to simulate nonlinear flow and transport in unsaturated fractured porous media using the 1D–2D discrete fracture matrix (DFM) approach (Younes et al., 2022b).

Figure 9 Concentration distribution with the hybrid MFE (a) and upwind MFE (b) methods for the transport problem with low dispersion in variably-saturated porous medium.

6 Conclusion

MFE is a robust numerical method well adapted for diffusion problems on heterogeneous domains and unstructured meshes. When applied to transport equations, the MFE solution can exhibit strong unphysical oscillations due to the hyperbolic nature of advection. Upwind schemes can be used to avoid such oscillations, although they introduce some numerical diffusion. In this work, we developed an upwind scheme that does not require any approximation for the upwind concentration. The method can be seen as a combination of an upwind edge/face centered FV method with the lumped formulation of the hybrid MFE method. It ensures continuity of both advective and dispersive fluxes between adjacent elements and allows to maintain the time derivative continuous, which facilitates employment of high-order time integration methods via the method of lines (MOL) for nonlinear problems.

Numerical simulations for the transport in a saturated porous medium show that the standard hybrid MFE method can generate unphysical oscillations due to the hyperbolic nature of advection. These unphysical oscillations are completely avoided with the new upwind MFE scheme. The simulation of the problem of contaminant transport in a variably-saturated porous medium shows that only the upwind MFE scheme provides a stable solution. The results point out the robustness of the developed upwind MFE scheme when combined with the MOL for solving nonlinear transport problems.