A Generalized Finite Volume Method for Density Driven Flows in Porous Media

Abstract

In this article, we consider a time evolution equation for solute transport, coupled with a pressure equation in space dimension 2. For the numerical discretization, we combine the generalized finite volume method SUSHI on adaptive meshes with a time semi-implicit scheme. In the first part of this article, we present numerical simulations for two problems: a rotating interface between fresh and salt water and a well-known test case proposed by Henry. In the second part, we also introduce heat transfer and perform simulations for a system from the documentation of the software SEAWAT.

1. Introduction

In this article, we present numerical simulations for density driven flows involved in the production of lithium batteries.

This work was performed in the context of an exploratory CNRS project on the numerical simulation of variable density flows. More precisely, the objective of the project was related to the exploitation of lithium deposits in salt lakes, also known as “salars”. In recent years, lithium has become a strategic element for developed countries, as it is the basic element of lithium-ion batteries used in hybrid and electric vehicles. Therefore, its production is of great interest to all major groups involved in the automotive industry and to the suppliers of these groups.

From a mathematical point of view, we study a system of equations that describes the interaction between flow and transport in a porous medium in which the density $\rho$ increases strictly with respect to the concentration u of the transported species. Specifically, the equations governing density-dependent transport are Darcy’s law (1), the continuity equation for the fluid (2), and the continuity equation for the solute (3), which are given as

