Modeling of the Saltwater Intrusion Using the Level Set Method. Application to Henry’s 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]:
(1)
where
,
,
,
and
are respectively the specific storage coefficient, the hydraulic head, the hydraulic conductivity, the source term and the Darcy velocity in the different phases
(
represents the freshwater phase;
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
which bounds an open domain
. This interface is defined at a point
of a domain
at time t, like the zero isocontours of a scalar function
(Figure 1).
In 2D the evolution interface
separating the two sub-domains
and
is defined by:
(2)
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,
, by using the sign of the level set function, i.e.:
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:
(3)
where d is the signed distance function.
The initial position
of the interface being known, its evolution is calculated using a given velocity field
. 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
is:
(4)
with a given initial condition
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 property i.e. if
, 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:
(5)
where
represents a fictitious reinitialization time of the function
, d represents the level set function at fictitious time
, and the
function represents the sign of the level set function
. This function is smoothed around the interface and is defined as follow:
(6)
where
and
are the steps of space discretization. At each time step
, 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
:
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:
(7)
and
(8)
with
a given function and
the unit normal vector at the interface going from
to
.
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:
(9)
where
is the surface tension coefficient,
the unit normal vector to
going from
to
,
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
, one has:
For
(no extraction, no infiltration to simplify), one obtains:
(10)
The expression of the Darcy velocity according to the pressure at the interface gives:
According to Laplace’s law, the discontinuity of the pressure at the interface is taken into account by the relation:
and according to the approximation of the Ghost Fluid Method (GFM) applied to the CSF method [41], one has:
Putting this expression in the previous relation, one obtains the expression of the velocity at the interface:
(11)
Equations (10) and (11) give:
Thus, the equation modeling the flow of saline intrusion in a porous medium is
(12)
The physical properties
and
are expressed using the level set function
, following the equations:
(13)
where the smoothed Heaviside function is defined by:
(14)
and the Dirac function by:
(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]:
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
and
a family of rectangles, such as:
with
,
and
the cells or controls volumes. One notices:
for
, and
for
, and
with
and
the coordinates of the centers of the cells
.
and
It is assumed that space steps
and
are constants.
3.1. Discretization of the Level Set Equation
Spatial discretization: the discretization of Equation (8) using the finite volume method is given:
(16)
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
, then using the upwind scheme one gets:
The calculation of
(respectively
) using the WENO5 scheme is done as follows:
(17)
with
(respectively
) the interpolation of the flux on each sub-stencils and r the number of points in sub-stencils,
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
represents the number of points in each sub-stencil, we have:
(18)
where the nonlinear weights
are defined as follow:
(19)
where
are the linear weights,
the smoothness indicators used for the weight calculations and
the small constant allowing to avoid that the denominator of
let be equal zero.
The parameters
and
are defined below (Equation (24) [42] ). Consequently, one has:
![]()
Figure 2. Stencil for
: presentation of the sub-stencils
making up the stencils.
![]()
Figure 3. Stencil for
: presentation of the sub-stencils
making up the stencils.
with
The linear weights of the flux
,
. And
are weighted so that
. For this, one takes
,
and
[42].
The success of WENO schemes lies in the choice of standardized weights (
). For consistency and stability of the scheme, these coefficients must satisfy the relations:
(20)
where
represents the mean value at the points of indices
[28].
Now the flux
for each sub-stencil are obtained as follows:
(21)
The expressions
are calculated based on the sub-stencils.
,
,
,
and
, are the mean values of
at points
.
For sub-stencils associated with
, i.e.
,
and
, one has:
(22)
For sub-stencils associated with
, i.e.
,
and
, one has:
(23)
The expressions
represent the parameters allowing to measure the regularity of the function on each sub-stencil. They are given by:
(24)
To obtain the flux
on the stencil
, one uses the same process that for the flux
.
Time discretization: Let L be a discrete operator,
the time step. The time discretization of the transport Equation (17) using the explicit Euler scheme gives:
with
Thus, the third-order TVD Runge-Kutta scheme, is:
(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
for all
and
an
, one gets:
(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
;
b) for each time step:
● calculation of the geometric properties of the interface such as the normal vector
and the curvature
, as well as the Heaviside function
, the Dirac function
, the storage coefficient
and the density
;
● resolution of the flow equation;
● calculation of the velocity knowing that
, with
the Darcy velocity and
the porosity of the domain
;
● resolution of the level set equation, using the velocity
obtained in step 3;
● application of the reinitialization algorithm to the level set function
if the distance property is not respected i.e. if
;
● update of the level set function
and of the hydraulic head H;
c) end if the initialized time is equal to the final time i.e.
, 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
, and Figure 4 allows to adopt the following initialization:
(27)
At initial condition, the interface is represented by the vertical axis
. So, the interface between freshwater and saltwater is given by:
For the boundary conditions of the level set function, one adopts the choice of homogeneous Neumann conditions at bounds
and
, i.e.
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
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
used by Lee and Cheng, as well as Henry. Some studies did not take into account this porosity and use the value
. Therefore, the solution obtained by Lee and Cheng, and Henry (with
) cannot be compared to the solution obtained by Pinder and Cooper, Segol et al., Frind and Huyakorn (with
).
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
.
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
) and the tip (intersection of the interface with the bound
) 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
and the hydraulic head H at time
days obtained using a reinitialization with the ENO1 scheme (number of iterations
). (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
days. (a) Result obtained by the level set method; (b) Interfaces obtained with the ENO1 scheme for
and
.
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
and the uncoupled solution of Simpson and Clement coincides exactly at their intersections with the bounds of the domain (axes
and
) and over the interval
. 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
(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
at the left boundary side of the domain
by
. 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
.
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.