Modeling of the Saltwater Intrusion Using the Level Set Method. Application to Henry’s Problem

Joachna Meya Loua-Bouayi, Christian Tathy, Adolphe Kimbonguila Manounou

Laboratoire Mécanique Energétique et Ingénierie, Ecole Nationale Supérieure Polytechnique, Université Marien Ngouabi, Brazzaville, Congo.

DOI: 10.4236/cweee.2022.111002                

Abstract

The salt intrusion phenomenon is caused by overexploitation of aquifers in coastal areas. This physical phenomenon has been the subject of numerous studies and numerous methods have been proposed, with the aim of protecting the quality of the water in these aquifers. This work proposes a two-dimensional saline intrusion model using the sharp interface approach and the level set method. It consists of a parabolic equation modeling the underground flow and a hyperbolic Equation (the level set equation) which makes it possible to track the evolution of the interface. High-order numerical schemes such as the space scheme WENO5 and the third-order time scheme TVD-RK were used for the numerical resolution of the hyperbolic equation. To limit the tightening of the contour curves of the level set function, the redistanciation or reinitialization algorithm proposed by Sussma et al. (1994) was used. To ensure the effectiveness and reliability of the proposed method, two tests relating to the standard Henry problem and the modified Henry problem were performed. Recall that Henry’s problem uses the variable density modeling approach in a confined and homogeneous aquifer. By comparing the results obtained by the level set method with reinitialization (LSMR) and those obtained by Henry (1964), and by Simpson and Clement (2004), we see in the two test cases that the level set method reproduces well the toe, the tip and the behaviour of the interface. These results correspond to the results obtained by Abarca for Henry’s problem with constant dispersion coefficients. The results obtained with LSMR, reproduced the interface with a slight spacing compared to those obtained by Henry. According to Abarca (2006), this spacing is due to the absence of the longitudinal and transversal dispersion coefficients in the model.

Keywords

Confined Aquifer, Seawater Intrusion, Level Set Method, Sharp Interface, Henry Problem

1. Introduction

The issue of the reserves and access to water is one of the major problem’s humanity would face for decades to come. Today one estimates that one inhabitant out of five of the planets does not have access to sufficient water, and one out of three to good quality water [1]. Although it is abundant on the planet, 97.2% of the water on earth is salty and therefore unfit for human consumption, and 2.8% of the water on earth is freshwater, distributed as follows: 68% in glaciers and snows; 29.9% in subsoils; 0.3% in surface freshwater (river, lakes, humidity, atmosphere, …) and the rest in water frozen in the ground. However, out of 0.3% of surface freshwater, 80% evaporate permanently, and the rest is difficult to access. This quantity of freshwater is continuously renewed by precipitation [2].

Groundwater is the main source of water for household consumption, irrigation and the food industry. Its storage in the layers of the subsoil, sometimes at very great depths, preserves its quality. In addition, it does not require large investments for its treatment, as is the case for water collected from the surface.

Because of the use of pesticides and fertilizers in agricultural land, this resource becomes fragile and vulnerable to pollution [3]. In addition, changes in hydrological regimes, rising of sea levels and anthropogenic activities, cause the transfer of pollutants into the subsoil, and damage the quality and quantity of the resource [4] [5]. This can make it unfit for consumption, even inaccessible to the most vulnerable populations. However, decontamination is possible. Unfortunately, it can be very costly.

Many regions of the world exploit groundwater for daily water requirements of households, industry and agriculture. Several studies have shown risks of degradation of the quality of the aquifers to overexploitation [6] [7]. In the case of coastal zones, the overexploitation of groundwater and sea level rise due to global warming can cause infiltration of oceanic saltwater into coastal aquifers, leading to the physical phenomenon of saline water or marine intrusion. The study carried out by Sherif [8] on different scenarios for adding pumping wells in the Nile Delta aquifer region in Egypt, shows that by reducing pumping wells in the zones of high-risk, it is possible to decelerate and minimize the process of saline intrusion into the aquifers. Thus, to ensure optimal use of freshwater and to monitor marine intrusion in coastal aquifers, Rifai et al. [9], Essaid [10], Ogata [11], Hubbert [12], Sark [13], have proposed mathematical models to track the interface between freshwater and saltwater, to strengthen management tools for the operators of aquifers and to preserve the quality of the resource.

The phenomenon of saline intrusion has been the subject of numerous studies, and numerous formulations have been proposed; let us quote among others the work carried out by Bardon-Ghyben [14], Herzberg [15], Hubbert [12], Henry [16], Bear and Verruijt [17], etc. Two large families of models have been developed. The first is the approach with a sharp or abrupt interface which considers two immiscible fluids with an interface separating them. The second is the variable density or density-dependent approach based on the salt transport equation and flow equation of miscible fluids.

The sharp interface approach assumes that the two fluids are immiscible; hence the hypothesis of the existence of an interface between freshwater and saltwater [12] [16]. This hypothesis depends on the characteristics of the aquifer and particle movements due to tides or recharge fluctuations [18]. In the case of this approach, it is assumed that the thickness of the transition zone due to the diffusion of saltwater in freshwater is relatively small compared to the size of the aquifer. These models have proven to be robust and reliable in the absence of important exterior forcing terms. In this context, one can necessarily cite the work carried out by Essaid [10], Fetter [19], Mualem and Bear [20], etc.

The variable density approach is based on the existence of a transition zone that contains a mixture of freshwater and saltwater [9] [11]. It involves a transport equation of convection-dispersion of salt in the aquifer. The variable density approach is the most realistic describing the phenomenon of saline intrusion, because the two fluids (freshwater and saltwater) are miscible. Mathematically, the equations developed in this approach take the form of a system of nonlinear strongly coupled, parabolic partial differential equations, which makes it heavier and more costly in computations both from the analytical and numerical point of view [21] [22]. In addition, in the case of an unconfined aquifer, this approach presents difficulties to define the zone of desaturation in the saturated and unsaturated zone with water.

The saltwater intrusion has been the subject of comparative studies using numerical, analytical and experimental formulations; let us cite for example the work carried out by Rifai et al. [9], Ogata [11] and Cheng and Ouazar [6] which present studies of this problem in different countries and cities, as well as numerical methods developed to solve them, including ModFlow, Seawat and Sharp softwares. The studies carried out by Mory [23] and Abudawia [24], show that it is also possible to combine the two different approaches to take into account certain assumptions neglected in the two previous approaches. This approach was deduced from the phase-field theory, introduced by Allen-Cahn to describe the transition phenomena between two fluids.

The use of the approach with a sharp interface imposes to take into account the conditions at the interface. The flow equations are either solved respectively in the freshwater phase and the saltwater phase [25], either solved in the freshwater phase assuming that the saltwater is hydrostatic [26]. The freshwater/saltwater interface is deduced by considering the pressure continuity between the two phases. Models based on the sharp interface approach do not describe the behavior of the transition zone, but they give information on the movement of the salt bevel (toe).

The problem solving in the abrupt interface approach uses twice the flow equation in the different phases, the objective of this work is to use only one flow equation modeling two phases (freshwater and saltwater) thanks to a level set function, solution of a first-order hyperbolic transport equation, whose zero-contour line describes the freshwater/saltwater interface.

To do that, one develops a two-dimensional method to solve the saline intrusion problem in a coastal aquifer, using the sharp interface approach based on the level set method, [27] - [31]. Implementation of the level set method requires the use of robust and accurate numerical schemes in order to avoid significant mass loss.

2. Formulation of the Problem

2.1. The Equations Modeling the Salt Intrusion Problem

The equations governing the salt intrusion phenomenon are based on the association of the continuity equation with Darcy’s Law. In the case of a confined aquifer, these equations are presented as follows [7] [17] [32] [33]:

$\{\begin{array}{l}{S}_{\alpha }{\rho }_{\alpha }\frac{\partial H_{\alpha }}{\partial t}+\nabla \cdot \left({\rho }_{\alpha }{v}_{d_{\alpha }}\right)={\rho }_{\alpha }Q\text{with}\alpha =f,s\\ v_{d_{\alpha }}=-K_{\alpha }\nabla H_{\alpha }\end{array}$ (1)

where $S_{\alpha }$, $H_{\alpha }$, $K_{\alpha }$, $Q$ and $v_{d_{\alpha }}$ are respectively the specific storage coefficient, the hydraulic head, the hydraulic conductivity, the source term and the Darcy velocity in the different phases $\alpha$ ( $\alpha =f$ represents the freshwater phase; $\alpha =s$ represents the saltwater phase).

Given the difficulties listed above for the variable density approach, on the implementation of analytical and numerical methods, it may be easier to use the sharp interface approach. Unfortunately, the use of interface continuity conditions in the sharp interface approach, gives coupled and nonlinear flow equations, therefore quite difficult to resolve.

In the tracking interface method used in this work, the pressure continuity condition at the interface between the two phases is taken into account in Darcy’s law through the surface tension coefficient.

2.2. The Level Set Method

Introduced by Osher and Sethian [30], the level set method was developed for the problems of tracking interfaces represented by closed curves or surfaces. It has been used in several works and has shown its effectiveness for the problems of two-phase flows of which one wishes to track the interface between the two phases [34] [35]. In addition, it has potential in terms of possible extensions, improvements and perspectives [36].

The basic idea of the level set method is to build a scalar function (level set), whose zero contour line or zero isocontours is the interface that one seeks to describe, i.e. to define an interface $\Gamma$ which bounds an open domain ${\Omega }^{-}$. This interface is defined at a point $x$ of a domain $\Omega ={\Omega }^{-}\cup \Gamma \cup {\Omega }^{+}$ at time t, like the zero isocontours of a scalar function $\varphi \left(x,t\right)$ (Figure 1).

In 2D the evolution interface $\Gamma$ separating the two sub-domains ${\Omega }^{-}$ and ${\Omega }^{+}$ is defined by:

$\Gamma \left(t\right)=\left\{x\in {ℝ}^2/\varphi \left(x,t\right)=0\right\}.$ (2)

The level set function $\varphi$ is negative in ${\Omega }^{-}$, positive outside i.e. in ${\Omega }^{+}$ and zero on $\Gamma \left(t\right)$.

For two immiscible fluids, the domain $\Omega$ can also be written as the union of three sub-domains, $\Omega ={\Omega }_f\cup {\Omega }_s\cup \Gamma$, by using the sign of the level set function, i.e.:

$\{\begin{array}{l}\varphi \left(x,t\right)<0\text{for}x\in {\Omega }_f={\Omega }^{-}\\ \varphi \left(x,t\right)>0\text{for}x\in {\Omega }_s={\Omega }^{+}\\ \varphi \left(x,t\right)=0\text{for}x\in \Gamma \end{array}$

The level set function can be perceived as an infinite variety of isocontours. When the distance d is zero, the isocontour is the interface itself. In this the isocontour is defined as:

$\{\begin{array}{l}\varphi \left(x,t\right)=d\left(x,x_{\Gamma }\right)\text{for}x\in {\Omega }_s,x_{\Gamma }\in \Gamma \\ \varphi \left(x,t\right)=-d\left(x,x_{\Gamma }\right)\text{for}x\in {\Omega }_f,x_{\Gamma }\in \Gamma \\ \varphi \left(x,t\right)=0\text{for}x\in \Gamma \end{array}$ (3)

where d is the signed distance function.

The initial position ${\varphi }_o$ of the interface being known, its evolution is calculated using a given velocity field $v$. These velocities can be functions of space, time, geometry of the interface or physics of the external environment. Thus, the transport equation making it possible to describe the evolution of the interface subjected to a velocity field $v$ is:

Figure 1. Presentation of the $\Omega$ domain.

${\partial }_t\varphi +\nabla \cdot \left(v\varphi \right)=0$ (4)

with a given initial condition $\varphi \left(x,t=0\right)={\varphi }_o\left(x\right)$ and some boundary conditions.

The level set method has several advantages, among which the taking into account of the topology and the easy calculation of geometric characteristics (curvature, normal). However, it also suffers from a few disadvantages namely, the non-conservation of the mass due to the numerical errors on the resolution of the transport equation, and the spacing or tightening of the level lines due to the sheared velocity fields resulting in the loss of accuracy of the geometric characteristics of the interface i.e. $\Vert \nabla \varphi \left(x,t\right)\Vert \ne 1$.

The distance property $\Vert \nabla \varphi \left(x,t\right)\Vert =1$ is important for signed distance functions, it is generally the basis of the level set method. In the case of loss of this property i.e. if $\Vert \nabla \varphi \left(x,t\right)\Vert \ne 1$, the method then becomes imprecise. This is why it is necessary to impose a condition on the evolution of $\varphi$ in $\Omega$ called a reinitialization or redistanciation algorithm, that is to say $\Vert \nabla \varphi \left(x,t\right)\Vert =1$, in order to correct the level lines, without the interface moving.

This reinitialization algorithm was developed by Sussman, Smereka and Osher [37]. Its purpose is to iteratively correct the position of the level lines, starting from the single zero level line which is an interface. This algorithm is based on the following evolution equation:

$\{\begin{array}{l}\frac{\partial d\left(x,\tau \right)}{\partial \tau }=\text{sign}\left(\varphi \left(x,t\right)\right)\left(1-\Vert d\left(x,\tau \right)\Vert \right)\\ d\left(x,\tau =0\right)=\varphi \left(x,t\right)\end{array}$ (5)

where $\tau$ represents a fictitious reinitialization time of the function $\varphi$, d represents the level set function at fictitious time $\tau$, and the $\text{sign}(\cdot )$ function represents the sign of the level set function $\varphi$. This function is smoothed around the interface and is defined as follow:

$\text{sign}\left(\varphi \right)=\{\begin{array}{l}-1\text{if}\varphi <-\Delta x\\ \frac{\varphi }{\sqrt{{\varphi }^2+\min {\left(\Delta x,\Delta z\right)}^2}}\text{if}\left|\varphi \right|\le \Delta x\\ 1\text{if}\varphi >\Delta x\end{array}$ (6)

where $\Delta x$ and $\Delta z$ are the steps of space discretization. At each time step $\Delta t$, the Equation (5) is initialized and resolved according to a fictitious time $\tau$, until the steady-state is reached. At the end of this iterative process, the property of distance function is found for the variable d; one then reinitializes the level set function in all the domain $\Omega$ :

$\forall x\in \Omega ,\varphi \left(x,t\right)=d\left(x,\tau ={\tau }_{final}\right)$

2.3. Coupling of the Flow Equation with the Level Set Method

The hypothesis of two immiscible fluids leads to the notion of interface, which is the locus of discontinuity for some variables. So, the interface we want to locate by an appropriate method is a zone of discontinuities. Hence the need to add additional conditions in the model, which connect the two fluids at the interface, called jump conditions at the interface. These jump conditions require implementation of efficient numerical tools to connect the solutions at the interface between the two phases. The method used here for such connection is the Continuum Surface Force (CSF) method also called Delta Function Method (DFM), proposed by Brackbill et al. [38] and described below. It models the surface tension between two immiscible fluids. This method does not really vanish the discontinuities, but it smooths them artificially by using the Heaviside and Dirac functions.

The main advantage of the CSF method is to construct a single flow equation whose parameters depend on the level set function for which the sign makes it possible to distinguish the different freshwater/saltwater phases. This method is robust and ease to implement. However, with such an approach, instability of the interface can occur when the surface tension forces are very strong [39].

Thus, the CSF method will make it possible to associate the flow equation with the transport equation making it possible to follow the interface in time and space. In the case of a confined aquifer, these equations are given as follows:

$\{\begin{array}{l}{S}_{\alpha }{\rho }_{\alpha }\frac{\partial H_{\alpha }}{\partial t}+\nabla \cdot \left({\rho }_{\alpha }{v}_{d_{\alpha }}\right)={\rho }_{\alpha }Q\text{with}\alpha =f,s\\ v_{d_{\alpha }}=-K_{\alpha }\nabla H_{\alpha }\end{array}$ (7)

and

$\{\begin{array}{l}\frac{\partial \varphi }{\partial t}+\nabla \cdot \left(v\varphi \right)=0\text{in}\Omega \times \left[0,T\right]\\ \varphi \left(x,t=0\right)={\varphi }_o\left(x\right)\text{in}\Omega \\ \nabla \varphi \cdot n=0\text{in}\partial \Omega \times \left[0,T\right]\end{array}$ (8)

with ${\varphi }_o$ a given function and $n$ the unit normal vector at the interface going from ${\Omega }_s$ to ${\Omega }_f$.

In fluid mechanics, each fluid has its physical properties. So, in the case of two liquid phases, some physical properties such as densities, storage coefficients and dynamic viscosities can be different.

Initially, the water molecules are linked to each other by keeping a hydraulic balance. When one introduces an interface between two fluids, this equilibrium is broken and this rupture produces a surface force that ensures the energy balance of molecular cohesion. This force still called capillary force or surface tension is the basis of the CSF method. According to Boliveau et al. [40], this force is transformed to a volume force in the region near the interface via the Delta function:

$\stackrel{\to }{F}=\sigma \mathcal{C}n\delta =\sigma \mathcal{C}\left(\varphi \right)\delta \left(\varphi \right)\nabla \varphi$ (9)

where $\sigma$ is the surface tension coefficient, $n$ the unit normal vector to $\Gamma$ going from ${\Omega }_s$ to ${\Omega }_f$, $\mathcal{C}$ the curvature of $\Gamma$ and $\delta$ the Dirac function. The calculation of the normal and the curvature is done using the following relations:

$\{\begin{array}{l}n=\frac{\nabla \varphi }{\Vert \nabla \varphi \Vert }\\ \mathcal{C}\left(\varphi \right)=\nabla \cdot n=\nabla \cdot \left(\frac{\nabla \varphi }{\Vert \nabla \varphi \Vert }\right)\end{array}$

The flow equation in different phases is given by the relation (7). By integrating this equation on $\Omega$, one has:

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_s}{S}_s}{\rho }_s\frac{\partial H_s}{\partial t}\text{d}\Omega +{\displaystyle {\int }_{{\Omega }_f}{S}_f}{\rho }_f\frac{\partial H_f}{\partial t}\text{d}\Omega +{\displaystyle {\int }_{\partial {\Omega }_s}\left({\rho }_s{v}_{d_s}\right)}\cdot n_s\text{d}s+{\displaystyle {\int }_{\partial {\Omega }_f}({\rho }_s{v}_{d_f})}\cdot n_f\text{d}s\\ ={\displaystyle {\int }_{\Omega }{\rho }_{ϵ}}\left(\varphi \right)Q\text{d}\Omega \end{array}$

Taking into account that $n=n_s=-n_f$, one has:

$\begin{array}{l}{\displaystyle {\int }_{\Omega }{S}_{ϵ}}\left(\varphi \right){\rho }_{ϵ}\left(\varphi \right)\frac{\partial H}{\partial t}\text{d}\Omega +{\displaystyle {\int }_{\partial \Omega }\left({\rho }_{ϵ}\left(\varphi \right)v_d\right)\cdot n\text{d}s}+{\displaystyle {\int }_{\Gamma }\left({\rho }_s{v}_{d_s}\right)}\cdot n\text{d}s-{\displaystyle {\int }_{\Gamma }\left({\rho }_f{v}_{d_f}\right)}\cdot n\text{d}s\\ ={\displaystyle {\int }_{\Omega }{\rho }_{ϵ}}\left(\varphi \right)Q\text{d}\Omega \end{array}$

For $Q=0$ (no extraction, no infiltration to simplify), one obtains:

$\begin{array}{l}{\displaystyle {\int }_{\Omega }{S}_{ϵ}}\left(\varphi \right){\rho }_{ϵ}\left(\varphi \right)\frac{\partial H}{\partial t}\text{d}\Omega +{\displaystyle {\int }_{\Omega }\nabla }\cdot \left({\rho }_{ϵ}\left(\varphi \right)v_d\right)\text{d}\Omega \\ +{\displaystyle {\int }_{\Gamma }\left({\rho }_s{v}_{d_s}\right)}\cdot n\text{d}s-{\displaystyle {\int }_{\Gamma }\left({\rho }_f{v}_{d_f}\right)}\cdot n\text{d}s=0.\end{array}$ (10)

The expression of the Darcy velocity according to the pressure at the interface gives:

$\begin{array}{l}{\displaystyle {\int }_{\Gamma }\left({\rho }_s{v}_s\right)}\cdot n\text{d}s-{\displaystyle {\int }_{\Gamma }({\rho }_f{v}_f)}\cdot n\text{d}s\\ =-{\displaystyle {\int }_{\Gamma }\left[{\rho }_{s}K\nabla \left(\frac{P_s}{{\rho }_{s}g}+z\right)\right]}\cdot n\text{d}s+{\displaystyle {\int }_{\Gamma }\left[{\rho }_{f}K\nabla \left(\frac{P_f}{{\rho }_{f}g}+z\right)\right]}\cdot n\text{d}s\\ =-{\displaystyle {\int }_{\Gamma }\left(\frac{K}{g}\nabla \left(P_s\right)+\frac{K}{g}\rho \left(\varphi \right)\nabla z\right)}\cdot n\text{d}s+{\displaystyle {\int }_{\Gamma }\left(\frac{K}{g}\nabla \left(P_f\right)+\frac{K}{g}\rho \left(\varphi \right)\nabla z\right)}\cdot n\text{d}s\\ =-{\displaystyle {\int }_{\Gamma }\left(\frac{K}{g}\nabla \left(P_s-P_f\right)\right)}\cdot n\text{d}s\end{array}$

According to Laplace’s law, the discontinuity of the pressure at the interface is taken into account by the relation:

$P_s-P_f=\sigma \mathcal{C}$

and according to the approximation of the Ghost Fluid Method (GFM) applied to the CSF method [41], one has:

$\nabla P_s-\nabla P_f=\sigma \mathcal{C}\nabla \mathcal{H}=\sigma \mathcal{C}\delta \nabla \varphi .$

Putting this expression in the previous relation, one obtains the expression of the velocity at the interface:

$\begin{array}{l}{\displaystyle {\int }_{\Gamma }\left({\rho }_s{v}_s\right)}\cdot n\text{d}s-{\displaystyle {\int }_{\Gamma }({\rho }_f{v}_f)}\cdot n\text{d}s\\ =-{\displaystyle {\int }_{\Gamma }\left(\frac{K\sigma \mathcal{C}\left(\varphi \right){\delta }_{ϵ}\left(\varphi \right)\nabla \varphi }{g}\right)}\cdot n\text{d}s=-{\displaystyle {\int }_{\Omega }\nabla }\cdot \left(\frac{K\sigma \mathcal{C}\left(\varphi \right){\delta }_{ϵ}\left(\varphi \right)\nabla \varphi }{g}\right)\text{d}\Omega \end{array}$ (11)

Equations (10) and (11) give:

${\displaystyle {\int }_{\Omega }{S}_{ϵ}}\left(\varphi \right){\rho }_{ϵ}\left(\varphi \right)\frac{\partial H}{\partial t}\text{d}\Omega +{\displaystyle {\int }_{\Omega }\nabla }\cdot \left(\rho \left(\varphi \right)v_d\right)\text{d}\Omega -{\displaystyle {\int }_{\Omega }\nabla }\cdot \left(\frac{K\sigma \mathcal{C}\left(\varphi \right){\delta }_{ϵ}\left(\varphi \right)\nabla \varphi }{g}\right)\text{d}\Omega =0$

Thus, the equation modeling the flow of saline intrusion in a porous medium is

$S_{ϵ}\left(\varphi \right){\rho }_{ϵ}\left(\varphi \right)\frac{\partial H}{\partial t}=\nabla \cdot \left({\rho }_{ϵ}\left(\varphi \right)K\nabla H+\frac{K\sigma \mathcal{C}\left(\varphi \right){\delta }_{ϵ}\left(\varphi \right)\nabla \varphi }{g}\right)$ (12)

The physical properties ${\rho }_{ϵ}\left(\varphi \right)$ and $S_{ϵ}\left(\varphi \right)$ are expressed using the level set function $\varphi$, following the equations:

$\{\begin{array}{l}{\rho }_{ϵ}\left(\varphi \right)={\rho }_f+\left({\rho }_s-{\rho }_f\right){\mathcal{H}}_{ϵ}\left(\varphi \right)\\ S_{ϵ}\left(\varphi \right)=S_f+\left(S_s-S_f\right){\mathcal{H}}_{ϵ}\left(\varphi \right)\end{array}$ (13)

where the smoothed Heaviside function is defined by:

${\mathcal{H}}_{ϵ}\left(\varphi \right)=\{\begin{array}{l}0\varphi <-ϵ\hfill \\ 1/2\left(1+\frac{\varphi }{ϵ}+\frac{1}{\pi }\sin \left(\frac{\pi \varphi }{ϵ}\right)\right)\left|\varphi \right|\le ϵ\hfill \\ 1\varphi >ϵ\hfill \end{array}$ (14)

and the Dirac function by:

${\delta }_{ϵ}\left(\varphi \right)=\frac{\text{d}{\mathcal{H}}_{ϵ}\left(\varphi \right)}{\text{d}\varphi }=\{\begin{array}{l}0\left|\varphi \right|>ϵ\\ \frac{1}{2ϵ}\left(1+\frac{1}{\pi }\cos \left(\frac{\pi \varphi }{ϵ}\right)\right)\left|\varphi \right|\le ϵ\end{array}$ (15)

In these equation $ϵ$ is a smoothing parameter that defines the fictitious thickness of the interface. This thickness must be chosen by respecting the following expression [42]:

$ϵ=\beta \times \Delta x\text{with}1\le \beta \le 2.$

3. Numerical Resolution

The equation of transport (8) is of the hyperbolic type and delicate to solve. Indeed, its approximation by some numerical methods often presents oscillations and numerical scattering. In addition, the problem of loss of mass which knows the level set method, can lead to instability, and therefore to inaccuracy on the exact position of the interface owing to these difficulties, the use of numerical schemes sufficiently robust and accurate is essential to solving the transport Equation (8). Thus, in this work are adopted the five order Weighted Essentially Non-Oscillatory (WENO) scheme for the spatial discretization as well as the third-order Total Variation Diminishing (TVD) Runge-Kutta (RK) scheme for time discretization, has proved their worth in numerous works [27] [28] [30] [43], etc.

The WENO scheme is a weighted version of the Essentially Non-Oscillatory (ENO) scheme. Introduced by Harten et al. [44], the ENO scheme is based on a polynomial interpolation to calculate the fluxes at the boundaries of the cells of space discretization, where the polynomial is the most regular [45]. This allows to obtain a high order of accuracy boundaries even in the vicinity of discontinuities and avoids oscillations. The construction of the polynomials requires a set of neighbouring points called sub-stencils. The number of points used determines the order of the scheme used and the sub-stencils constitute the stencil of the scheme [42].

The WENO scheme uses the same idea, but allocated weight to each point according to the degree of regularity of the solution. The difference between ENO and WENO schemes is that the ENO scheme implies the choice of stencils where the discrete solution is as regular as possible, while the WENO scheme uses a linear combination of stencils by weighting them according to the regularity of the solution. This increases the order of convergence of the method [42].

One considers a mesh of a rectangle $\Omega =\left[a;b\right]\times \left[c;d\right]$ and ${\Omega }_{ik}$ a family of rectangles, such as:

${\Omega }_{ik}=\left[x_{i-1/2};x_{i+1/2}\right]\times \left[z_{k-1/2};z_{k+1/2}\right]$

with $1\le i\le N_x+1$, $1\le k\le N_z$ and ${\Omega }_{ik}$ the cells or controls volumes. One notices:

$\left|{\Omega }_i\right|=\Delta x_i=x_{i+1/2}-x_{i-1/2}$ for $i=1,2,\cdots ,N_x+1$, and $x_i=a+\left(i-1\right)\Delta x$

$\left|{\Omega }_k\right|=\Delta z_k=z_{k+1/2}-z_{k-1/2}$ for $k=1,2,\cdots ,N_z$, and $z_k=c+\left(k-1\right)\Delta z+\Delta z/2$

with $x_i$ and $z_k$ the coordinates of the centers of the cells ${\Omega }_{ik}$.

$\Delta x=\text{max}\left\{\Delta x_i,i=1,\cdots ,N_x+1\right\}$ and $\Delta z=\text{max}\left\{\Delta z_k,k=1,\cdots ,N_z\right\}$

It is assumed that space steps $\Delta x$ and $\Delta z$ are constants.

3.1. Discretization of the Level Set Equation

Spatial discretization: the discretization of Equation (8) using the finite volume method is given:

$\frac{\partial {\varphi }_{ik}}{\partial t}\left(t\right)=-\left(\frac{F_{i+1/2,k}\left(t\right)-F_{i-1/2,k}\left(t\right)}{\Delta x}+\frac{F_{i,k+1/2}\left(t\right)-F_{i,k-1/2}\left(t\right)}{\Delta z}\right)$ (16)

where $F_{i+1/2,k}$ is the approximation of the numerical flux at the point $\left(x_{i+1/2};z_k\right)$. The choice to approximate numerical flux as a function of discrete unknowns determines the numerical scheme. In this study, one uses the scheme WENO5. Since the flux is $F=u\varphi$, then using the upwind scheme one gets:

$\{\begin{array}{l}{F}_{i+1/2,k}\left(t\right)\approx u_{i+1/2,k}{\varphi }_{i+1/2,k}\left(t\right),\text{and}\{\begin{array}{l}\text{if}{u}_{i+1/2,k}>0,\text{wetake}{\varphi }_{i+1/2,k}\left(t\right)={\varphi }_{i+1/2,k}^{-}\left(t\right)\\ \text{if}{u}_{i+1/2,k}<0,\text{wetake}{\varphi }_{i+1/2,k}\left(t\right)={\varphi }_{i+1/2,k}^{+}\left(t\right)\end{array}\\ F_{i-1/2,k}\left(t\right)\approx u_{i-1/2,k}{\varphi }_{i-1/2,k}\left(t\right),\text{and}\{\begin{array}{l}\text{if}{u}_{i-1/2,k}>0,\text{wetake}{\varphi }_{i-1/2,k}\left(t\right)={\varphi }_{i-1/2,k}^{-}\left(t\right)\\ \text{if}{u}_{i-1/2,k}<0,\text{wetake}{\varphi }_{i-1/2,k}\left(t\right)={\varphi }_{i-1/2,k}^{+}(\; t\; )\end{array}\end{array}$

The calculation of ${\varphi }_{i\pm 1/2,k}^{\pm }$ (respectively ${\varphi }_{i,k\pm 1/2}^{\pm }$ ) using the WENO5 scheme is done as follows:

${\varphi }_{i\pm 1/2,k}^{\pm }\left(t\right)={\displaystyle \underset{j=0}{\overset{r-1}{\sum }}{w}_j^{\left(r\right)}}{\varphi }_{i\pm 1/2,k}^{\left(j\right)\pm }\left(t\right);{\varphi }_{i,k\pm 1/2}^{\pm }\left(t\right)={\displaystyle \underset{j=0}{\overset{r-1}{\sum }}{w}_j^{\left(r\right)}}{\varphi }_{i,k\pm 1/2}^{(j)\pm }\left(t\right)$ (17)

with ${\varphi }_{i\pm 1/2,k}^{\left(h\right)}$ (respectively ${\varphi }_{i,k\pm 1/2}^{\left(j\right)}$ ) the interpolation of the flux on each sub-stencils and r the number of points in sub-stencils, $w_j^{\left(r\right)}$ is the nonlinear weights.

The stencils are defined according to the sign of velocity. For a WENO5 scheme, each sub-stencil contains three points as shown in Figure 2 & Figure 3.

Thus, if $r=3$ represents the number of points in each sub-stencil, we have:

$\{\begin{array}{l}\text{if}{u}_{i+1/2,k}>0\text{and}{u}_{i-1/2,k}>0\{\begin{array}{l}{\varphi }_{i+1/2,k}^{-}={\omega }_o^{\left(3\right)}{\varphi }_{i+1/2,k}^{\left(o\right)-}+{\omega }_1^{\left(3\right)}{\varphi }_{i+1/2,k}^{\left(1\right)-}+{\omega }_2^{\left(3\right)}{\varphi }_{i+1/2,k}^{\left(2\right)-}\\ {\varphi }_{i-1/2,k}^{-}={\omega }_o^{\left(3\right)}{\varphi }_{i-1/2,k}^{\left(o\right)-}+{\omega }_1^{\left(3\right)}{\varphi }_{i-1/2,k}^{\left(1\right)-}+{\omega }_2^{\left(3\right)}{\varphi }_{i-1/2,k}^{\left(2\right)-}\end{array}\\ \text{if}{u}_{i+1/2,k}<0\text{and}{u}_{i-1/2,k}<0\{\begin{array}{l}{\varphi }_{i+1/2,k}^{+}={\omega }_o^{\left(3\right)}{\varphi }_{i+1/2,k}^{\left(o\right)+}+{\omega }_1^{\left(3\right)}{\varphi }_{i+1/2,k}^{\left(1\right)+}+{\omega }_2^{\left(3\right)}{\varphi }_{i+1/2,k}^{\left(2\right)+}\\ {\varphi }_{i-1/2,k}^{+}={\omega }_o^{\left(3\right)}{\varphi }_{i-1/2,k}^{\left(o\right)+}+{\omega }_1^{\left(3\right)}{\varphi }_{i-1/2,k}^{\left(1\right)+}+{\omega }_2^{\left(3\right)}{\varphi }_{i-1/2,k}^{\left(2\right)+}\end{array}\end{array}$ (18)

$\{\begin{array}{l}\text{if}{u}_{i,k+1/2}>0\text{and}{u}_{i,k-1/2}>0\{\begin{array}{l}{\varphi }_{i,k+1/2}^{-}={\omega }_o^{\left(3\right)}{\varphi }_{i,k+1/2}^{\left(o\right)-}+{\omega }_1^{\left(3\right)}{\varphi }_{i,k+1/2}^{\left(1\right)-}+{\omega }_2^{\left(3\right)}{\varphi }_{i,k+1/2}^{\left(2\right)-}\\ {\varphi }_{i,k-1/2}^{-}={\omega }_o^{\left(3\right)}{\varphi }_{i,k-1/2}^{\left(o\right)-}+{\omega }_1^{\left(3\right)}{\varphi }_{i,k-1/2}^{\left(1\right)-}+{\omega }_2^{\left(3\right)}{\varphi }_{i,k-1/2}^{\left(2\right)-}\end{array}\\ \text{if}{u}_{i,k+1/2}<0\text{and}{u}_{i,k-1/2}<0\{\begin{array}{l}{\varphi }_{i,k+1/2}^{+}={\omega }_o^{\left(3\right)}{\varphi }_{i,k+1/2}^{\left(o\right)+}+{\omega }_1^{\left(3\right)}{\varphi }_{i,k+1/2}^{\left(1\right)+}+{\omega }_2^{\left(3\right)}{\varphi }_{i,k+1/2}^{\left(2\right)+}\\ {\varphi }_{i,k-1/2}^{+}={\omega }_o^{\left(3\right)}{\varphi }_{i,k-1/2}^{\left(o\right)+}+{\omega }_1^{\left(3\right)}{\varphi }_{i,k-1/2}^{\left(1\right)+}+{\omega }_2^{\left(3\right)}{\varphi }_{i,k-1/2}^{\left(2\right)+}\end{array}\end{array}$

where the nonlinear weights $w_j^{\left(r\right)}$ are defined as follow:

${\omega }_j^{\left(r\right)}=\frac{{\tilde{\omega }}_j}{{\displaystyle \underset{l=0}{\overset{r-1}{\sum }}{\tilde{\omega }}_l}},\text{for}j=0,\cdots ,r-1and{\tilde{\omega }}_j=\frac{{\alpha }_j^{\left(r\right)}}{{\left(ϵ+{\beta }_j^{\pm }\right)}^2}$ (19)

where ${\alpha }_j^{\left(r\right)}$ are the linear weights, ${\beta }_j$ the smoothness indicators used for the weight calculations and $ϵ>0$ the small constant allowing to avoid that the denominator of ${\tilde{\omega }}_j$ let be equal zero.

The parameters ${\alpha }_j^{\left(r\right)}$ and ${\beta }_j^{\pm }$ are defined below (Equation (24) [42] ). Consequently, one has:

Figure 2. Stencil for ${\varphi }_{i+1/2,k}$ : presentation of the sub-stencils $S_i^{\pm }$ making up the stencils.

Figure 3. Stencil for ${\varphi }_{i-1/2,k}$ : presentation of the sub-stencils $S_i^{\pm }$ making up the stencils.

${\omega }_o^{\left(3\right)}=\frac{{\tilde{\omega }}_o}{{\tilde{\omega }}_o+{\tilde{\omega }}_1+{\tilde{\omega }}_2};{\omega }_1^{\left(3\right)}=\frac{{\tilde{\omega }}_1}{{\tilde{\omega }}_o+{\tilde{\omega }}_1+{\tilde{\omega }}_2};{\omega }_2^{\left(3\right)}=\frac{{\tilde{\omega }}_2}{{\tilde{\omega }}_o+{\tilde{\omega }}_1+{\tilde{\omega }}_2}$

with

${\tilde{\omega }}_o=\frac{{\alpha }_o^{\left(3\right)}}{{\left(ϵ+{\beta }_o^{\pm }\right)}^2};{\tilde{\omega }}_1=\frac{{\alpha }_1^{\left(3\right)}}{{\left(ϵ+{\beta }_1^{\pm }\right)}^2};{\tilde{\omega }}_2=\frac{{\alpha }_2^{\left(3\right)}}{{\left(ϵ+{\beta }_2^{\pm }\right)}^2}.$

The linear weights of the flux ${\alpha }_o^{\left(3\right)}$, ${\alpha }_1^{\left(3\right)}$. And ${\alpha }_2^{\left(3\right)}$ are weighted so that ${\alpha }_o^{\left(3\right)}+{\alpha }_1^{\left(3\right)}+{\alpha }_2^{\left(3\right)}=1$. For this, one takes ${\alpha }_o^{\left(3\right)}=0.3$, ${\alpha }_1^{\left(3\right)}=0.6$ and ${\alpha }_2^{\left(3\right)}=0.1$ [42].

The success of WENO schemes lies in the choice of standardized weights ( ${\omega }_j^{\left(r\right)}$ ). For consistency and stability of the scheme, these coefficients must satisfy the relations:

${\omega }_j^{\left(r\right)}\ge 0,{\displaystyle \underset{j=0}{\overset{r-1}{\sum }}{\omega }_j^{\left(r\right)}}=1\text{and}ϵ={10}^{-6}\max \left\{v_j^2\right\}+{10}^{-99}.$ (20)

where $v_j$ represents the mean value at the points of indices $\left\{\left(i-2,k\right);\left(i-1,k\right);\left(i,k\right);\left(i+1,k\right);\left(i+2,k\right)\right\}$ [28].

Now the flux ${\varphi }_{i+1/2,k}^{\left(j\right)\pm }$ for each sub-stencil are obtained as follows:

$\{\begin{array}{l}{\varphi }_{i+1/2,k}^{\left(o\right)\pm }\left(t\right)=\frac{1}{3}{v}_1^{\pm }-\frac{7}{6}{v}_2^{\pm }+\frac{11}{6}{v}_3^{\pm }\text{for the sub-stencil 1}\\ {\varphi }_{i+1/2,k}^{\left(1\right)\pm }\left(t\right)=-\frac{1}{6}{v}_2^{\pm }+\frac{5}{6}{v}_3^{\pm }+\frac{1}{3}{v}_4^{\pm }\text{for the sub-stencil 2}\\ {\varphi }_{i+1/2,k}^{\left(2\right)\pm }\left(t\right)=\frac{1}{3}{v}_3^{\pm }+\frac{5}{6}{v}_4^{\pm }-\frac{1}{6}{v}_5^{\pm }\text{for the sub-stencil 3}\end{array}$ (21)

The expressions $v_j^{\pm }$ are calculated based on the sub-stencils. $v_1^{\pm }$, $v_2^{\pm }$, $v_3^{\pm }$, $v_4^{\pm }$ and $v_5^{\pm }$, are the mean values of $\varphi$ at points $\left\{\left(i-2,k\right);\left(i-1,k\right);\left(i,k\right);\left(i+1,k\right);\left(i+2,k\right)\text{and}\left(i+3,k\right)\right\}$.

For sub-stencils associated with ${\varphi }_{i+1/2,k}^{-}$, i.e. $S_1^{-}$, $S_2^{-}$ and $S_3^{-}$, one has:

$\begin{array}{c}{v}_1^{-}={\varphi }_{i-2,k},v_2^{-}={\varphi }_{i-1,k},v_3^{-}={\varphi }_{ik},v_4^{-}={\varphi }_{i+1,k},v_5^{-}={\varphi }_{i+2,k}\end{array}$ (22)

For sub-stencils associated with ${\varphi }_{i+1/2,k}^{+}$, i.e. $S_1^{+}$, $S_2^{+}$ and $S_3^{+}$, one has:

$\begin{array}{c}{v}_1^{+}={\varphi }_{i+3,k},v_2^{+}={\varphi }_{i+2,k},v_3^{+}={\varphi }_{i+1,k},v_4^{+}={\varphi }_{ik},v_5^{+}={\varphi }_{i-1,k}\end{array}$ (23)

The expressions ${\beta }_j$ represent the parameters allowing to measure the regularity of the function on each sub-stencil. They are given by:

$\{\begin{array}{l}{\beta }_o^{\pm }=\frac{13}{12}{\left(v_1^{\pm }-2v_2^{\pm }+v_3^{\pm }\right)}^2+\frac{1}{4}{\left(v_1^{\pm }-4v_2^{\pm }+3v_3^{\pm }\right)}^2\\ {\beta }_1^{\pm }=\frac{13}{12}{\left(v_2^{\pm }-2v_3^{\pm }+v_4^{\pm }\right)}^2+\frac{1}{4}{\left(v_2^{\pm }-v_4^{\pm }\right)}^2\\ {\beta }_2^{\pm }=\frac{13}{12}{\left(v_3^{\pm }-2v_4^{\pm }+v_5^{\pm }\right)}^2+\frac{1}{4}{\left(3v_3^{\pm }-4v_4^{\pm }+v_5^{\pm }\right)}^2\end{array}$ (24)

To obtain the flux ${\varphi }_{i,k\pm 1/2}\left(t\right)$ on the stencil $\left\{\left(i,k-2\right);\left(i,k-1\right);\left(i,k\right);\left(i,k+1\right);\left(i,k+2\right)\right\}$, one uses the same process that for the flux ${\varphi }_{i\pm 1/2,k}\left(t\right)$.

Time discretization: Let L be a discrete operator, $\Delta t$ the time step. The time discretization of the transport Equation (17) using the explicit Euler scheme gives:

${\varphi }_{ik}^{n+1}={\varphi }_{ik}^n+\Delta tL(\; \varphi \; n\; )$

with

$L\left({\varphi }^n\right)=-\left(\frac{F_{i+1/2,k}^n-F_{i-1/2,k}^n}{\Delta x}+\frac{F_{i,k+1/2}^n-F_{i,k-1/2}^n}{\Delta z}\right)$

Thus, the third-order TVD Runge-Kutta scheme, is:

$\{\begin{array}{l}{\varphi }^{\left(1\right)}={\varphi }^n+\Delta tL\left({\varphi }^n\right)\\ {\varphi }^{\left(2\right)}=\frac{3}{4}{\varphi }^n+\frac{1}{4}{\varphi }^{\left(1\right)}+\frac{1}{4}\Delta tL\left({\varphi }^1\right)\\ {\varphi }^{\left(n+1\right)}=\frac{1}{3}{\varphi }^n+\frac{2}{3}{\varphi }^{\left(2\right)}+\frac{2}{3}\Delta tL\left({\varphi }^2\right)\end{array}$ (25)

3.2. Discretization of the Flow Equation

The flow Equation (12) is of the parabolic type. One uses the method of volumes finite in space and the explicit Euler scheme in time for its resolution. Therefore, by integrating it on ${\Omega }_{ik}\times \left[t_n,t_{n+1}\right]$ for all $1\le i\le N_x+1$ and $1\le k\le N_z$ an $n\ge 0$, one gets:

$\begin{array}{c}{H}_{ik}^{n+1}=H_{ik}^n+\frac{\Delta t\left(\frac{2{\rho }_{ϵ}\left({\varphi }_{ik}^n\right){\rho }_{ϵ}\left({\varphi }_{i+1,k}^n\right)}{{\rho }_{ϵ}\left({\varphi }_{ik}^n\right)+{\rho }_{ϵ}\left({\varphi }_{i+1,k}^n\right)}\right)}{\Delta x^2{S}_{ϵ}\left({\varphi }_{ik}^n\right){\rho }_{ϵ}\left({\varphi }_{ik}^n\right)}\left[K\left(H_{i+1,k}^n-H_{i,k}^n\right)+\sigma \left(\frac{\mathcal{C}\left({\varphi }_{i+1,k}^n\right)+\mathcal{C}\left({\varphi }_{i,k}^n\right)}{2}\right)\left(\frac{{\delta }_{ϵ}\left({\varphi }_{i+1,k}^n\right)+{\delta }_{ϵ}\left({\varphi }_{i,k}^n\right)}{2}\right)\left({\varphi }_{i+1,k}^n-{\varphi }_{i,k}^n\right)\right]\\ -\frac{\Delta t\left(\frac{2{\rho }_{ϵ}\left({\varphi }_{ik}^n\right){\rho }_{ϵ}\left({\varphi }_{i-1,k}^n\right)}{{\rho }_{ϵ}\left({\varphi }_{ik}^n\right)+{\rho }_{ϵ}\left({\varphi }_{i-1,k}^n\right)}\right)}{\Delta x^2{S}_{ϵ}\left({\varphi }_{ik}^n\right){\rho }_{ϵ}\left({\varphi }_{ik}^n\right)}\left[K\left(H_{i,k}^n-H_{i-1,k}^n\right)+\sigma \left(\frac{\mathcal{C}\left({\varphi }_{i-1,k}^n\right)+\mathcal{C}\left({\varphi }_{i,k}^n\right)}{2}\right)\left(\frac{{\delta }_{ϵ}\left({\varphi }_{i-1,k}^n\right)+{\delta }_{ϵ}\left({\varphi }_{i,k}^n\right)}{2}\right)\left({\varphi }_{i,k}^n-{\varphi }_{i-1,k}^n\right)\right]\\ +\frac{\Delta t\left(\frac{2{\rho }_{ϵ}\left({\varphi }_{ik}^n\right){\rho }_{ϵ}\left({\varphi }_{i,k+1}^n\right)}{{\rho }_{ϵ}\left({\varphi }_{ik}^n\right)+{\rho }_{ϵ}\left({\varphi }_{i,k+1}^n\right)}\right)}{\Delta z^2{S}_{ϵ}\left({\varphi }_{ik}^n\right){\rho }_{ϵ}\left({\varphi }_{ik}^n\right)}\left[K\left(H_{i,k+1}^n-H_{i,k}^n\right)+\sigma \left(\frac{\mathcal{C}\left({\varphi }_{i,k+1}^n\right)+\mathcal{C}\left({\varphi }_{i,k}^n\right)}{2}\right)\left(\frac{{\delta }_{ϵ}({\varphi }_{i,k+1}^n)+{\delta }_{ϵ}\left({\varphi }_{i,k}^n\right)}{2}\right)\left({\varphi }_{i,k+1}^n-{\varphi }_{i,k}^n\right)\right]\end{array}$

$-\frac{\Delta t\left(\frac{2{\rho }_{ϵ}\left({\varphi }_{ik}^n\right){\rho }_{ϵ}\left({\varphi }_{i,k-1}^n\right)}{{\rho }_{ϵ}\left({\varphi }_{ik}^n\right)+{\rho }_{ϵ}\left({\varphi }_{i,k-1}^n\right)}\right)}{\Delta z^2{S}_{ϵ}\left({\varphi }_{ik}^n\right){\rho }_{ϵ}\left({\varphi }_{ik}^n\right)}\left[K\left(H_{i,k}^n-H_{i,k-1}^n\right)+\sigma \left(\frac{\mathcal{C}\left({\varphi }_{i,k-1}^n\right)+\mathcal{C}\left({\varphi }_{i,k}^n\right)}{2}\right)\left(\frac{{\delta }_{ϵ}\left({\varphi }_{i,k-1}^n\right)+{\delta }_{ϵ}\left({\varphi }_{i,k}^n\right)}{2}\right)\left({\varphi }_{i,k}^n-{\varphi }_{i,k-1}^n\right)\right]$ (26)

Thus, the algorithm used for the resolution of the saline intrusion problem is as follows:

a) initialization of the interface using the function level set $\varphi$ ;

b) for each time step:

● calculation of the geometric properties of the interface such as the normal vector $n$ and the curvature $\mathcal{C}$, as well as the Heaviside function ${\mathcal{H}}_{ϵ}$, the Dirac function ${\delta }_{ϵ}$, the storage coefficient $S_{ϵ}$ and the density ${\rho }_{ϵ}$ ;

● resolution of the flow equation;

● calculation of the velocity knowing that $U=v_d/\Phi$, with $v_d$ the Darcy velocity and $\Phi$ the porosity of the domain $\Omega$ ;

● resolution of the level set equation, using the velocity $U$ obtained in step 3;

● application of the reinitialization algorithm to the level set function $\varphi$ if the distance property is not respected i.e. if $\Vert \nabla \varphi \left(x,t\right)\Vert \ne 1$ ;

● update of the level set function $\varphi$ and of the hydraulic head H;

c) end if the initialized time is equal to the final time i.e. $t=t_f$, otherwise return to step b).

4. Numerical Results: Application to Henry Problem

In order to test the efficiency of the method developed, one will use it on Henry’s problem [46] which models the saline intrusion in a confined aquifer using the density-dependent groundwater flow approach.

Henry’s problem has been taken up by many authors, who have solved it either by the quasi-analytical method, either by numerical methods withs finite difference, finite volume or finite element schemes (Segol et al. [32], Frind [47], Huyakorn et al. [48], Voss and Souza [49], Croucher and d’Osullivan [50] ).

This problem is one of the most widely used to validate saline intrusion models in development. Henry’s problem is illustrated in Figure 4. The parameters and boundary conditions used are presented in Table 1.

Figure 4. Illustration of Henry’s problem.

Table 1. Parameters used for Henry’s problem.

The interface is represented by the zero level curve of the level set function $\varphi$, and Figure 4 allows to adopt the following initialization:

${\varphi }_o\left(x,z\right)=x-200.$ (27)

At initial condition, the interface is represented by the vertical axis $x=200$. So, the interface between freshwater and saltwater is given by:

$\varphi \left(x=200,z,t=0\right)={\varphi }_o\left(x=200,z\right)=0.$

For the boundary conditions of the level set function, one adopts the choice of homogeneous Neumann conditions at bounds $z=0$ and $z=100$, i.e.

$\frac{\partial \varphi }{\partial n}\left(x,0,t\right)=\frac{\partial \varphi }{\partial n}\left(x,100,t\right)=0.$

Table 2 gives the values of the surface tension coefficient.

One notices that the results obtained by Lee and Cheng [52], and Henry [46] do not correspond with any other results obtained in other works such as Pinder and Cooper, Segol et al. [32], Frind [47], etc. In addition, the parameters used by the last two studies have been the subject of several discussions. However, Croucher and O’sullivan [50], and Voss and Souza [49] explain that the dimensional parameter noted $b=\left(\Phi D\right)/U$ influences the position of the interface between freshwater and saltwater. Taking into account the porosity in the advection-dispersive equation, one gets the value of $b=0.1$ used by Lee and Cheng, as well as Henry. Some studies did not take into account this porosity and use the value $b=0.035$. Therefore, the solution obtained by Lee and Cheng, and Henry (with $b=0.1$ ) cannot be compared to the solution obtained by Pinder and Cooper, Segol et al., Frind and Huyakorn (with $b=0.035$ ).

The influence of the different parameters used for the Henry problem has been the subject of numerous studies. These have shown that modification of parameters such as freshwater flow, dispersion coefficient, the coupling of the densities can influence the progression of the interface and the toe [49] [50] [51]. In this work one calls standard solution, the result obtained by Henry (1964), and Henry’s modified solution, the result obtained by Simpson and Clement by

Table 2. Parameters used for Henry’s problem.

modifying the parameters such as dimensional parameters a and b influenced by flux U and dispersion coefficient D [53]. In addition, Simpson and Clement made a comparative study between the coupled solution (assuming that the density depends on the concentration of salt in each fluid) and the uncoupled solution (ignoring density variability within the domain). It emerged that the coupled solution obtained by Simpson and Clement coincides exactly with Henry’s solution (1964), and will be called Henry’s solution in this work.

Several techniques have been proposed for the improvements of the level set method such as the Accurate Conservative Level Set (ACLS) method, the coupling of the different tracking interface methods (VOF-LS), the use of reinitialization algorithm, etc.

In this study, the reinitialization algorithm is used to improve the results. In addition, two numerical schemes ENO1 and ENO2 are used to resolve the reinitialization.

One presents in the next the numerical results obtained by using two tests relating to the standard Henry problem and Henry modified problem.

● Test 1: The test aims to compare the solution of the level set method with reinitialization (LSMR), with Henry’s standard solution and Simpson and Clement’s uncoupled solution, using the parameters given in Table 1. The LSMR represents a sharp interface approach associated with the level set method with reinitialization algorithm using the ENO1 numerical scheme for spatial discretization. The reinitialization algorithm can use a given number it of iterations for the fictitious time. In this case, one uses $it=100$.

Figure 5(a) shows the permanent solutions of the level set function $\varphi$ and Figure 5(b) the hydraulic head obtained after 6000 days. One observes in Figure 5(b) a slight modification of the curve of level 2 of the hydraulic head compared to that obtained by Huyakorn [48] and Hamidi et al. [33] in the variation density approach. But this corresponds to the results obtained by Abarca [54] on the Henry problem in the case of the uncoupled model with a variable coefficient of dispersion (no diffusion), and in the case of constant dispersion coefficient (i.e. pure diffusion).

Figure 6(a) proposes a comparison of a semi-analytical solution obtained by Henry [46], an uncoupled problem solution obtained by Simpson and Clement [53] and the solution obtained in this study with LSMR.

One sees that LSMR reproduces fairly well the toe or salt bevel (intersection of the interface with the substrant i.e. the bound $z=0$ ) and the tip (intersection of the interface with the bound $z=100\text{m}$ ) obtained by Henry in the case of variable dispersion coefficient. The interface obtained with the LSMR method possesses a

(a)(b)

Figure 5. Level set function $\varphi$ and the hydraulic head H at time $t=6000$ days obtained using a reinitialization with the ENO1 scheme (number of iterations $it=100$ ). (a) Level Set function. (b) Hydraulic head.

concavity that moves it away from the interface obtained by Henry which presents an inflection point [54]. This distance is justified by the absence of the dispersion coefficient [54]. On the other hand, the LSMR solution obtained reproduces quite well the behavior of the interface in the case of the uncoupled solution of Simpson and Clement with a slight spacing. This result is undoubtedly due to the number of iterations used in the reinitialization algorithm which affects the narrowing of the contour lines. This solution also corresponds to Henry’s solution presented by Abarca in the variable density case [54].

(a)(b)

Figure 6. Comparison of the interfaces obtained with LSMR with those obtained by Henry [46] and Simpson and Clement [53] at time $t=6000$ days. (a) Result obtained by the level set method; (b) Interfaces obtained with the ENO1 scheme for $it=100$ and $it=200$.

Figure 6(b) takes back the experiment made in Figure 6(a), but with a modification of the number it of iterations in the reinitialization algorithm. Here, the interface obtained with the level set method for $it=200$ and the uncoupled solution of Simpson and Clement coincides exactly at their intersections with the bounds of the domain (axes $z=0$ and $z=100\text{m}$ ) and over the interval $z\in \left[0,60\right]$. Elsewhere one still observes a gap of the two interfaces as in Figure 6(a). This gap can be justified by the absence of the coefficient of dispersion

Figure 7. Comparison of the LMSR ENO2 solution to the modified solutions obtained by Henry and Simpson and Clement.

allowing to obtain the inflection point. [54] such solutions have been obtained by the softwares Sharp and Sutra [54], etc.

Note that in the Henry problem and Simpson and Clement case, the salt concentration is represented by $c=0.5$ (0.5 isochlor distribution).

● Test 2: This second test consists in comparing the LSMR ENO2 solution to the modified solution of Henry and Simpson and Clement (in uncoupled case). The modified solution of Henry and Simpson and Clement is obtained by replacing the value of freshwater flow $U=6.6\times {10}^{-3}\text{m}/\text{d}$ at the left boundary side of the domain $\Omega$ by $U=3.3\times {10}^{-3}\text{m}/\text{d}$. The aim is to show the sensitivity and the influence of parameters on the position of the salt bevel (toe) and the behavior of the interface. According to Figure 6(b), one notices that the number of iterations in the reinitialization algorithm affects the position of the interface between freshwater and saltwater. Thus, this test aims to study the influence of numerical scheme in this algorithm, by replacing the numerical scheme ENO1 by the scheme ENO2 with a number of iterations $it=100$.

Figure 7 shows that the LSMR solution reproduces quite the tip (salt wedge) of Henry’s modified solution. It represents quite the point of inflection and the behavior of the interface in the case of the modified solution by Simpson and Clement (uncoupled case). Elsewhere one observes the same behaviour as in Figure 6(a).

5. Conclusions

This work consisted of solving a two-dimensional sharp interface saline water intrusion by using level set functions. The level set method has been used in many two-phase or multiphasic problems; it makes it possible to track the evolution of the different phases and to capture the interface separating them more particularly when these interfaces are represented by closed curves.

The LSMR method developed in this work does not take into account the dispersion of the salt and the closeness of the interface. One notices that the results obtained are similar to those obtained by Abarca et al. with constant dispersion coefficient (otherwise called Henry’s solution with molecular diffusion) [54], and by Simpson and Clement in the uncoupled case. Thus, the results obtained clearly show that the method proposed can track and capture the position of tip and toe with accuracy, and the behaviour and evolution of the freshwater/saltwater interface.

Acknowledgements

The authors would like to thank African Union Scholarship Scheme to sponsor the PhD studies of the first of the them under grant Mwalimu Nyerere (Ref: HRST/ED/606). They also thank Mr. Yawo Kobara for his remarks and observations.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Bear, J., Cheng, H., Sorek, S., Ouazar, D. and Herrera (1999) Seawater Intrusion in Coastal Aquifers Concepts, Methods and Practices. Kluwer Academic Publishers, Dordrecht. https://doi.org/10.1007/978-94-017-2969-7
[2]	Werner, A., Bakker, M., Post, V., Alexander, V., Lu, C., Behzad, A.A., Craig, S. and Barry, D. (2013) Seawater Intrusion Processes, Investigation and Management: Recent Advances and Future Challenges. Advances in Water Resources, 51, 3-26. https://doi.org/10.1016/j.advwatres.2012.03.004
[3]	Mohamed, B., Margaret, S., Ramachandran, G. and Shelat, K. (2014) Vulnerability of Agriculture, Water and Fisheries to Climate Change.
[4]	Alix, L. (2006) Current Variations in Sea Level: Observations and Climatic Causes. Laboratory of Geophysics and Oceanographic Space, Toulouse Meteorology, 59.
[5]	Solomon, S., Qin, D., et al. (2007) The Scientific Elements. Contribution of Working Group I to Fourth Assessment Report of the Intergovernmental Panel on Climate Change. GIES, Climatic Changes.
[6]	Cheng, A. and Ouazar, D. (2004) Coastal Aquifer Management. Monitoring, Modeling, and Case Studies. https://doi.org/10.1201/9780203493496
[7]	Bear, J. and Cheng, A. (2010) Modeling Groundwater Flow and Contaminant Transport. Part of the Theory and Applications of Transport in Porous Media Book Series (TATP, Volume 23). https://doi.org/10.1007/978-1-4020-6682-5
[8]	Sherif, M. (1999) The Nile Delta Aquifer in Egypt, Chapter 17. In: Seawater Intrusion in Coastal Aquifers: Concepts, Methods and Practices, Book Series: Theory and Application of Transport in Porous Media, Kluwer Academic Publishers, Dordrecht, Vol. 14, 559-590. https://doi.org/10.1007/978-94-017-2969-7_17
[9]	Rifai, M., et al. (1956) Dispersion Phenomena in Laminar Flow through Porous Media. University of California Institute of Engineering, San Diego, Progress Report 2.
[10]	Essaid, H. (1986) A Comparison of the Coupled Fresh Water—Salt Water Flow and the Ghyben-Herzberg Sharp Interface Approaches to Modeling of Transient Behavior in Coastal Aquifer Systems. Journal of Hydrology, 86, 169-193. https://doi.org/10.1016/0022-1694(86)90012-0
[11]	Ogata, A. (1970) Theory of Dispersion in Granular Medium. U.S. Geological Survey Professional Paper 411-I. https://doi.org/10.3133/pp411I
[12]	Hubbert, M. (1940) The Theory of Ground Water Motion. Journal of Geology, 48, 785-944. https://doi.org/10.1086/624930
[13]	Sark, M. (1980) Seawater Intrusion and Potential Yield of Aquifers in the Soquel Aptos Area, Santa Cruz Country, California. Geological Survey, Water Resources Investigations, 80-84.
[14]	Badon-Ghyben (1888) Nota in verband met de voorgenomen putboring nabij Amsterdam (Notes on the Probable Results of Well Drilling near Amsterdam). Tijdschrift van het Koninklijk Institut van Ingenieurs, The Hague, 8-22.
[15]	Herzberg, A. (1901) Die Wasserversorgung einiger Nordseebader (The Water Supply of Parts of the North Sea Coast in Germany). Journal Gasbeleuchtung u. Wasserversorgung, 44, 815-819.
[16]	Henry, H. (1959) Salt Intrusion into Fresh-Water Aquifers. Journal of Geophysical Research, 64, 1911-1919. https://doi.org/10.1029/JZ064i011p01911
[17]	Bear, J. and Verruijt, A. (1987) Modeling Groundwater Flow and Pollution. Theory and Applications of Transport in Porous Media. https://doi.org/10.1007/978-94-009-3379-8
[18]	Cooper, H. (1959) A Hypothesis Concerning the Dynamic Balance of Fresh Water and Salt Water in a Coastal Aquifer. Journal Geophysics Resources, 64, 461-467. https://doi.org/10.1029/JZ064i004p00461
[19]	Fetter, C. (1972) Position of the Saline Water Interface beneath Oceanic Islands. Water Resources Research, 8, 1307-1315. https://doi.org/10.1029/WR008i005p01307
[20]	Mualem, Y. and Bear, J. (1974) The Shape of the Interface in Steady Flow in a Stratified Aquifer. Water Resources Research, 10, 1207-1215. https://doi.org/10.1029/WR010i006p01207
[21]	Amirat, Y., Hamdache, K. and Ziani, A. (1996) Mathematical Analysis for Compressible Miscible Displacement Models in Porous Media. Mathematical Models and Methods in Applied Science, 6, 729-747. https://doi.org/10.1142/S0218202596000316
[22]	Choquet, C. (2010) Parabolic and Degenerate Parabolic Models for Pressure-Driven Transport Problems. Mathematical Models and Methods in Applied Science, 20, 543-566. https://doi.org/10.1142/S0218202510004337
[23]	Mory, M. (2015) Mixed Net/Diffuse Interface Approach for Saline Intrusion Problems in Basement: Modeling, Mathematical Analysis and Digital Illustrations. PhD Dissertation, Rochelle University, Rochelle.
[24]	Abudawia, A. (2015) Numerical Analysis of a Finite Element Approximation for an Intrusion Model Saline in Coastal Aquifers. Ph.D. Thesis, Littoral Côte d’Opale University, Dunkerque.
[25]	Shamir, U. and Dagan, G. (1971) Motion If the Seawater Interface in Costal Aquifers: A Numerical Solution. Water Resources Research, 7, 644-657. https://doi.org/10.1029/WR007i003p00644
[26]	Anderson, M. (1976) Unsteady Groundwater Flow beneath Strip Oceanic Islands. Water Resources Research, 20, 640-644. https://doi.org/10.1029/WR012i004p00640
[27]	Sethian, J. (1996) Level Set Methods. Evolving Interface in Geometry, Fluid Mechanics, Computer Vision, and Materials Science. Cambridge Monographs on Applied and Computational Mathematics.
[28]	Osher, S. and Fedkiw, R. (2003) Level Set Methods and Dynamic Implicit Surfaces. Applied Mathematical Sciences. Springer, Berlin. https://doi.org/10.1007/b98879
[29]	Zhu, J. and Sethian, J. (1992) Projection Methods Coupled to Level Set Interface Techniques. Journal of computational Physics, 102, 128-138. https://doi.org/10.1016/S0021-9991(05)80011-7
[30]	Osher, S. and Sethian, J. (1988) Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations. Journal of Computational Physics, 79, 12-49. https://doi.org/10.1016/0021-9991(88)90002-2
[31]	Narkbuakaew, W., Nagahashi, H., Aoki, K. and Kubota, Y. (2014) An Efficient Liver-Segmentation System Based on a Level-Set Method and Consequent Processes. Journal of Biomedical Science and Engineering, 7, 994-1004. https://doi.org/10.4236/jbise.2014.712097
[32]	Segol, G., Pinder, G. and William, G. (1975) A Galerkin-Finite Element Technique for Calculating the Transient Position of the Saltwater Front. Water Resources Research, 11, 343-347. https://doi.org/10.1029/WR011i002p00343
[33]	Hamidi, M., Reza, S. and Yazdi, S. (2006) Numerical Modeling of Seawater Intrusion in Coastal Aquifer Using Finite Volume Unstructured Mesh Method. Proceedings of the 9th WSEAS International Conference on Applied Mathematics, Istanbul, 27-29 May 2006, 572-579.
[34]	Sethian, J. and Smereka, P. (2003) Level Set Methods for Fluid Interfaces. Annual Review of Fluid. Mechanics, 35, 341-372. https://doi.org/10.1146/annurev.fluid.35.101101.161105
[35]	Garzon, M., Sethian, J., et al. (2000) Numerical Simulation of Non-Viscous Liquid Pinch-Off Using a Coupled Level Set Boundary Integral Method. Journal of Computational Physics, 228, 6079-6106. https://doi.org/10.1016/j.jcp.2009.04.048
[36]	Thibaut, M. (2007) Development of a Level Set Method for Interface Monitoring, Application to Liquid Jet Rupture. Ph.D. Thesis, Rouen University, Rouen.
[37]	Sussman, M., Smereka, P. and Osher, S. (1994) A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flows. Journal of Computational Physics, 114, 146-159. https://doi.org/10.1006/jcph.1994.1155
[38]	Brackbill, J., Kothe, D. and Zemach, C. (1992) A Continuum Method for Modeling Surface Tension. Journal of Computational Physics, 100, 335-354. https://doi.org/10.1016/0021-9991(92)90240-Y
[39]	Sefollahi, M. and Shirani, E. (2004) An Improved Method for Calculation of Interface Pressure Force in PLICVOF Methods. Abdus Salam International Centre for Theoretical Physics, Trieste.
[40]	Boliveau, A., Fortin, A. and Demay, Y. (1998) A Two-Dimensional Numerical Method for the Deformation of Drops with Surface Tension. International Journal of Computational Fluid Dynamics, 10, 225-240. https://doi.org/10.1080/10618569808961687
[41]	Popinet, S. (2018) Numerical Models of Surface Tension. Annual Review of Fluid Mechanics, 50, 49-75. https://doi.org/10.1146/annurev-fluid-122316-045034
[42]	Chi-Wang (1997) Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws. Institute for Computer Application in Science and Engineering. NASA Langley Research Center, Hampton, Report, 97-65.
[43]	Sébastien, T. (2004) Development of an Interface Monitoring Method. Flow Applications Diphasic. Ph.D. Thesis, Rouen University, Rouen.
[44]	Harten, A., Engquist, B., Osher, S. and Chakravarthy, S. (1987) Uniformly High-Order Accurate Essentially Non-Oscillatory Schemes III. Hydrology. Journal of Computational Physics, 71, 231-303. https://doi.org/10.1016/0021-9991(87)90031-3
[45]	Biswarup, B. and Kumar, R. (2007) Accuracy Preserving ENO and WENO Schemes Using Novel Smoothness Measurement. Abdus Salam International Center for Theoretical Physics, Trieste.
[46]	Henry, H. (1964) Effects of Dispersion on Salt Encroachment in Coastal Aquifers, Sea Water in Coastal Aquifers. U.S. Geological Survey. Water Supply, 1613-C.
[47]	Frind, E. (1982) Simulation of Long-Term Transient Density-Dependent Transport in Groundwater. Advances in Water Resources, 5, 73-88. https://doi.org/10.1016/0309-1708(82)90049-5
[48]	Huyakorn, P., Andersen, P., Mercer, J. and White, H. (1987) Saltwater Intrusion in Aquifers: Development and Testing of a Three-Dimensional Finite Element Model. Water Resources Research, 23, 293-312. https://doi.org/10.1029/WR023i002p00293
[49]	Voss, C. and Souza, W. (1987) Modeling a Regional Aquifer Containing a Narrow Transition between Freshwater and Saltwater Using Solution Transport Simulation: Part I Theory and Methods. Geological Survey, U.S., 493514.
[50]	Croucher, A. and O’Sullivan, M. (1995) The Henry Problem for Saltwater Intrusion. Water Resources Research, 31, 1809-1814. https://doi.org/10.1029/95WR00431
[51]	Thomas, A., Eldho, T. and Rastogi, A. (2016) Simulation of Seawater Intrusion in Coastal Confined Aquifer Using a Point Collocation Method Based Meshfree Model. Journal of Water Resource and Protection, 8, 534-549. https://doi.org/10.4236/jwarp.2016.84045
[52]	Lee, C. and Cheng, R. (1974) On Seawater Encroachment in Coastal Aquifers. Water Resources Research, 10, 1039-1043. https://doi.org/10.1029/WR010i005p01039
[53]	Simpson, M. and Clement, P. (2004) Improving the Worthiness of the Henry Problem as a Benchmark for Density-Dependent Groundwater Flow Models. Water Resources Research, 40, W01504. https://doi.org/10.1029/2003WR002199
[54]	Abarca, E., Carrera, J., Xavier, S.V. and Dentz, M. (2006) Anisotropic Dispersive Henry Problem. Advances in Water Resources, 30, 913-926. https://doi.org/10.1016/j.advwatres.2006.08.005