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

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 re-lating 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


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.

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]: are respectively the specific storage coefficient, the hydraulic head, the hydraulic conductivity, the source term and the Darcy velocity in the different phases α ( f α = represents the freshwater phase; 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.

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 de-  (Figure 1). In 2D the evolution interface Γ separating the two sub-domains − Ω and + Ω is defined by: The level set function φ is negative in − Ω , positive outside i.e. in + Ω and zero on ( ) For two immiscible fluids, the domain Ω can also be written as the union of three sub-domains, 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: , 0 for where d is the signed distance function.
The initial position 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 Ω domain.
with a given initial condition ( x 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.
( ) The distance property ( ) is important for signed distance functions, it is generally the basis of the level set method. In the case of loss of this , the method then becomes imprecise. This is why it is necessary to impose a condition on the evolution of φ in Ω called a reinitialization or redistanciation algorithm, that is to say ( ) , 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: where τ represents a fictitious reinitialization time of the function φ , d represents the level set function at fictitious time τ , and the ( ) sign ⋅ function represents the sign of the level set function φ . This function is smoothed around the interface and is defined as follow: the Equation (5) is initialized and resolved according to a fictitious time τ , 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 Ω :

Coupling of the Flow Equation with the Level Set Method
The hypothesis of two immiscible fluids leads to the notion of interface, which is Computational Water, Energy, and Environmental Engineering 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: with o φ a given function and n the unit normal vector at the interface going from s Ω to 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: where σ is the surface tension coefficient, n the unit normal vector to Γ going from s Ω to f Ω ,  the curvature of Γ and δ the Dirac function. The calculation of the normal and the curvature is done using the following relations: The flow equation in different phases is given by the relation (7). By integrating this equation on Ω , one has: Taking into account that s f = = − n n n , one has: For 0 Q = (no extraction, no infiltration to simplify), one obtains: The expression of the Darcy velocity according to the pressure at the interface gives: Equations (10) and (11) give: Thus, the equation modeling the flow of saline intrusion in a porous medium is The physical properties where the smoothed Heaviside function is defined by: and the Dirac function by: 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]:

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 ; ; and ik Ω the cells or controls volumes. One notices: It is assumed that space steps x ∆ and z ∆ are constants.

Discretization of the Level Set Equation
Spatial discretization: the discretization of Equation (8) using the finite volume method is given: where is the approximation of the numerical flux at the point ( ) 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φ = , then using the upwind scheme one gets:   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 3 r = represents the number of points in each sub-stencil, we have: where the nonlinear weights ( ) r j w are defined as follow:   ; The linear weights of the flux ( ) The expressions j v ± are calculated based on the sub-stencils.

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

Parameters
Values [33], [51] Porosity  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 dimen- ).
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 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 100 it = .
Figure 5(a) shows the permanent solutions of the level set function φ 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 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].  . 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 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 0.5 c = (0.5 isochlor distribution). . 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 100 it = . 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).

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.