$$\{\begin{array}{}\mathrm{(1)}& \mathbf{V}=-{\displaystyle \frac{k}{\nu }}\left(\nabla P+\rho \left(u\right)g\nabla z\right)\mathrm{(2)}& \theta \frac{\partial \rho \left(u\right)}{\partial t}+\nabla ·\left(\mathbf{V}\rho \left(u\right)\right)={\rho }_s{Q}_s\mathrm{(3)}& \theta \frac{\partial \left(\rho \left(u\right)u\right)}{\partial t}+\nabla ·\left(\mathbf{V}\rho \left(u\right)u-\rho \left(u\right)K\nabla u\right)=u_s{\rho }_s{Q}_s\end{array}$$

in $D\times (0,T)$ together with suitable boundary conditions of the form

$$\left\{\begin{array}{cc}u=u_D(\mathbf{x},t)\hfill & \mathrm{on} D_D^u\times (0,T)\hfill \\ {\displaystyle \frac{\partial u}{\partial n}={\overline{u}}_N(\mathbf{x},t)}\hfill & \mathrm{on} D_N^u\times (0,T)\hfill \\ P=P_D(\mathbf{x},t)\hfill & \mathrm{on} D_D^p\times (0,T)\hfill \\ \mathbf{V}·\mathbf{n}={\overline{V}}_N(\mathbf{x},t)\hfill & \mathrm{on} D_N^p\times (0,T)\hfill \end{array}\right.$$

(4)

and the initial condition

$$u(\mathbf{x},0)=u_0\left(\mathbf{x}\right)$$

(5)

in D, with $\partial D=\overline{\partial D_D^u}\bigcup \overline{\partial D_N^u}=\overline{\partial D_D^P}\bigcup \overline{\partial D_N^P}$ where $\partial D_D^u$ and $\partial D_D^P$ correspond to Dirichlet boundary conditions and $\partial D_N^u$ and $\partial D_N^P$ correspond to Neumann boundary conditions. $\mathbf{V}$ is the flow velocity, P [Pa] is the pressure, and u [kg/m${}^3$] is the concentration of the species to be transported. The porosity $\theta$ is the ratio of voids to the total volume, D [m${}^2$/s] is the dispersion-diffusion tensor, k [m${}^2$] is the permeability, $\nu$ [kg/(m·s)] is the dynamic viscosity, $\rho$ [kg/m${}^3$] is the density, and $-g\nabla z$ [m/s${}^2$] is the gravity. $Q_s$ models the source or sink term of the fluid with density ${\rho }_s$, and $u_s$ is the normalized mass fraction of the source/sink term in reference to the maximum mass fraction. In addition, we apply a constitutive equation that expresses the relationship between the fluid density and the concentration of a dilute solution under isothermal conditions, and express it as

$$\rho \left(u\right)={\rho }_0(1+\overline{a}u)$$

(6)

where $\overline{a}={\rho }_{max}/{\rho }_0-1$ [1].

We refer to [2,3] for modeling aspects. As for numerical calculations, Hilhorst et al., 2012 [4] have applied the standard finite volume method to the system (1). However, in order to discretize the diffusion term by the standard finite volume method, a conforming mesh that satisfies the orthogonality condition must be used. The generalized finite volume method SUSHI [5] allows calculations on nonconforming meshes (see Definition 1 below). This is important in applying adaptive discretization methods on squares or cubes of different sizes. By using such an adaptive mesh, numerical calculations can be performed with a smaller number of volume elements. We propose a criterion for selecting the volume elements to be refined based on the discrete difference and we refer to the Section 2.2 for more details.

In order to apply the generalized finite volume method SUSHI, we set ${\displaystyle w:=w\left(u\right)=}$${\int }_0^u\rho \left(s\right)ds$, $m\left(w\right)=\rho \left(u\right)u$ and $\varrho \left(w\right)=\rho \left(u\right)$ and rewrite the system as

$$\left\{\begin{array}{c}\mathbf{V}=-{\displaystyle \frac{k}{\nu }}\left(\nabla P+\varrho \left(w\right)g(\nabla z)\right)\hfill \\ \theta {\displaystyle \frac{\partial \varrho \left(w\right)}{\partial t}}+\nabla ·\left(\mathbf{V}\varrho \left(w\right)\right)={\rho }_s{Q}_s\hfill \\ \theta \frac{\partial m\left(w\right)}{\partial t}+\nabla ·\left(\mathbf{V}m\left(w\right)\right)-\nabla ·(K\nabla w)=w_s{\rho }_s{Q}_s\hfill \end{array}\right.$$

(7)

in $D\times (0,T)$, where the unknown functions are now P and w. In the benchmarks, the initial condition and boundary condition of w are directly computed from the corresponding conditions on u. We consider two benchmarks, a problem involving a rotating interface and Henry’s problem. An advantage of finite volume methods is that the discrete integrals of the functions $\rho \left(u\right)$ and $\rho \left(u\right)u$ are numerically accurately preserved. We discuss the mass conservation in Section 3.4.

We then take heat transfer into consideration. Mathematically, we solve a system for the hydraulic head h, the solute concentration u, and the temperature $\mathsf{\Theta }$; note that the fluid density $\rho$ and the viscosity $\nu$ are given functions of h, u, and $\mathsf{\Theta }$.

The organization of this article is as follows. In Section 2, we introduce the time and space discretization and the generalized finite volume method SUSHI and present the method of mesh refinement. In Section 3, we present the numerical algorithm, which we apply. We present the results of numerical simulations for two problems: a problem for a rotating interface and a problem proposed by Henry. In Section 4, we consider the density driven flow problem coupled with heat transfer and apply a time semi-implicit scheme coupled to the SUSHI method for the space discretization. We perform simulations and compare our results with the results of SEAWAT. Our results are very close to the results obtained by SEAWAT.

2. Numerical Method and Adaptive Mesh

2.1. The Generalized Finite Volume Method SUSHI

We first introduce standard notations for the discretization in space and time.

Definition 1

(Discretization of D). Let D be a polyhedral open bounded connected subset of ${\mathbb{R}}^d$ and $\partial D=\overline{D}\setminus D$ its boundary. A discretization of D, denoted by $\mathcal{D}$, is defined as the triplet $\mathcal{D}=(\mathcal{T},\mathcal{E},\mathcal{P})$, where:

1. $\mathcal{T}$ is a finite family of non empty convex open disjoint subsets of D (the “control volumes”) such that $\overline{D}={\bigcup }_{p\in \mathcal{T}}\overline{p}$. For all $p\in \mathcal{T}$, we denote by $\partial p=\overline{p}\setminus p$ the boundary of p; and by $\left|p\right|$ its measure.

2. $\mathcal{E}$ is a finite family of disjoint subsets of $\overline{D}$ (the “interfaces”), such that, for all $\sigma \in \mathcal{E}$, σ is a nonempty open subset of a hyperplane of ${\mathbb{R}}^d$ and we denote its measure by $\left|\sigma \right|$. We assume that there exists a subset ${\mathcal{E}}_p$ of $\mathcal{E}$ such that $\partial p={\cup }_{\sigma \in {\mathcal{E}}_p}\overline{\sigma }$ for all $p\in \mathcal{T}$.

3. $\mathcal{P}$ is a family of points of D indexed by $\mathcal{T}$ and $\mathcal{E}$, denoted by $\mathcal{P}=({\left({\mathbf{x}}_p\right)}_{p\in \mathcal{T}},{\left({\mathbf{x}}_{\sigma }\right)}_{\sigma \in \mathcal{E}})$, such that for all $p\in \mathcal{T}$, ${\mathbf{x}}_p\in p$ and for all $\sigma \in \mathcal{E}$, ${\mathbf{x}}_{\sigma }\in \sigma$. p is assumed to be ${\mathbf{x}}_p$-star-shaped, which means that there holds the inclusion $[{\mathbf{x}}_p,\mathbf{x}]\subset p$ for all $\mathbf{x}\in p$.

For all $p\in \mathcal{T}$ and $\sigma \in {\mathcal{E}}_p$, we denote by ${\mathbf{n}}_{p,\sigma }$ the unit vector normal to $\sigma$ outward to p, by $d_{p,\sigma }$ the Euclidean distance between ${\mathbf{x}}_p$ and the hyperplane including $\sigma$ and by $V_{p,\sigma }$ the cone with vertex ${\mathbf{x}}_p$ and basic $\sigma$. For all $\sigma \in \mathcal{E}$, we define ${\mathcal{T}}_{\sigma }=\{p\in \mathcal{T}:\sigma \in {\mathcal{E}}_p\}$. We denote by ${\mathcal{E}}_{int}$ the ensemble of all the inner edges.

Definition 2

(Time discretization). We define the time step $\delta t=T/N$ with $N\in {\mathbb{N}}^{+}$. In this way, $(0,T)={\cup }_{n=0}^{N-1}(n\delta t,(n+1)\delta t)$ and we set $t_n=n\delta t$ for all $n=0,1,2,\cdots ,N$.

In order to present the application of the SUSHI method, we formally integrate the two last equations of (7) on the domain $p\times (t_{n-1},t_n)$ for each $p\in \mathcal{T}$ and $n=1,\cdots ,N$ to obtain

$$\begin{array}{cc}\hfill & \theta {\int }_p\left(\varrho \left(w(\mathbf{x},t_n)\right)-\varrho \left(w(\mathbf{x},t_{n-1})\right)\right)d\mathbf{x}\hfill \\ \hfill & +\sum _{\sigma \in {\mathcal{E}}_p}{\int }_{t_{n-1}}^{t_n}{\int }_{\sigma }\mathbf{V}\varrho \left(w\right)·{\mathbf{n}}_{p,\sigma }d\gamma dt={\int }_{t_{n-1}}^{t_n}{\int }_p{\rho }_s{Q}_{s}d\mathbf{x}dt,\hfill \\ \hfill & \theta {\int }_p\{m\left(w(\mathbf{x},t_n)\right)-m\left(w(\mathbf{x},t_{n-1})\right)\}d\mathbf{x}+\sum _{\sigma \in {\mathcal{E}}_p}{\int }_{t_{n-1}}^{t_n}{\int }_{\sigma }\mathbf{V}m\left(w\right)·{\mathbf{n}}_{p,\sigma }d\gamma dt\hfill \\ \hfill & -\sum _{\sigma \in {\mathcal{E}}_p}{\int }_{t_{n-1}}^{t_n}{\int }_{\sigma }K\nabla w·{\mathbf{n}}_{p,\sigma }d\gamma dt={\int }_{t_{n-1}}^{t_n}{\int }_p{w}_s{\rho }_s{Q}_{s}d\mathbf{x}dt.\hfill \end{array}$$

We denote by $F_{p,\sigma }^D\left(w\right)$ the numerical flux that approximates the diffusion flux ${\displaystyle {\int }_{\sigma }-K\nabla w·{\mathbf{n}}_{p,\sigma }d\gamma }$. It is given by [5].

$$F_{p,\sigma }^D\left(w\right)=\sum _{{\sigma }^{\prime }\in {\mathcal{E}}_p}{A}_p^{\sigma {\sigma }^{\prime }}(w_p-w_{{\sigma }^{\prime }}),$$

(8)

where the matrices $A_p^{\sigma {\sigma }^{\prime }}$ are defined by

$$A_p^{\sigma {\sigma }^{\prime }}=\sum _{{\sigma }^{″}\in {\mathcal{E}}_p}{\mathbf{y}}^{{\sigma }^{″}\sigma }·K_{p,{\sigma }^{″}}{\mathbf{y}}^{{\sigma }^{″}{\sigma }^{\prime }}\mathrm{where}{K}_{p,{\sigma }^{″}}={\int }_{V_{p,{\sigma }^{″}}}Kd\mathbf{x},$$

and

$${\mathbf{y}}^{\sigma {\sigma }^{\prime }}=\left\{\begin{array}{cc}{\displaystyle \frac{\left|\sigma \right|}{\left|p\right|}{\mathbf{n}}_{p,\sigma }+\frac{\sqrt{d}}{d_{p,\sigma }}\left(1-\frac{\left|\sigma \right|}{\left|p\right|}{\mathbf{n}}_{p,\sigma }·({\mathbf{x}}_{\sigma }-{\mathbf{x}}_p)\right){\mathbf{n}}_{p,\sigma }}\hfill & \mathrm{if}\sigma ={\sigma }^{\prime }\hfill \\ {\displaystyle \frac{\left|{\sigma }^{\prime }\right|}{\left|p\right|}{\mathbf{n}}_{p,{\sigma }^{\prime }}-\frac{\sqrt{d}}{d_{p,\sigma }\left|p\right|}\left|{\sigma }^{\prime }\right|{\mathbf{n}}_{p,{\sigma }^{\prime }}·({\mathbf{x}}_{\sigma }-{\mathbf{x}}_p){\mathbf{n}}_{p,\sigma }}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(9)

where d is the space dimension.

Next, we denote by $F_{p,\sigma }^C(P,w)$ the approximation of the term ${\displaystyle {\int }_{\sigma }\mathbf{V}·{\mathbf{n}}_{p,\sigma }d\gamma }$. We apply an idea in [4] and set

$$F_{p,\sigma }^C(P,w)=\frac{k}{\nu }\left(\sum _{{\sigma }^{\prime }\in {\mathcal{E}}_p}\overline{A_p^{\sigma {\sigma }^{\prime }}}(P_p-P_{{\sigma }^{\prime }})-\left|\sigma \right|\varrho \left(w_{\sigma }\right)g\nabla z·{\mathbf{n}}_{p,\sigma }\right),$$

for all $\sigma \in {\mathcal{E}}_p$, where ${\displaystyle \overline{A_p^{\sigma {\sigma }^{\prime }}}=\sum _{{\sigma }^{″}\in {\mathcal{E}}_p}{\mathbf{y}}^{{\sigma }^{″}\sigma }·\overline{K_{p,{\sigma }^{″}}}{\mathbf{y}}^{{\sigma }^{″}{\sigma }^{\prime }}}$ with ${\displaystyle \overline{K_{p,{\sigma }^{″}}}={\int }_{V_{p,{\sigma }^{″}}}1d\mathbf{x}}$. As a result, we have the following approximations:

$${\int }_{\sigma }\mathbf{V}\varrho \left(w\right)·{\mathbf{n}}_{p,\sigma }d\gamma \approx F_{p,\sigma }^C(P,w)\varrho \left(w_{\sigma }\right),{\int }_{\sigma }\mathbf{V}m\left(w\right)·{\mathbf{n}}_{p,\sigma }d\gamma \approx F_{p,\sigma }^C(P,w)m\left(w_{\sigma }\right),$$

which will be discussed in more details in Section 3.2.

2.2. Adaptive Mesh Refinement

An essential feature of the generalized finite volume method SUSHI is that it can be applied on non-matching meshes. This allows us to apply an adaptive mesh method combining square or cubic elements of different sizes in the numerical tests. We refine the mesh in regions where the variation of the unknown function is large, and we merge the elements in regions where the variation of the unknown function is small. This reduces the number of elements and edges and saves CPU time. After each refinement or merge, we calculate the values of discrete unknown functions of the next time step on the new mesh.

In many articles, the value of the discrete gradients is chosen as a refinement criterion. In this article, we present a criterion based on the discrete difference of unknowns. The advantages of this method are: (1) it is easy to implement, (2) we can easily manage the number of unknowns as well as the refinement-remerge areas. Moreover, in some cases, the discrete gradient and discrete difference give the same results; for example, when the mesh size is close to uniform in the refinement area. However, this criterion is not based on a space error indicator.

We define the maximum value of the discrete differences between the unknown at the center and edges of the cell: ${\displaystyle U_p=\underset{\sigma \in {\mathcal{E}}_p}{max}|u_p-u_{\sigma }|}$, where u is an arbitrary unknown. Next, all elements are sorted in increasing order of $U_p$. If element q satisfies

$${\displaystyle (U_q-\underset{p\in \mathcal{T}}{min}{U}_p)}/{\displaystyle (\underset{p\in \mathcal{T}}{max}{U}_p-\underset{p\in \mathcal{T}}{min}{U}_p)}>\alpha ,$$

where $\alpha$ is a fixed number such that $\alpha \in (0,1)$, we added it to the “refined-list”. This means that the refined list will contain elements for which the unknowns undergo the largest discrete difference. The scalar $\alpha$ controls the number of refined elements (the length of the refined list). For example, if $\alpha$ = 0.75, then about a quarter of the elements in the current mesh will be added to the refined-list. If $\alpha$ is close to 1, the refined-list contains the elements with high discrete difference and the list will be short. If the element q satisfies

$${\displaystyle (U_q-\underset{p\in \mathcal{T}}{min}{U}_p)}/{\displaystyle (\underset{p\in \mathcal{T}}{max}{U}_p-\underset{p\in \mathcal{T}}{min}{U}_p)}<\beta ,$$

where $0<\beta <\alpha$, we added it to the “coarse-list”. Similarly, $\beta$ controls the length of the coarse-list. If $\beta$ tends to 0, the coarse-list will be empty.

The simulation starts with a uniform square mesh. After each time step, we calculate the corresponding refined list and coarse list. Each element that belongs to the refined-list will be divided into four small squares. Neighbor elements that belong to the coarse-list will be merged. We refer to [6] for detailed interpolation methods and we apply linear interpolation to assign values to new unknowns in each refined mesh. In practice, u corresponds to the salt concentration, except at the end of this paper where u is taken one time as the temperature.

3. Density Driven Flows in Porous Media

3.1. Initial Condition on P

In system (7), when u is considered as known, the second equation is an elliptic equation of P. Thus, we do not need any initial condition for P. We refer to [7] for a discussion of the compatibility relation between the pressure P and the density $\rho$.

However, in the discretized form of the problem, we will need to impose an initial condition for P, which we do by solving the elliptic equation $\nabla ·\left(\mathbf{V}\varrho \left(w_0\right)\right)=0$ together with the given boundary conditions provided that a Dirichlet boundary condition is imposed on a part of the boundary (Henry’s problem); otherwise, the initial pressure would not be uniquely defined.

In the rotating interface problem, due to the Neumann boundary condition on P, the solution P is not unique. This leads us to add an extra term $ϵ\partial P(\mathbf{x},t)/\partial t$ to the second equation of system (7) with $ϵ$ small. Therefore the equation for P becomes parabolic so that we have to impose an initial condition; we impose that the initial condition $P(\mathbf{x},0)$ satisfies that $P(\mathbf{x},0)=0$ at ${\mathbf{x}}_0=(0,0)$ and that $\nabla ·\left(\mathbf{V}\varrho \left(w_0\right)\right)=0$ for all $\mathbf{x}\in D\setminus {\mathbf{x}}_{\mathbf{0}}$ together with the homogenous Neumann boundary condition on $\partial D$. Then, the initial pressure is uniquely defined.

3.2. Numerical Scheme

The following discrete spaces is associated with the mesh

$$X_{\mathcal{D}}=\{({\left(v_p\right)}_{p\in \mathcal{T}},{\left(v_{\sigma }\right)}_{\sigma \in \mathcal{E}}),v_p\in \mathbb{R},v_{\sigma }\in \mathbb{R}\}.$$

We first discuss the discretization of the initial condition. As w is a function of u, the initial condition of the scheme for w is computed from the initial condition on u by

$${\displaystyle w_p^0=\frac{1}{\left|p\right|}{\int }_{p}w\left(u_0\left(\mathbf{x}\right)\right)d\mathbf{x},w_{\sigma }^0=\frac{1}{\left|\sigma \right|}{\int }_{\sigma }w(u_0\left(\mathbf{x}\right))d\gamma .}$$

As mentioned in Section 3.1, if $P_0\left(x\right)$ is given, we use a similar scheme to obtain the discrete initial condition

$${\displaystyle P_p^0=\frac{1}{\left|p\right|}{\int }_p{P}_0\left(\mathbf{x}\right)d\mathbf{x},P_{\sigma }^0=\frac{1}{\left|\sigma \right|}{\int }_{\sigma }{P}_0\left(\mathbf{x}\right)d\gamma .}$$

We numerically solve the equations described in the Section 3.1 to obtain the discrete initial condition for P in $X_{\mathcal{D}}$ for the rotating interface problem and Henry’s problem.

Next, we present a semi-implicit finite volume scheme for system (7), after plugging in the Darcy’s law into the 2 last equations. For each $n\in \{1,\cdots ,N\}$:

We suppose that $w^{n-1}$ and $w^{n-2}$ are already known and search for $P^n\in X_D$ such that

$$\left\{\begin{array}{}\mathrm{(10)}& \frac{\theta \left|p\right|}{\delta t}\left(\varrho \left(w_p^{n-1}\right)-\varrho \left(w_p^{n-2}\right)\right)+\sum _{\sigma \in {\mathcal{E}}_p}\{F_{p,\sigma }^C(P^n,w^{n-1})·\varrho \left(w_{\sigma }^{n-1}\right)\}={\displaystyle \frac{{\left(Q_{\rho }\right)}_p^n}{\delta t}}& \mathrm{for}\mathrm{all}p\in \mathcal{T}\mathrm{(11)}& \sum _{p\in {\mathcal{T}}_{\sigma }}\{F_{p,\sigma }^C(P^n,w^{n-1})·\varrho \left(w_{\sigma }^{n-1}\right)\}=0& \mathrm{for}\mathrm{all}\sigma \in {\mathcal{E}}_{int}\mathrm{(12)}& F_{P,\sigma }^C(P^n,w^{n-1})=\left|\sigma \right|{\overline{v}}_N({\mathbf{x}}_{\sigma },t_n)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_N^P\mathrm{(13)}& P_{\sigma }^n=P_D({\mathbf{x}}_{\sigma },t_n)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_D^P.\end{array}\right.$$

Then, knowing $P^n$ and $w^{n-1}$, we search for $w^n\in X_D$ such that

$$\left\{\begin{array}{}\mathrm{(14)}& \frac{\theta \left|p\right|}{\delta t}\left(m\left(w_p^n\right)-m\left(w_p^{n-1}\right)\right)+\sum _{\sigma \in {\mathcal{E}}_p}{F}_{p,\sigma }^D\left(w^n\right)+\sum _{\sigma \in {\mathcal{E}}_p}\{F_{p,\sigma }^C(P^n,w^{n-1})·m\left({\tilde{w}}_{p,\sigma }^n\right)\}={\displaystyle \frac{Q_p^n}{\delta t}}& \mathrm{for}\mathrm{all}p\in \mathcal{T}\mathrm{(15)}& \sum _{p\in {\mathcal{T}}_{\sigma }}\{F_{p,\sigma }^D\left(w^n\right)+F_{p,\sigma }^C(P^n,w^{n-1})·m\left({\tilde{w}}_{p,\sigma }^n\right)\}=0& \mathrm{for}\mathrm{all}\sigma \in {\mathcal{E}}_{int}\mathrm{(16)}& F_{p,\sigma }^D\left(w^n\right)=-\left|\sigma \right|D({\mathbf{x}}_{\sigma },t_n)\frac{\partial w}{\partial u}({\mathbf{x}}_{\sigma },t_n)·{\overline{u}}_N({\mathbf{x}}_{\sigma },t_n)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_N^u\mathrm{(17)}& w_{\sigma }=w\left(u_D({\mathbf{x}}_{\sigma },t_n)\right)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_D^u.\end{array}\right.$$

where ${\left(Q_{\rho }\right)}_p^n={\int }_{t_{n-1}}^{t_n}{\int }_p{\rho }_s{Q}_{s}d\mathbf{x}dt$ and $Q_p^n={\int }_{t_{n-1}}^{t_n}{\int }_p{w}_s{\rho }_s{Q}_{s}d\mathbf{x}dt$. We denote by $\partial D_D^P$ and $\partial D_D^u$ the parts of the boundary corresponding to the Dirichlet boundary conditions and by $\partial D_N^P$ and $\partial D_N^u$ the parts of the boundary corresponding to the Neumann boundary conditions. When $n=1$, we omit the term $\varrho \left(w_p^{n-1}\right)-\varrho \left(w_p^{n-2}\right)$ in the Equations (10). The Equations (11) and (15) are derived from the local conservation of the fluxes on the interior edges. We apply an upwind scheme for the convection term in (14), namely we define ${\tilde{w}}_{p,\sigma }$ according to the upwind scheme

$$m\left({\tilde{w}}_{p,\sigma }^n\right)=\left\{\begin{array}{cc}m\left(w_p^n\right)\hfill & \mathrm{if}{F}_{p,\sigma }^u(P^n,w^{n-1})>0\hfill \\ m\left(w_{\sigma }^n\right)\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

Moreover, we apply Newton’s method to calculate $w^n$ because in the discrete Equation (14), $m\left(w^n\right)$ is a nonlinear function of the unknown $w^n$.

3.3. The Rotating Interface Problem

The rotating interface problem involves a zero source term $Q_s=0$ together with the boundary conditions,

$$\left\{\begin{array}{cc}{\displaystyle \frac{\partial u}{\partial n}=0}\hfill & \mathrm{on} D\times (0,T)\hfill \\ \mathbf{V}·\mathbf{n}=0\hfill & \mathrm{on} D\times (0,T)\hfill \end{array}\right.$$

(18)

where $D=(0,100)\times (0,100){\mathrm{m}}^2$ and $T=500$ days.

As already mentioned in Section 3.1, we add the term $ϵ\partial P(\mathbf{x},t)/\partial t$ in the continuity equation for the fluid, which makes the model slightly compressible,

$$ϵ\frac{\partial P}{\partial t}+\theta \frac{\partial \varrho \left(w\right)}{\partial t}+\nabla ·\left(\mathbf{V}\varrho \left(w\right)\right)=0.$$

Then, the discretized form of the equation for P is given by:

$$\begin{array}{cc}\hfill \frac{ϵ\left|p\right|}{\delta t}(P_p^n-P_p^{n-1})& +\frac{\theta \left|p\right|}{\delta t}\left(\varrho \left(w_p^{n-1}\right)-\varrho \left(w_p^{n-2}\right)\right)\hfill \\ & +\sum _{\sigma \in {\mathcal{E}}_p}{F}_{p,\sigma }^C\left(P^n,m\left(w^{n-1}\right)\right)\varrho \left(w_{p,\sigma }^{n-1}\right)=0.\hfill \end{array}$$

(19)

The initial condition is such that the left half domain contains sea water and the right half domain contains fresh water, which corresponds to $u(\mathbf{x},0)=1$ if $x<L/2$ and $u(\mathbf{x},0)=0$ otherwise, where x is the x-coordinate, as presented in Figure 1. The initial condition for P has been discussed in Section 3.1.

The parameter values are defined in

Table 1. Parameters for the rotating interface problem.

Parameters	Value	Unit
Permeability k	$3.1\times {10}^{-12}$	m${}^2$
Dispersion-diffusion tensor K	$3.3\times {10}^{-6}$	m${}^2$/s
Gravity g	9.81	m/s${}^2$
Porosity $\theta$	0.5	-
Viscosity $\nu$	${10}^{-3}$	kg/(m·s)
Density of pure water ${\rho }_0$	${10}^3$	kg/m${}^3$
$\overline{a}$	0.3	-
$ϵ$	${10}^{-5}$	-

. Note that with the homogenous Neumann boundary condition (18), the conservation of mass holds as indicated in (21). The conservation of mass is discussed in the Section 3.4 below.

The salt water with higher concentration diffuses towards the bottom under the influence of the gravity. Meanwhile the fresh water diffuses towards the top. The interface between the sea water and fresh water evolves from vertical direction to horizontal direction. After about 300 days, the fluid is in equilibrium and the evolution stops. Figure 2 shows the flux field $\mathbf{V}$ at the initial time and at time $t=200$ days. One can observe the interface motion in Figure 3. The refined elements are located in the neighborhood of the interface and follow the movement of the interface motion.

We start with a uniform mesh of $5\times 5$ volume elements. An adaptive mesh is applied at each time step, taking into account the variations of the concentration u. We present the numerical results for the concentration in Figure 3. Since the mesh is coarse at the beginning, we choose $\alpha =0.75$ to create a finer mesh. Additionally, we choose $\beta =0.2$ to balance refinement and merge. At $t=50$ days, the number of elements is large enough, so we that set $\alpha =0.85$ to reduce the number of volume elements to be refined and set $\beta =0.1$. Since the evolution is slow from $t=250$ to the end, we choose $\alpha =0.95$ and $\beta \simeq 0$.

Table 2. Number of volume elements (NoE) and Number of unknowns (NoU) by using an adaptive mesh in the rotating interface test between fresh and sea water.

Time (Days)	NoE	NoU	Time (Days)	NoE	NoU
0	120	420
10	355	1123	100	934	2882
20	499	1567	200	1192	3660
30	583	1825	300	1309	4013
50	597	2171	500	1378	4221

shows the number of volume elements of the adaptive mesh at different times. In the early stage, the adaptive mesh has a small number of elements. At the end, the refined mesh in the interface region is approximately a uniform mesh of $40\times 40$.

In Figure 4, we compare the CPU times with different meshes with respect to the time in the physical flow process. We consider three uniform meshes with respectively 100 ($10\times 10$), 400 ($20\times 20$) and 1600 ($40\times 40$) volume elements and an adaptive mesh. We observe that using an adaptive mesh saves CPU time comparing to the uniform mesh $40\times 40$ case at the same flow process time.

3.4. Mass Conservation

We discuss the conservation of mass in this subsection. We first check the mass conservation in the partial differential equation system. We integrate the transport Equation (3) on $D\times (0,t)$ for a given t, in the case that $Q_s=0$:

$$\begin{array}{cc}\hfill & {\int }_D\theta \left(\rho \left(u(\mathbf{x},t)\right)u(\mathbf{x},t)-\rho \left(u_0\left(\mathbf{x}\right)\right)u_0\left(\mathbf{x}\right)\right)d\mathbf{x}\hfill \\ \hfill +& {\int }_0^t{\int }_D\nabla ·\left(\mathbf{V}\rho \left(u(\mathbf{x},\tau )\right)u(x,\tau )\right)d\mathbf{x}d\tau -{\int }_0^t{\int }_D\nabla ·(K\nabla w\left(u(\mathbf{x},\tau )\right))d\mathbf{x}d\tau =0,\hfill \end{array}$$

so that

$$\begin{array}{cc}\hfill & \theta \left({\int }_D\rho \left(u(\mathbf{x},t)\right)u(\mathbf{x},t)d\mathbf{x}-{\int }_D\rho \left(u_0\left(\mathbf{x}\right)\right)u_0\left(\mathbf{x}\right)d\mathbf{x}\right)\hfill \\ \hfill =& -{\int }_0^t{\int }_{\partial D}\rho \left(u(\mathbf{x},\tau )\right)u(\mathbf{x},\tau )\mathbf{V}·\mathbf{n}d\gamma d\tau +{\int }_0^t{\int }_{\partial D}K\rho \left(u(\mathbf{x},\tau )\right)\nabla u·\mathbf{n}d\gamma d\tau .\hfill \end{array}$$

(20)

We recall that that $\mathbf{V}·\mathbf{n}=0$ and $\partial u/\partial \mathbf{n}=0$ on $\partial D$, which we substitute in (20) to deduce

$${\int }_D\rho \left(u(\mathbf{x},t)\right)u(\mathbf{x},t)d\mathbf{x}={\int }_D\rho \left(u_0\left(\mathbf{x}\right)\right)u_0\left(\mathbf{x}\right)d\mathbf{x},$$

so that the system possesses the mass conservation property. Next, we check the mass is also conserved by the numerical scheme. We add up the semi-implicit discrete Equations (14) on all volume elements $p\in \mathcal{T}$, and set $Q_s=0$ to obtain: -4.6cm0cm

$$\sum _{p\in \mathcal{T}}\theta \left|p\right|\left(m\left(w_p^n\right)-m\left(w_p^{n-1}\right)\right)+\delta t\sum _{p\in \mathcal{T}}\sum _{\sigma \in {\mathcal{E}}_p}\left(F_{p,\sigma }^D\left(w^n\right)+F_{p,\sigma }^C\left(p^n,m\left(w^{n-1}\right)\right)m\left({\tilde{w}}_{p,\sigma }^n\right)\right)=0.$$

We define the total mass at the time $t_n$ by ${\displaystyle \sum _{p\in \mathcal{T}}\left|p\right|m\left(w_p^n\right)}$. Using the local conservation of the discrete fluxes on the interior edges (14) and the homogenous Neumann boundary condition, we deduce that:

$$\sum _{p\in \mathcal{T}}\left|p\right|m\left(w_p^n\right)=\sum _{p\in \mathcal{T}}\left|p\right|m\left(w_p^0\right),$$

(21)

which shows the fact that the mass is conserved in our algorithm. We check in

Table 3. Relative error of the total mass obtained by the generalized finite volume method SUSHI.

Time (Days)	Relative Error of Mass	Time (Days)	Relative Error of Mass
10	$2.73\times {10}^{-16}$	100	$1.50\times {10}^{-15}$
20	$5.47\times {10}^{-16}$	200	$6.83\times {10}^{-16}$
30	$2.73\times {10}^{-16}$	300	$3.83\times {10}^{-15}$
50	$8.20\times {10}^{-16}$	500	$6.01\times {10}^{-15}$

that, in the simulation, the mass is indeed conserved. The total mass at the initial time is $6.5\times {10}^6$.

3.5. Test Case Proposed by Henry

Henry’s problem describes the case that the sea water front advances in a confined aquifer which is initially filled with fresh water. It is one of the most standard tests for variable density flow in groundwater.

Mathematically, Henry’s problem [4] is defined as System (1) in space domain D together with the initial conditions such that $u(\mathbf{x},0)=0$ and that the pressure $P(\mathbf{x},0)$ satisfies $\nabla ·\left(\mathbf{V}\varrho \left(w\right(\mathbf{x},0\left)\right)\right)=0$ for all $\mathbf{x}\in D$. The source term is such that $Q_s=0$. Suppose that y is y-coordinate; then the boundary conditions are given by

$$\left\{\begin{array}{cc}u=0\hfill & \mathrm{on} {\mathsf{\Gamma }}_2\times (0,T)\hfill \\ u=1\hfill & \mathrm{on} {\mathsf{\Gamma }}_3^b\times (0,T)\hfill \\ {\displaystyle \frac{\partial u}{\partial n}=0}\hfill & \mathrm{on} ({\mathsf{\Gamma }}_1\cup {\mathsf{\Gamma }}_4\cup {\mathsf{\Gamma }}_3^h)\times (0,T)\hfill \\ P={\rho }_{0}g(\alpha (1-y)-y)\hfill & \mathrm{on} ({\mathsf{\Gamma }}_3^h\cup {\mathsf{\Gamma }}_3^b)\times (0,T)\hfill \\ \mathbf{V}·\mathbf{n}=0\hfill & \mathrm{on} ({\mathsf{\Gamma }}_1\cup {\mathsf{\Gamma }}_4)\times (0,T)\hfill \\ \mathbf{V}·\mathbf{n}=V_0\hfill & \mathrm{on} {\mathsf{\Gamma }}_2\times (0,T).\hfill \end{array}\right.$$

We consider the space domain $D=(0,1)\times (0,1)m^2$, and $T=0.05$ days. Figure 5 shows the configuration of the boundary conditions. The parameters are given in

Table 4. Parameters for Henry’s problem.

Parameters	Value	Unit
Permeability k	$1.02\times {10}^{-9}$	m${}^2$
Dispersion-diffusion tensor K	$6.6\times {10}^{-6}$	m${}^2$/s
Gravity g	9.81	m/s${}^2$
Porosity $\theta$	0.3	-
Viscosity $\nu$	${10}^{-3}$	kg/(m·s)
Density of pure water ${\rho }_0$	${10}^3$	kg/m${}^3$
$\overline{a}$	0.025	-
$V_0$	$-6.6\times {10}^{-5}$	m/d

.

Salt water invades from the right side, and fresh water with density ${\rho }_0$ flows in from the left side at a constant rate. Therefore, the concentration near the coast increases in time. The interface between freshwater and saltwater flowing from the lower right corner has a large slope at first (left in Figure 6), but gradually decreases as the salinity enters the domain (right in Figure 6). The simulation in Figure 7 shows how salt water invades the confined aquifer.

We start with a uniform mesh with $10\times 10$ square elements. The variation of the concentration u stays in a region around the boundary $x=1$, i.e., ${\mathsf{\Gamma }}_3$. We refine the mesh according to the variation of u and choose $\alpha =0.8$ to keep the refine-list long enough at each time step and set $\beta =0.15$ to keep the simulation process stable.

Table 5. Number of volume elements (NoE) and number of unknowns (NoU) by using an adaptive mesh in Henry’s problem.

Time (Days)	NoE	NoU	Time (Days)	NoE	NoU
0	120	420
$0.001$	262	844	$0.01$	469	1477
$0.0025$	307	983	$0.025$	586	1834
$0.003$	328	1046	$0.03$	649	2025
$0.005$	331	1053	$0.05$	673	2099

presents the number of elements of the adaptive mesh at various times. In the adaptive mesh, the smallest volume elements have the same diameter as the volume elements in the uniform mesh $40\times 40$. The number of elements is one of the reasons why the CPU time by using adaptive mesh is smaller than the CPU time corresponding to the uniform mesh $40\times 40$. The evolutions of the CPU times with respect to the time in the flow process in different meshes are shown in Figure 8. As a conclusion, by using an adaptive mesh, we keep the same precision over the high variation region and we save CPU time, especially in the case of Henry’s problem.

4. Density Driven Flow Coupled with Heat Transfer

In this section, we also introduce heat transfer. We perform computations for a test problem proposed in the SEAWAT documentation [8].

SEAWAT is a computer program for simulating the transport of solutes and heat [8]. It is a combination of MODFLOW and MT3DMS, which calculates the flow with MODFLOW and simulates solute transport with MT3DMS. SEAWAT can use both explicit and implicit time schemes to couple flow and solute transport.

The SEAWAT code is based on the finite difference method to solve systems of equations for variable density flow. MODFLOW and MT3DMS are both based on the finite difference method using a cell centered grid. The dependent variables obtained by the finite difference method represent the mean value in each cell (assuming that it is given at the center of the cell). Note that SEAWAT, MODFLOW, and MT3DMS use the same block center grid.

We propose a single code by using the generalized finite volume method SUSHI. The structure is simpler because it deals with one code rather than combining two separate codes with each other. Let us quote an article by [9] about the study of related problems using the Voronoi box-based finite volume method.

4.1. Partial Differential Equation System

4.1.1. Variable Density Groundwater Equation

We consider a model in SEAWAT and present below the case of a single species. The space domain is given by the rectangle $D=(0,L)\times (0,H)$. We apply the following variable density groundwater flow equation

$${\displaystyle S_s\rho \frac{\partial h}{\partial t}+\theta \frac{\partial \rho }{\partial u}\frac{\partial u}{\partial t}+\nabla ·\left(\mathbf{V}\rho \right)=q_s{\rho }_s}$$

(22)

in $D\times (0,T)$, where h [m] is the equivalent fresh water head, $S_s$ [1/m] is the specific storage, defined as the water volume released from the storage per unit decline of h, $\rho =\rho (h,u,\mathsf{\Theta })$ [kg/m${}^3$] is the fluid density, $\theta$ is the porosity, $q_s$ [d${}^{-1}$] is a source/sink of the fluid with density ${\rho }_s$. Equation (12) is a generalized form of Equation (1). The velocity $\mathbf{V}$ [m/d] is given by Darcy’s law:

$$\mathbf{V}(h,\nu ,\rho )=-{\displaystyle \frac{{\nu }_0}{\nu }}{\mathbf{K}}_0\left(\nabla h+{\displaystyle \frac{\rho -{\rho }_0}{{\rho }_0}}\nabla z\right),$$

(23)

where $\nu$ [kg/(m·s)] is the dynamic viscosity and ${\nu }_0$ is the reference viscosity, ${\mathbf{K}}_0$ [m/d] is the hydraulic conductivity tensor, ${\rho }_0$ is the density when the fluid is at the reference concentration $u_0$ [kg/m${}^3$] and the reference temperature ${\mathsf{\Theta }}_0$ [${}^{\circ }$C]. Additionally, z [m] is the elevation, such that $\nabla z={(0,1)}^T$.

4.1.2. The Fluid Density $\rho$

We refer to Hughes and Sanford [10] for the following equation for the transport of the solute species,

$$\rho (h,u,\mathsf{\Theta })={\rho }_0+{\displaystyle \frac{\partial \rho }{\partial u}}(u-u_0)+{\displaystyle \frac{\partial \rho }{\partial \mathsf{\Theta }}}(\mathsf{\Theta }-{\mathsf{\Theta }}_0)+{\displaystyle \frac{\partial \rho }{\partial P}}(P-P_0).$$

The concentration of the species can affect the fluid density. We reformulated here the term ${\displaystyle \frac{\partial \rho }{\partial P}}(P-P_0)$ by using the height of a water column l of density ${\rho }_0$. The variable l is related to the pressure by $P=l{\rho }_{0}g,$ where $l=h-z$. We rewrite the formula for the fluid density $\rho$ as a function of the concentration u, the temperature $\mathsf{\Theta }$ and the hydraulic head h as

$$\rho (h,u,\mathsf{\Theta })={\rho }_0+{\displaystyle \frac{\partial \rho }{\partial u}}(u-u_0)+{\displaystyle \frac{\partial \rho }{\partial \mathsf{\Theta }}}(\mathsf{\Theta }-{\mathsf{\Theta }}_0)+{\displaystyle \frac{\partial \rho }{\partial l}}(h-h_0).$$

(24)

An illustration of of the fluid density $\rho$ is given in Figure 9.

4.1.3. The Dynamic Fluid Viscosity $\nu$

In [3], the dynamic viscosity is considered as a function of temperature and the solute concentration is given. We ignore the dependence of the viscosity on the fluid pressure. A general formula for fluid viscosity as a function of concentration and temperature is given by

$$\nu (u,\mathsf{\Theta })={\nu }_0+{\displaystyle \frac{\partial \nu }{\partial u}}(u-u_0)+{\displaystyle \frac{\partial \nu }{\partial \mathsf{\Theta }}}(\mathsf{\Theta }-{\mathsf{\Theta }}_0).$$

In many temperature ranges, the linear approximation does not adequately represent the effect of temperature on dynamic viscosity. For this reason, we have introduced an alternative formula for dynamic viscosity, namely

$$\nu (u,\mathsf{\Theta })={\displaystyle \frac{\partial \nu }{\partial u}}(u-u_0)+{\nu }_{\mathsf{\Theta }}\left(\mathsf{\Theta }\right).$$

There are many ways to express ${\nu }_{\mathsf{\Theta }}\left(\mathsf{\Theta }\right)$ as a function of temperature. Here, the following relationship between viscosity and temperature is used.

$${\nu }_{\mathsf{\Theta }}\left(\mathsf{\Theta }\right)=A_1·A_2^{\left({\displaystyle \frac{A_3}{\mathsf{\Theta }+A_4}}\right)},$$

where $A_1$, $A_2$, $A_3$ and $A_4$ are positive constants (cf. C.I. Voss [11]). It yields the following formula for the dynamic viscosity:

$$\nu (u,\mathsf{\Theta })={\displaystyle \frac{\partial \nu }{\partial u}}(u-u_0)+A_1·A_2^{\left({\displaystyle \frac{A_3}{\mathsf{\Theta }+A_4}}\right)}.$$

(25)

An illustration of the fluid viscosity $\nu$ is given in Figure 10.

4.1.4. The Solute Transport Equation and the Heat Transport Equation

The solute transport in groundwater is described by an advection–dispersion equation. We apply the general form

$${\displaystyle (1+{\displaystyle \frac{{\rho }_b{K}_D^u}{\theta }})\frac{\partial \left(\theta u\right)}{\partial t}+\nabla ·\left(\mathbf{V}u\right)-\nabla ·\left((\theta D_m^u+\mathbf{a}·\mathbf{V})\nabla u\right)=q_s{u}_s}$$

(26)

in $D\times (0,T)$, where ${\rho }_b$ [kg/m³] is the bulk density, $K_D^u$ [m³/kg] is the distribution coefficient for salinity, $D_m^u$ [m²/d] is the diffusion coefficient and $\mathbf{a}$ [m] is the dispersivity tensor. We remark that the dispersion tensor is defined by

$$\mathbf{a}·\mathbf{V}:=(a_L-a_T)\frac{{\mathbf{V}}^T\otimes \mathbf{V}}{\left|\mathbf{V}\right|}+a_T\left|\mathbf{V}\right|\mathbf{I},$$

(27)

where $a_L$ and $a_T$ are given in

Table 6. Parameters of the model problem.

Parameter	Value	Unit
Specific storage $S_s$	$1.00\times {10}^{-5}$	1/m
Porosity $\theta$	$0.35$	−
Reference viscosity ${\nu }_0$	$0.001$	kg/(m·s)
Reference density ${\rho }_0$	1000	kg/m${}^3$
Bulk density ${\rho }_b$	$1761.5$	kg/m${}^3$
Reference concentration $u_0$	0	kg/m${}^3$
Reference temperature ${\mathsf{\Theta }}_0$	25	${}^{\circ }$C
Horizontal hydraulic conductivity $K_{h0}$	10	m/d
Vertical hydraulic conductivity $K_{v0}$	$0.1$	m/d
Longitudinal dispersivity $a_L$	1	m
Transverse dispersivity $a_T$	$0.1$	m
Diffusion coefficient $D_m^u$	$1.00\times {10}^{-10}$	m${}^2$/d
Bulk thermal diffusivity $D_m^{\mathsf{\Theta }}$	$0.150309621$	m${}^2$/d
Distribution coefficient for concentration $K_D^u$	0	m${}^3$/kg
Distribution coefficient for temperature $K_d^{\mathsf{\Theta }}$	$2.00\times {10}^{-4}$	m${}^3$/kg
$A_1$	$2.394\times {10}^{-5}$	−
$A_2$	10	−
$A_3$	$248.37$	−
$A_4$	$133.15$	−
$\delta \nu /\delta u$	$1.92\times {10}^{-6}$	m${}^2$/d
$\delta \rho /\delta h$	$4.46\times {10}^{-3}$	kg/m${}^4$
$\delta \rho /\delta u$	$0.7$	−
$\delta \rho /\delta \mathsf{\Theta }$	$-0.375$	kg/(m${}^3·{}^{\circ }$C)

.

Next, we apply the heat transport equation in Thorne et al., 2006 [12].

$${\displaystyle (1+{\displaystyle \frac{{\rho }_b{K}_d^{\mathsf{\Theta }}}{\theta }})\frac{\partial \left(\theta \mathsf{\Theta }\right)}{\partial t}+\nabla ·\left(\mathbf{V}\mathsf{\Theta }\right)-\nabla ·\left((\theta D_m^{\mathsf{\Theta }}+\mathbf{a}·\mathbf{V})\nabla \mathsf{\Theta }\right)=q_s{\mathsf{\Theta }}_s}$$

(28)

in $D\times (0,T)$, where $K_d^{\mathsf{\Theta }}$ [m³/kg] is the temperature distribution coefficient and $D_m^{\mathsf{\Theta }}$ [m²/d] is the bulk thermal diffusivity. We numerically solve the system of Equations (22), (23), (26) and (28) together with the boundary conditions

$$\left\{\begin{array}{cc}h=h_D(\mathbf{x},t)\hfill & \mathrm{on} D_D^h\times (0,T)\hfill \\ \mathbf{V}·\mathbf{n}={\overline{q}}_N(\mathbf{x},t)\hfill & \mathrm{on} D_N^h\times (0,T)\hfill \\ u=u_D(\mathbf{x},t)\hfill & \mathrm{on} D_D^u\times (0,T)\hfill \\ {\displaystyle \frac{\partial u}{\partial n}={\overline{u}}_N(\mathbf{x},t)}\hfill & \mathrm{on} D_N^u\times (0,T)\hfill \\ \mathsf{\Theta }={\mathsf{\Theta }}_D(\mathbf{x},t)\hfill & \mathrm{on} D_D^{\mathsf{\Theta }}\times (0,T)\hfill \\ {\displaystyle \frac{\partial \mathsf{\Theta }}{\partial n}={\overline{\mathsf{\Theta }}}_N(\mathbf{x},t)}\hfill & \mathrm{on} D_N^{\mathsf{\Theta }}\times (0,T)\hfill \end{array}\right.$$

(29)

where $\partial D=\overline{\partial D_D^h}\bigcup \overline{\partial D_N^h}=\overline{\partial D_D^u}\bigcup \overline{\partial D_N^u}=\overline{\partial D_D^{\mathsf{\Theta }}}\bigcup \overline{\partial D_N^{\mathsf{\Theta }}}$ and where $\partial D_D^h$, $\partial D_D^u$ and $\partial D_D^{\mathsf{\Theta }}$ correspond to Dirichlet boundary conditions and $\partial D_N^h$, $\partial D_N^u$ and $\partial D_N^{\mathsf{\Theta }}$ correspond to Neumann boundary conditions for the hydraulic head h, the concentration u, and the temperature $\mathsf{\Theta }$, respectively. We denote initial conditions by

$$\left\{\begin{array}{cc}h(\mathbf{x},0)=h_{ini}\left(\mathbf{x}\right)\hfill & \mathrm{in}D\hfill \\ u(\mathbf{x},0)=u_{ini}\left(\mathbf{x}\right)\hfill & \mathrm{in}D\hfill \\ \mathsf{\Theta }(\mathbf{x},0)={\mathsf{\Theta }}_{ini}\left(\mathbf{x}\right)\hfill & \mathrm{in}D.\hfill \end{array}\right.$$

(30)

4.2. Numerical Approximation

We remark that from [8], the derivatives $\frac{\partial \rho }{\partial h}$, $\frac{\partial \rho }{\partial u}$, $\frac{\partial \rho }{\partial \mathsf{\Theta }}$ in the Equation (24) and $\frac{\partial \nu }{\partial u}$ in the Equation (25) are constants. We set ${\theta }_u=\theta +{\rho }_b{K}_D^u$ and ${\theta }_{\mathsf{\Theta }}=\theta +{\rho }_b{K}_d^{\mathsf{\Theta }}$ and present below a numerical scheme corresponding to the system (22)–(28) with the boundary condition (29) and the initial conditions (30).

4.2.1. Numerical Approximation of the Fluxes

We plug Darcy’s law (23) into (22) and integrate over the volume element p for each $p\in \mathcal{T}$ to obtain

$$\begin{array}{cc}\hfill & S_s{\int }_p\rho {\displaystyle \frac{\partial h}{\partial t}}d\mathbf{x}+{\int }_p\theta \frac{\partial \rho }{\partial u}\frac{\partial u}{\partial t}d\mathbf{x}-\sum _{\sigma \in {\mathcal{E}}_p}{\int }_{\sigma }\rho {\displaystyle \frac{{\nu }_0}{\nu }}{\mathbf{K}}_0\nabla h·{\mathbf{n}}_{p,\sigma }d\gamma \hfill \\ & -\sum _{\sigma \in {\mathcal{E}}_p}{\int }_{\sigma }\rho {\displaystyle \frac{{\nu }_0}{\nu }}{\mathbf{K}}_0{\displaystyle \frac{\rho -{\rho }_0}{{\rho }_0}}\nabla z·{\mathbf{n}}_{p,\sigma }d\gamma \hfill \\ \hfill =& {\int }_p{q}_s{\rho }_{s}d\mathbf{x}.\hfill \end{array}$$

(31)

The diffusion flux ${\displaystyle -{\int }_{\sigma }\rho {\displaystyle \frac{{\nu }_0}{\nu }}{\mathbf{K}}_0\nabla h·{\mathbf{n}}_{p,\sigma }d\gamma }$ is approximated by ${\rho }_{\sigma }{F}_{p,\sigma }^h\left(h\right)$, such that

$${\rho }_{\sigma }{F}_{p,\sigma }^h\left(h\right)={\rho }_{\sigma }{\displaystyle \frac{{\nu }_0}{{\nu }_{\sigma }}}\sum _{{\sigma }^{\prime }\in {\mathcal{E}}_p}{\tilde{A}}_p^{\sigma {\sigma }^{\prime }}(h_p-h_{{\sigma }^{\prime }}),$$

(32)

for $h\in X_{\mathcal{D}}$ where ${\tilde{A}}_p^{\sigma {\sigma }^{\prime }}$ are defined as

$${\tilde{A}}_p^{\sigma {\sigma }^{\prime }}=\sum _{{\sigma }^{″}\in {\mathcal{E}}_p}{\mathbf{y}}^{{\sigma }^{″}\sigma }·\left({\int }_{V_{p,{\sigma }^{″}}}{\mathbf{K}}_{0}d\mathbf{x}\right){\mathbf{y}}^{{\sigma }^{″}{\sigma }^{\prime }},$$

and ${\mathbf{y}}^{\sigma {\sigma }^{\prime }}$ is given as in (9).

We integrate Equations (26) and (28) over the volume element p for each $p\in \mathcal{T}$, to obtain

$$\begin{array}{cc}\hfill & {\int }_p{\theta }_u\frac{\partial u}{\partial t}d\mathbf{x}+\sum _{\sigma \in {\mathcal{E}}_p}{\int }_{\sigma }\mathbf{V}u·{\mathbf{n}}_{p,\sigma }d\gamma -\sum _{\sigma \in {\mathcal{E}}_p}{\int }_{\sigma }(\theta D_m^u+\mathbf{a}·\mathbf{q})\nabla u·{\mathbf{n}}_{p,\sigma }d\gamma ={\int }_p{q}_s{u}_{s}d\mathbf{x},\hfill \\ \hfill & {\int }_p{\theta }_{\mathsf{\Theta }}\frac{\partial \mathsf{\Theta }}{\partial t}d\mathbf{x}+\sum _{\sigma \in {\mathcal{E}}_p}{\int }_{\sigma }\mathbf{V}\mathsf{\Theta }·{\mathbf{n}}_{p,\sigma }d\gamma -\sum _{\sigma \in {\mathcal{E}}_p}{\int }_{\sigma }(\theta D_m^{\mathsf{\Theta }}+\mathbf{a}·\mathbf{q})\nabla \mathsf{\Theta }·{\mathbf{n}}_{p,\sigma }d\gamma ={\int }_p{q}_s{\mathsf{\Theta }}_{s}d\mathbf{x}.\hfill \end{array}$$

We then define the numerical fluxes $F_{p,\sigma }^u\left(u\right)$ and $F_{p,\sigma }^{\mathsf{\Theta }}\left(\mathsf{\Theta }\right)$ to approximate the fluxes ${\displaystyle -{\int }_{\sigma }(\theta D_m^u+\mathbf{a}·\mathbf{q})\nabla u·{\mathbf{n}}_{p,\sigma }d\gamma }$ and ${\displaystyle -{\int }_{\sigma }(\theta D_m^{\mathsf{\Theta }}+\mathbf{a}·\mathbf{q})\nabla \mathsf{\Theta }·{\mathbf{n}}_{p,\sigma }d\gamma }$, respectively. In view of the SUSHI method, we present the discrete fluxes $F_{p,\sigma }^u\left(u\right)$ and $F_{p,\sigma }^{\mathsf{\Theta }}\left(\mathsf{\Theta }\right)$ as:

$$F_{p,\sigma }^u\left(u\right)=\sum _{{\sigma }^{\prime }\in {\mathcal{E}}_p}{\overline{A}}_p^{\sigma {\sigma }^{\prime }}(u_p-u_{{\sigma }^{\prime }}),$$

(33)

where

$${\overline{A}}_p^{\sigma {\sigma }^{\prime }}=\sum _{{\sigma }^{″}\in {\mathcal{E}}_p}{\mathbf{y}}^{{\sigma }^{″}\sigma }·\left({\int }_{V_{p,{\sigma }^{″}}}(\theta D_m^u+\mathbf{a}·\mathbf{q})d\mathbf{x}\right){\mathbf{y}}^{{\sigma }^{″}{\sigma }^{\prime }}$$

and

$$F_{p,\sigma }^{\mathsf{\Theta }}\left(\mathsf{\Theta }\right)=\sum _{{\sigma }^{\prime }\in {\mathcal{E}}_p}{\widehat{A}}_p^{\sigma {\sigma }^{\prime }}({\mathsf{\Theta }}_p-{\mathsf{\Theta }}_{{\sigma }^{\prime }}),$$

(34)

where

$${\widehat{A}}_p^{\sigma {\sigma }^{\prime }}=\sum _{{\sigma }^{″}\in {\mathcal{E}}_p}{\mathbf{y}}^{{\sigma }^{″}\sigma }·\left({\int }_{V_{p,{\sigma }^{″}}}(\theta D_m^{\mathsf{\Theta }}+\mathbf{a}·\mathbf{q})d\mathbf{x}\right){\mathbf{y}}^{{\sigma }^{″}{\sigma }^{\prime }}.$$

4.2.2. Numerical Scheme

The discrete initial conditions are given by

$${\displaystyle h_p^0=\frac{1}{\left|p\right|}{\int }_p{h}_{ini}\left(\mathbf{x}\right)d\mathbf{x},u_p^0=\frac{1}{\left|p\right|}{\int }_p{u}_{ini}\left(\mathbf{x}\right)d\mathbf{x},{\mathsf{\Theta }}_p^0=\frac{1}{\left|p\right|}{\int }_p{\mathsf{\Theta }}_{ini}\left(\mathbf{x}\right)d\mathbf{x}.}$$

Next, we present the discretized problem. For all $n\in \{1,...,N\}$:

We suppose that $h^{n-1}$, $u^{n-1}$, $u^{n-2}$ and ${\mathsf{\Theta }}^{n-1}$ are already known and search for $h^n\in X_D$ such that

$$\left\{\begin{array}{}\mathrm{(35)}& S_s\left|p\right|{\rho }_p^{n-1}(h_p^n-h_p^{n-1})+\theta \left|p\right|{\displaystyle \frac{\partial \rho }{\partial u}}(u_p^{n-1}-u_p^{n-2})+\delta t\sum _{\sigma \in {\mathcal{E}}_p}{F}_{p,\sigma }^h\left(h^n\right){\rho }_{\sigma }^{n-1}-\delta t\sum _{\sigma \in {\mathcal{E}}_p}\frac{{\nu }_0}{{\nu }_{\sigma }^{n-1}}\frac{{\rho }_{\sigma }^{n-1}-{\rho }_0}{{\rho }_0}{\rho }_{\sigma }^{n-1}{\mathbf{K}}_0\nabla z·{\mathbf{n}}_{p,\sigma }\left|\sigma \right|={\rho }_s{Q}_p^n& \mathrm{for}\mathrm{all}p\in \mathcal{T}\mathrm{(36)}& \sum _{p\in {\mathcal{T}}_{\sigma }}\left\{F_{p,\sigma }^h\left(h^n\right){\rho }_{\sigma }^{n-1}\right\}=0& \mathrm{for}\mathrm{all}\sigma \in {\mathcal{E}}_{int}\mathrm{(37)}& h_{\sigma }^n=h_D({\mathbf{x}}_{\sigma },t_n)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_D^h\mathrm{(38)}& V_{\sigma }^n=\left|\sigma \right|{\overline{V}}_N({\mathbf{x}}_{\sigma },t_n)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_N^h.\end{array}\right.$$

where ${\displaystyle Q_p^n={\int }_{t^{n-1}}^{t^n}{\int }_p{q}_{s}d\mathbf{x}dt}$, $F_{p,\sigma }^h\left(v\right)$, ${\rho }_{\sigma }^{n-1}$, and ${\nu }_{\sigma }^{n-1}$ are computed as (30), (40) and (41), respectively. When $n=1$, we omit the term $(u_p^{n-1}-u_p^{n-2})$ in Equation (35). The velocity $V_{\sigma }^n$ approximates ${\int }_{\sigma }\mathbf{V}·\mathbf{n}d\gamma \left(\mathbf{x}\right)$ at time $t=t_n$ and is given by

$$V_{\sigma }^n=F_{p,\sigma }^h\left(h^n\right)-{\displaystyle \frac{{\nu }_0}{{\nu }_{\sigma }^{n-1}}}{\displaystyle \frac{{\rho }_{\sigma }^{n-1}-{\rho }_0}{{\rho }_0}}{\mathbf{K}}_0\nabla z·{\mathbf{n}}_{p,\sigma }\left|\sigma \right|,$$

(39)

where

$${\rho }_{\sigma }^n={\rho }_0+{\displaystyle \frac{\partial \rho }{\partial u}}(u_{\sigma }^n-u_0)+{\displaystyle \frac{\partial \rho }{\partial \mathsf{\Theta }}}({\mathsf{\Theta }}_{\sigma }^n-{\mathsf{\Theta }}_0)+{\displaystyle \frac{\partial \rho }{\partial l}}(h_{\sigma }^n-h_0),$$

(40)

and

$${\nu }_{\sigma }^n={\displaystyle \frac{\partial \nu }{\partial u}}(u_{\sigma }^n-u_0)+A_1·A_2^{\left({\displaystyle \frac{A_3}{{\mathsf{\Theta }}_{\sigma }^n+A_4}}\right)}.$$

(41)

Knowing $V^n$, $u^{n-1}$, we search for $u^n\in X_D$ such that

$$\left\{\begin{array}{}\mathrm{(42)}& {\theta }_u\left|p\right|(u_p^n-u_p^{n-1})+\delta t\sum _{\sigma \in {\mathcal{E}}_p}{V}_{\sigma }^n{\tilde{u}}_{p,\sigma }^n+\delta t\sum _{\sigma \in {\mathcal{E}}_p}{F}_{p,\sigma }^u\left(u^n\right)=u_s{Q}_p^n& \mathrm{for}\mathrm{all}p\in \mathcal{T}\mathrm{(43)}& \sum _{p\in {\mathcal{T}}_{\sigma }}\left\{V_{\sigma }^n{\tilde{u}}_{p,\sigma }^n+F_{p,\sigma }^u\left(u^n\right)\right\}=0& \mathrm{for}\mathrm{all}\sigma \in {\mathcal{E}}_{int}\mathrm{(44)}& u_{\sigma }^n=u_D({\mathbf{x}}_{\sigma },t_n)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_D^u\mathrm{(45)}& F_{p,\sigma }^u\left(u^n\right)=-(\theta D_m^u+\mathbf{a}·\frac{V_{\sigma }^n}{\left|\sigma \right|}{\mathbf{n}}_{p,\sigma })\left|\sigma \right|{\overline{u}}_N({\mathbf{x}}_{\sigma },t_n)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_N^u.\end{array}\right.$$

Knowing $V^n$, ${\mathsf{\Theta }}^{n-1}$, we search for ${\mathsf{\Theta }}^n\in X_D$ such that

$$\left\{\begin{array}{}\mathrm{(46)}& {\theta }_{\mathsf{\Theta }}\left|p\right|({\mathsf{\Theta }}_p^n-{\mathsf{\Theta }}_p^{n-1})+\delta t\sum _{\sigma \in {\mathcal{E}}_p}{V}_{\sigma }^n{\tilde{\mathsf{\Theta }}}_{p,\sigma }^n+\delta t\sum _{\sigma \in {\mathcal{E}}_p}{F}_{p,\sigma }^{\mathsf{\Theta }}\left({\mathsf{\Theta }}^n\right)={\mathsf{\Theta }}_s{Q}_p^n& \mathrm{for}\mathrm{all}p\in \mathcal{T}\mathrm{(47)}& \sum _{p\in {\mathcal{T}}_{\sigma }}\left\{V_{\sigma }^n{\tilde{\mathsf{\Theta }}}_{p,\sigma }^n+F_{p,\sigma }^{\mathsf{\Theta }}\left({\mathsf{\Theta }}^n\right)\right\}=0& \mathrm{for}\mathrm{all}\sigma \in {\mathcal{E}}_{int}\mathrm{(48)}& {\mathsf{\Theta }}_{\sigma }^n={\mathsf{\Theta }}_D({\mathbf{x}}_{\sigma },t_n)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_D^{\mathsf{\Theta }}\mathrm{(49)}& F_{p,\sigma }^{\mathsf{\Theta }}\left({\mathsf{\Theta }}^n\right)=-(\theta D_m^{\mathsf{\Theta }}+\mathbf{a}·\frac{V_{\sigma }^n}{\left|\sigma \right|}{\mathbf{n}}_{p,\sigma })\left|\sigma \right|{\overline{\mathsf{\Theta }}}_N({\mathbf{x}}_{\sigma },t_n)& \mathrm{for}\mathrm{all}{\mathbf{x}}_{\sigma }\in \partial D_N^{\mathsf{\Theta }}.\end{array}\right.$$

where ${\tilde{u}}_{p,\sigma }^n$ and ${\tilde{\mathsf{\Theta }}}_{p,\sigma }^n$ are given by the upwind scheme

$${\tilde{u}}_{p,\sigma }^n=\left\{\begin{array}{cc}{u}_p^n\hfill & \mathrm{if}{V}_{\sigma }^n>0\hfill \\ u_{\sigma }^n\hfill & \mathrm{otherwise}\hfill \end{array}\right.\mathrm{and}{\tilde{\mathsf{\Theta }}}_{p,\sigma }^n=\left\{\begin{array}{cc}{\mathsf{\Theta }}_p^n\hfill & \mathrm{if}{V}_{\sigma }^n>0\hfill \\ {\mathsf{\Theta }}_{\sigma }^n\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

for each time step n.

The discrete fluxes $F_{P,\sigma }^u\left(u\right)$ and $F_{p,\sigma }^{\mathsf{\Theta }}\left(\mathsf{\Theta }\right)$ are defined by (33) and (34), respectively. The three Equations (36), (43) and (47) are the conservation of the discrete fluxes on the interior edges.

4.3. Numerical Tests

We consider the model problem developed in the SEAWAT [8] that is a two-dimensional cross section of a confined coastal aquifer initially filled with seawater at temperature 5 ${}^{\circ }$C, which corresponds to $h_{ini}=0$, $u_{ini}=35$ and ${\mathsf{\Theta }}_{ini}=5$. Along the left boundary, warm freshwater with a temperature of $25{}^{\circ }\mathrm{C}$ is injected into the coastal aquifer to represent the flow from the inland area. Warm freshwater flows to the right and flows into a vertical ocean boundary. At the ocean boundary, hydrostatic pressure conditions based on the fluid density calculated from the salinity of seawater at 5 ${}^{\circ }$C were set. No flow boundary conditions are set for the upper and lower boundaries. Because of the gravity and the no flow boundary conditions the velocities differ at the top, bottom, and middle of the aquifer. At the bottom part, the flow is slow, so the concentration and temperature values do not change much. At the top part, the flow is fast, and the concentration and temperature values change rapidly.

The boundary condition on the sea water boundary is given by $h(x=L,z)=(H-z)·\frac{{\rho }_{seawater}-{\rho }_0}{{\rho }_0}$, where ${\rho }_{seawater}$ is the seawater density and the elevation z is such that $0\le z\le H$. The average flux velocity on the boundary $(x=0,z)$ is given by ${\overline{V}}_N=\frac{Q}{H}$ with $Q=10m^2/d$. We have that $q_s=0$, and the parameters are given in Table 6.

The adaptive mesh is efficient in this problem, since the variations of the concentration and of the temperature only take place in their interface areas, respectively. We present numerical simulations in three cases. For test cases 1 and 2 below, we neglect the heat transfer equation and consider the mesh refinement to be based on the variation of the concentration u (Figure 11 and Figure 12). In test case 3, the interfaces of the concentration and of the temperature do not advance with the same speed, so that we perform a first computation with the refinement based on the variation of the concentration (Figure 13) and a second time computation with the refinement based on the variation of the temperature in order to present the results on the temperature profiles (Figure 14).

We start from a uniform $40\times 20$ square mesh to discretize the domain $(0,L)\times (0,H)$ where $L=2000$ and $H=1000$. We fix $\alpha =0.8$ and $\beta =0.2$ to keep the simulation process stable; we observe that the mesh is mostly refined in the neighborhood of the interface.

4.3.1. Test Cases 1 and 2 without Heat Transfer

In the first 2 test cases, we consider different expressions of the fluid density $\rho$ and suppose that the fluid viscosity $\nu$ is constant. We neglect the heat transfer equation and consider the system

$$\left\{\begin{array}{c}{S}_s\rho \frac{\partial h}{\partial t}+\theta \frac{\partial \rho }{\partial u}\frac{\partial u}{\partial t}+\nabla ·\left(\mathbf{V}\rho \right)=q_s{\rho }_s,\hfill \\ (1+{\displaystyle \frac{{\rho }_b{K}_D^u}{\theta }})\frac{\partial \left(\theta u\right)}{\partial t}+\nabla ·\left(\mathbf{V}u\right)-\nabla ·\left((\theta D_m^u+\mathbf{a}·\mathbf{V})\nabla u\right)=q_s{u}_s.\hfill \end{array}\right.$$

(50)

We apply the Equations (35), (39) and (42) for the numerical simulations.

Test case 1: We suppose that the fluid density only depends on the concentration u and can be rewritten as

$$\rho ={\rho }_0+\frac{\partial \rho }{\partial u}(u-u_0).$$

Figure 11 shows the evolution of the concentration u in time in test case 1.

Test case 2: In test case 2, we suppose that the density $\rho$ is a linear function of concentration u and of temperature $\mathsf{\Theta }$:

$$\rho ={\rho }_0+\frac{\partial \rho }{\partial u}(u-u_0)+\frac{\partial \rho }{\partial \mathsf{\Theta }}(\mathsf{\Theta }-{\mathsf{\Theta }}_0).$$

(51)

${\displaystyle \frac{\partial \rho }{\partial \mathsf{\Theta }}}$ equals to $-0.375$ and it indicates that the fluid density will decrease as the temperature increases.We suppose that the temperature $\mathsf{\Theta }$ is a linear function of the concentration u:

$$\begin{array}{cc}\hfill \mathsf{\Theta }& =\left(\frac{{\mathsf{\Theta }}_0-{\mathsf{\Theta }}_{ocean}}{u_0-u_{ocean}}\right)u+({\mathsf{\Theta }}_0-{\mathsf{\Theta }}_{ocean})\hfill \\ & =-\frac{4}{7}u+20,\hfill \end{array}$$

where ${\mathsf{\Theta }}_0={25}^{\circ }C$, ${\mathsf{\Theta }}_{ocean}=5^{\circ }C$, $u_0=0$ kg/m${}^3$ and $u_{ocean}=35$ kg/m${}^3$. Figure 12 shows the evolution in time of the concentration u in this test case.

We compare our results to those in SEAWAT [8] in the Figure 11 and Figure 12, and by a uniform piece of code, we obtain results which look very close to the results of SEAWAT. We then compare the numerical results in the two cases and find out that in test case 2, the front advancement speed is lower and that the transition front’s width is thinner due to the influence of the temperature.

4.3.2. Test Case 3

In test case 3, we take the heat transfer equation into account. In this case, the fluid density $\rho$ is given by (14) and the viscosity $\nu$ is given by (15). The transient motion of the concentration and that of the temperature are shown in the Figure 13 and Figure 14. In Figure 13, an adaptive mesh corresponding to the concentration variation is used. Meanwhile, we perform the simulation for a second time with an adaptive mesh corresponding to the temperature in Figure 14.

We observe that the concentration transition zone and temperature transition zone do not advance at the same speed; and the width of the transition fronts of the concentration u is thicker in test case 3 than in the previous 2 cases.

5. Conclusions

In this article, the generalized finite volume method SUSHI on an adaptive mesh is applied to simulate density-driven flow in a porous medium.

HydroExpert, the small hydrology company with whom we used to work, only used square mesh elements for the space discretization; these elements coincided with data coming from physical measurements. This led us to also apply a standard square mesh in space dimension 2. We proposed a criterion for selecting the volume elements to be refined based on the variability of selected geophysical quantities, such as the concentration. It turns out that the mesh is mostly refined along the interface between fresh water and sea water so that we save CPU time while maintaining the same accuracy in interface regions.

Refining any part of the domain will result in a non-matching mesh. The finite difference method or the standard finite volume method cannot be applied to the non-matching mesh, but the SUSHI method can be applied. Moreover, the SUSHI method can handle problems with anisotropic and inhomogeneous diffusion tensors.

In the simulations of the SEAWAT model, we proposed a single algorithm that is simpler than the SEAWAT software, which combines MODFLOW to compute the flow with MT3DMS to simulate the solute-transport. Taking all this into account, the method that we propose is an effective option for simulations related to density driven flow in porous media.