# Parametrical Study of Freshwater–Saltwater Interface Dynamic

Authors: Mody Lo; Daouda Sangare
Journal: Modelling and Simulation in Engineering
 (2023)

---

## Abstract

In this paper, a saltwater intrusion model, in view to study the dynamics of the interface between the saltwater and the freshwater in a coastal aquifer, is established. This dynamic is caused by an injection of saltwater and a freshwater pumping through a well located in a given position. From the flow model in each phase, we defined an appropriate hypothesis to obtain our global model based only on the height of the interface. The numerical simulation of our model led us to study the effect of the parameters and obtain some empirical laws of the pollution time versus the distance well-injection area and pollution time versus the pumping flow.


---

## Body

## 1. Introduction

Water in general, and particularly the freshwater, is a scarce commodity and extremely useful to the human life. Only 3% of the quantity of water in earth is fresh. The aquifers are, for most of the countries, a source of supply in freshwater. Coastal areas are generally the most condensed population regions in the world [1–3]. Concentrated populations in those regions result in increased demand for freshwater and accelerated groundwater pumping, leading to groundwater depletion, especially in arid and semi-arid regions. To meet water needs for agriculture, industry, and public water supplies, groundwater resources have been seriously over exploited in the last several decades [4]. Globally, fresh groundwater resources in coastal aquifers are significantly impacted by seawater intrusion [5]. Seawater intrusion in coastal aquifers is a common problem and is encountered, with different degrees, in almost all coastal aquifers. It is regarded as a natural process that might be accelerated or retarded by external factors such as increase or decrease in the groundwater pumping, irrigation and recharge practices, land use, and possible seawater rise due to the impacts of global warming. The seawater intrusion problem has been under investigation for well over a century [6, 7]. A comprehensive review on different aspects of seawater intrusion assessment, monitoring, and modeling is provided by Bear et al. [8]. Physically, seawater intrusion is a density-dependent problem [9–16]. Modeling a seawater intrusion process needs to couple groundwater flow equation with solute (salt) transport equation [17], since the solution of salt transport is based on the groundwater flow field, which is in turn affected by salt and density distribution in the groundwater field.Over the years, many results have been established. An analytical solution for the steady-state salt distribution in a confined aquifer has been proposed by Henry [18]. A few years after, Jacob Bear [19] lays the foundations of this modeling. Segol et al. [9] developed the first transient solution based on the velocity-dependent dispersion coefficient using the Galerkin finite element method to solve the set of non-linear partial differential equations describing the movement of a saltwater front in a coastal confined aquifer [20]. A remarkable work has been achieved by Ghyben and Hergberg, by succeeding to make a relation between the height of the freshwater above the sea level h and the height of the freshwater under the sea level H (H=40h). Recent results have been due to [21, 22]. In [23], the author has proposed a mixed approach between sharp interface and diffuse interface. Others authors have also achieved great works. That’s the case of Lèye et al. [24] who developed a model of saltwater intrusion in a coastal aquifer by taking into account the high hydrodynamic dispersion of the saltwater creating, so a wide transition area between the freshwater and the saltwater. A similar model has been studied by Hamidi and Yazdi [20].In this paper, the problem of the seawater intrusion into the aquifer is studied by modeling the medium; therefore, the aquifer as a porous medium in which there is a two-phase flow. The two phases are separated by a freshwater-saltwater interface, which is considered as a contact surface. Following large water mass movements and freshwater pumping through a well, this interface can move, i.e., the shape and the position of the interface can vary. Our objective in this paper is to model the dynamic of this interface into the aquifer by an injection of saltwater and a freshwater pumping through a well located in a given position. The global model based only on the height of the interface is obtained from the flow model in each phase and an appropriate hypothesis. A mathematical analysis of the global model is done before the numerical simulation by a finite element method. If the interface elevation, due to the freshwater pumping through the well, reaches a certain threshold, we said that the well is polluted and we take this pollution time. A parametrical study of this pollution time according to the flow pumping and to the well position variation is done. Some empirical laws are obtained.The paper is organized as follows. In Section 2, we present the global model describing the dynamic of the freshwater-saltwater interface. The resolution of this mathematical model is done in Section 3. We start it by a mathematical analysis before the numerical resolution. A parametrical study and the determination of the different empirical laws end this section. The last section is devoted to the conclusion and some perspectives.


## 2. Model Description

In the present study, conceptual, unconfined coastal aquifer is considered as shown by a schematic section in Figure1. We have two phases: the freshwater and the saltwater, which is the sea. Between those two phases, we have the interface. Like shown in Figure 1, there is an injection of saltwater by the sea and a freshwater pumping by a well; this phenomenon involves the dynamic of the interface. We set then, for each position of this interface, the head by zint. So we have just to study the dynamic of this head. For this, we need to know the governing equations of the flow in each phase; we call it local model, and then the global model will translate the dynamic of the interface to the freshwater pumping and the injection of the saltwater.Figure 1
Schematic section of an unconfined coastal aquifer.We consider a representative domainΩ⊂ℝ3 composed of Ωs, the saltwater phase, and Ωf, the freshwater phase. Between the two phases, we have the interface. To obtain our model, we consider the continuity equations coupled with Darcy’s equation, which means mass conservation laws for each phase (fresh and saltwater) coupled with the classical Darcy law for porous media. The two phases are the same fluid with different characteristics. The flow is then governed by the same laws. Therefore, we consider just one phase to obtain the model, and for the second phase, the model is obtained by similarity.


### 2.1. Governing Equations in Each Phase: Local Models

Letθs be the saltwater content, ρs be the density, and Us be the velocity. The continuity equation is given by
(1)divρsUs+∂ρsθs∂t=ρsqsp−qst,where qsp (respectively qst) is the provided mass flow (respectively the taken mass flow).The effective velocityUs of the flow is thus related to the pressure P through the so-called Darcy law
(2)Us=−kμ∇ps+ρsg∇zs,where ρs and μ are respectively the density and the viscosity of the fluid, k is the permeability of the soil, and g is the gravitational acceleration constant. Introducing the piezometric head hs defined by
(3)hs=psρsg+zs,where zs is the elevation of the considering particle, we write equation (2) as follows:
(4)Us=−K∇hs,where K=kρsg/μ is the hydraulic conductivity, which expresses the ability of the underground to conduct the fluid. The saltwater content, given by θs=Vs/V, where Vs is the volume of Ωs and V is the volume of Ω, verifies
(5)dθs=1V1−θsdVs.We assume that the solid matrix of the aquifer is non-deformable. However, to take into account the contribution of its compressibility effect in the specific storage, we will assume that the solid matrix is elastic, i.e., there is a linear relationship between the effective compressive stress and the strain. The state equation of the subdomainΩs is then given by
(6)∂θs∂t=VsV1−θsα∂ps∂t,where ps and α are respectively the pressure of saltwater and the specific coefficient of compressibility of the media. The partial derivative of the pressure ps according the time is given by the following equation: from equation (3), we have ps=ρsghs−zs, therefore
(7)∂ps∂t=ghs−zs∂ρs∂t+ρsg∂hs∂t.Combining the state equation of the saltwater given as follows:(8)∂ρs∂t=ρsβs∂ps∂t,and equation (7), we deduce this following relation:
(9)1+zs−hsρsβsg∂ps∂t=ρsg∂hs∂t,with βs the compressibility coefficient of the saltwater. And since zs−hsρsβsg≪1, we obtain
(10)∂ps∂t=ρsg∂hs∂t.Developing the continuity equation (1) of the saltwater, we obtain:
(11)ρsdivUs+Us⋅∇ρs+ρs∂θs∂t+θs∂ρs∂t=ρsqsp−qst.Sinceρs≠0 and does not depend on the space variable, with relations (6) and (8), equation (11) becomes
(12)divUs+VsV1−θsα∂ps∂t+θsβs∂ps∂t=qsp−qst.Considering equation (10) and setting the specific storage coefficient of ΩS(13)Ss=ρsθsgβs+VsV1−θsθsα,equation (12) becomes
(14)divUs+Ss∂hs∂t=qsp−qst.With Darcy’s equation in each phase, we obtain the following system in the saltwater phase(15)divUs⟶+Ss∂hs∂t=qsp−qstin0,T×Ωs,Us⟶=−Ks∇hs,,and by similarity the governing system for the freshwater movement
(16)divUf⟶+Sf∂hf∂t=qfp−qftin0,T×Ωf,Uf⟶=−Kf∇hf,,where Sf is given like in equation (13) replacing the “s” indice by “f”. Being in the same homogeneous media Ω, we take Ks=Kf=k. The specific storage coefficient of the saltwater and the freshwater, Ss and Sf, respectively, depends on θs and θf, respectively. To close the models (15) and (16), we consider the case Ss and Sf are constants.


### 2.2. The Global Model

In the previous section, we have local systems that govern the flow in each compartmentΩs and Ωf of our study domain Ω. Now, we need to find a valid system throughout the domain Ω, which means a global model of the phenomena. For that, we set global variables using indicator function, which is defined as follows:
(17)χax=1ifx∈Ωa0elsewhere,like
(18)h=χfhf+χshs,S=χfSf+χsSs,U=χfUf+χsUsqp=χfqfp+χsqsp,qt=χfqft+χsqst..In this work, freshwater and saltwater are considered. Off course those two fluids are miscible. Therefore, they are separated by a transition zone characterized by the variations of the salt concentration. Nevertheless, the thickness of the transition zone is small compared to dimensions of the aquifer. We then assume that an abrupt interface separates two distinct domains, one for the saltwater and one for the freshwater. This interface is then considered like a contact surface, which means there is continuity of the pressure at the interface.The piezometric head in the saltwater,hs and in the freshwater, hf are given respectively by hs=ps/ρsg+zsandhf=pf/ρfg+zf.Let a particle at the interface with elevationzint, means zs=zf=zint, with the continuity of the pressure at the interface, ps=pf, we have ρsghs−zint=ρfghf−zintmeanszint=ρs/ρs−ρfhs−ρf/ρs−ρfhf.We obtain the following relation(19)zint=−δhf+1+δhs,withδ=ρfρs−ρf.Only pumping freshwater throughout the pumping wellω is considered, and let qf, the pumping term. We assume that there are not providing freshwater and not taking saltwater, i.e., qfp=0 and qst=0. Let Γs, the part of the boundary of Ωs, where the saltwater injection is done and let qs, the injection term. The dynamics of the interface between the saltwater and the freshwater can be obtained by solving the local models given by equations (15) and (16) and the interface head in equation (19). But in this work, we find a global system that gives directly the dynamics of the interface. We can remark that the dynamics of the interface drive the dynamics of the salt wedge and vice versa. Thereby, considering systems (15) and (16) and the expression of zint given in equation (19), we have the following equation
(20)S∂zint∂t−kΔzint=−δqfp−qft+1+δqsp−qst,where S is given in equation (18). Like in the local models, we study only the case Ss=Sf=S=constant. Under the pumping of freshwater and the injection of saltwater, we obtain the following system governing the movement of the interface
(21)S∂zint∂t−kΔzint=δqfp1win0,T×Ω,∂zint∂n=1+δqston0,T×Γs,zint0,x,y=zint,0x,yinΩ..The next step is devoted to the resolution of the system (21) for studying the movement of the interface.


## 2.1. Governing Equations in Each Phase: Local Models

Letθs be the saltwater content, ρs be the density, and Us be the velocity. The continuity equation is given by
(1)divρsUs+∂ρsθs∂t=ρsqsp−qst,where qsp (respectively qst) is the provided mass flow (respectively the taken mass flow).The effective velocityUs of the flow is thus related to the pressure P through the so-called Darcy law
(2)Us=−kμ∇ps+ρsg∇zs,where ρs and μ are respectively the density and the viscosity of the fluid, k is the permeability of the soil, and g is the gravitational acceleration constant. Introducing the piezometric head hs defined by
(3)hs=psρsg+zs,where zs is the elevation of the considering particle, we write equation (2) as follows:
(4)Us=−K∇hs,where K=kρsg/μ is the hydraulic conductivity, which expresses the ability of the underground to conduct the fluid. The saltwater content, given by θs=Vs/V, where Vs is the volume of Ωs and V is the volume of Ω, verifies
(5)dθs=1V1−θsdVs.We assume that the solid matrix of the aquifer is non-deformable. However, to take into account the contribution of its compressibility effect in the specific storage, we will assume that the solid matrix is elastic, i.e., there is a linear relationship between the effective compressive stress and the strain. The state equation of the subdomainΩs is then given by
(6)∂θs∂t=VsV1−θsα∂ps∂t,where ps and α are respectively the pressure of saltwater and the specific coefficient of compressibility of the media. The partial derivative of the pressure ps according the time is given by the following equation: from equation (3), we have ps=ρsghs−zs, therefore
(7)∂ps∂t=ghs−zs∂ρs∂t+ρsg∂hs∂t.Combining the state equation of the saltwater given as follows:(8)∂ρs∂t=ρsβs∂ps∂t,and equation (7), we deduce this following relation:
(9)1+zs−hsρsβsg∂ps∂t=ρsg∂hs∂t,with βs the compressibility coefficient of the saltwater. And since zs−hsρsβsg≪1, we obtain
(10)∂ps∂t=ρsg∂hs∂t.Developing the continuity equation (1) of the saltwater, we obtain:
(11)ρsdivUs+Us⋅∇ρs+ρs∂θs∂t+θs∂ρs∂t=ρsqsp−qst.Sinceρs≠0 and does not depend on the space variable, with relations (6) and (8), equation (11) becomes
(12)divUs+VsV1−θsα∂ps∂t+θsβs∂ps∂t=qsp−qst.Considering equation (10) and setting the specific storage coefficient of ΩS(13)Ss=ρsθsgβs+VsV1−θsθsα,equation (12) becomes
(14)divUs+Ss∂hs∂t=qsp−qst.With Darcy’s equation in each phase, we obtain the following system in the saltwater phase(15)divUs⟶+Ss∂hs∂t=qsp−qstin0,T×Ωs,Us⟶=−Ks∇hs,,and by similarity the governing system for the freshwater movement
(16)divUf⟶+Sf∂hf∂t=qfp−qftin0,T×Ωf,Uf⟶=−Kf∇hf,,where Sf is given like in equation (13) replacing the “s” indice by “f”. Being in the same homogeneous media Ω, we take Ks=Kf=k. The specific storage coefficient of the saltwater and the freshwater, Ss and Sf, respectively, depends on θs and θf, respectively. To close the models (15) and (16), we consider the case Ss and Sf are constants.


## 2.2. The Global Model

In the previous section, we have local systems that govern the flow in each compartmentΩs and Ωf of our study domain Ω. Now, we need to find a valid system throughout the domain Ω, which means a global model of the phenomena. For that, we set global variables using indicator function, which is defined as follows:
(17)χax=1ifx∈Ωa0elsewhere,like
(18)h=χfhf+χshs,S=χfSf+χsSs,U=χfUf+χsUsqp=χfqfp+χsqsp,qt=χfqft+χsqst..In this work, freshwater and saltwater are considered. Off course those two fluids are miscible. Therefore, they are separated by a transition zone characterized by the variations of the salt concentration. Nevertheless, the thickness of the transition zone is small compared to dimensions of the aquifer. We then assume that an abrupt interface separates two distinct domains, one for the saltwater and one for the freshwater. This interface is then considered like a contact surface, which means there is continuity of the pressure at the interface.The piezometric head in the saltwater,hs and in the freshwater, hf are given respectively by hs=ps/ρsg+zsandhf=pf/ρfg+zf.Let a particle at the interface with elevationzint, means zs=zf=zint, with the continuity of the pressure at the interface, ps=pf, we have ρsghs−zint=ρfghf−zintmeanszint=ρs/ρs−ρfhs−ρf/ρs−ρfhf.We obtain the following relation(19)zint=−δhf+1+δhs,withδ=ρfρs−ρf.Only pumping freshwater throughout the pumping wellω is considered, and let qf, the pumping term. We assume that there are not providing freshwater and not taking saltwater, i.e., qfp=0 and qst=0. Let Γs, the part of the boundary of Ωs, where the saltwater injection is done and let qs, the injection term. The dynamics of the interface between the saltwater and the freshwater can be obtained by solving the local models given by equations (15) and (16) and the interface head in equation (19). But in this work, we find a global system that gives directly the dynamics of the interface. We can remark that the dynamics of the interface drive the dynamics of the salt wedge and vice versa. Thereby, considering systems (15) and (16) and the expression of zint given in equation (19), we have the following equation
(20)S∂zint∂t−kΔzint=−δqfp−qft+1+δqsp−qst,where S is given in equation (18). Like in the local models, we study only the case Ss=Sf=S=constant. Under the pumping of freshwater and the injection of saltwater, we obtain the following system governing the movement of the interface
(21)S∂zint∂t−kΔzint=δqfp1win0,T×Ω,∂zint∂n=1+δqston0,T×Γs,zint0,x,y=zint,0x,yinΩ..The next step is devoted to the resolution of the system (21) for studying the movement of the interface.


## 3. Resolutions

In this section, we solve theoretically and numerically our global model (21). We start by a mathematical analysis, where we show that our problem is well posed. Before ending this section by a parametrical study, we expose the numerical solution, which shows the interface dynamic under the freshwater pumping and saltwater injection effects.


### 3.1. Mathematical Analysis

We consider an open bounded domainΩ of ℝ3 describing the interface. The boundary of Ω, assumed C1 is denoted Γ. The time interval is 0,T, T being any non-negative real number, and we set ΩT=0,T×Ω. The space of real values functions that are square integral on Ω with respect to the Lebesgue measure dx is denoted by the relationship
(22)L2Ω=u:Ω⟶ℝsuchthat∫Ωu2xdx<∞.The spaceL2Ω is equipped with the norm
(23)uL2Ω=∫Ωux2dx1/2,and the scalar product
(24)u,vL2Ω=∫Ωuxvxdx.The Sobolev spaceH1Ω is denoted by the expression
(25)H1Ω=u:Ω⟶ℝsuchthatu∈L2Ω,∇u∈L2Ω3.The spaceH1Ω is equipped with the scalar product
(26)u,vH1Ω=∫Ωuxvx+∇ux∇vxdx,that induces the norm
(27)uH1Ω=∫Ωux2+∇ux2dx1/2.We can remark thatH1Ω is a particular case of the Sobolev space Wm,pΩ, m and p integers, where
(28)Wm,pΩ=u:Ω⟶ℝsuchthatu∈LpΩ,Dαu∈LpΩ3,∀α=α1,α2,α3∈ℕ3:∣α∣≤m,Dαu is given by Dαu=∂∣α∣u/∂x1α1∂x2α2∂x3α3,x=x1,x2,x3∈Ω,and∣α∣=α1+α2+α3.IndeedH1Ω=W1,2Ω and so H2Ω=W2,2Ω.We defineL20,T,H2Ω by the set of all function u such that ut,⋅∈H2Ω,∀t∈0,T and u⋅,x∈L20,T∀x∈Ω.The following theorem shows that our model (21), under some assumptions, has an unique solution.Theorem 1.
Forqd in L20,T×w, qsinL20,T×Γs, andzint,0inL2Ω, the system (21) admits a unique solution inL20,T,H2Ω. This solution that we denote byzintsatisfies∂zint/∂n∈L20,T,L2Γsand∂zint/∂t∈L20,T,L2Ω.For elements of proof and more documentation, see [25] or [26].


### 3.2. Numerical Simulation and Parametrical Study

In this section, using some numerical schemes, we solve the problem (21). Like the dynamic of the interface depends on the parameters of the model, a parametrical study ends this section.


#### 3.2.1. Numerical Simulation

For the numerical simulation of our model, we use theP1 Lagrange finite element to deal with the spatial discretization of the problem (21). For that, it is convenient to write the variational formulation of the problem. It consists to replace the equation of the problem by an equivalent formula, said variational. This formula is obtained by multiplying the equation by a test function to integrate. The main tools for the variational formulation are the Green formula [26–28]. We obtain the following integral equation
(29)k∫Ωzint∇vdX+S∫Ω∂∂tzint⋅vdX−k∫Γs1+δqsvdσ=∫Ωδqf1wvdX,zint0,x=zint,0x.The time operator∂/∂tzint is approximated by an implicit Euler scheme ∂zint/∂t=zintm+1−zintm/dt.The variational formulation is given as follows:(30)∫ΩSzintm+1v+dt⋅k∇zintm+1∇vdX−dt⋅k∫Γs1+δqsvdσ=∫ΩSzintmv+dt⋅δqf1wvdX,zint0,x=zint,0x.The software are FreeFem++ for the resolution of the model and Python for visualization of the obtained results. The aquifer is represented by a three-dimensionalΩ of size −4,4×−4,4×−2,35. The freshwater pumping is done in a circular well ω, which is centered at 0,0,34.5, and the saltwater injection is done at the face Γs=4×−4,4×−2,20. We use the P1 finite elements with a structured mesh. The initial state of our domain is illustrated in Figure 2, where we have the interface alone in Figure 2(a) and the interface and the two fluids (freshwater and saltwater) in Figure 2(b).Figure 2
Schematic representation of the aquifer at the initial state. Data: (a) the full interface; (b) the whole aquifer with the saltwater in green, the freshwater in blue, and the interface in red.
(a)(b)To start solving our model, the following flow parametersqf=1.5, k=S=1, and qs=−1 are considered. The different results plotted in Figure 3 show a movement of the interface under the freshwater pumping effects. We notice that this movement is much more visible at the level of the pumping well, hence this kind of bump on the interface. In Figure 3, we can see that the dynamic of the interface is not uniform and depends on which position we are according to the well. In Figures 3(a), 3(b), 3(c), 3(d), 3(e), and 3(f), we plot the interface respectively in the section y=−1 and y=−2. The bump elevation due to the freshwater pumping is more remarkable at y=−1 than at y=−2. Moreover, this bump elevation depends also on the pumping duration. We notice it in Figures 3(a), 3(b), and 3(c) in the position y=−1 and in Figures 3(d), 3(e), and 3(f), for y=−2 for times t=40,50,and60, respectively. We can remark that the pumping flow plays an important role in the interface dynamic. This effect can be shown if we consider the pollution time of the well. Indeed, we saw that when we pump the freshwater, the interface moves and we have a bump on the well level; this bump can grow up to reach one certain head. If this bump level is equal to 34,5, we said there is pollution and stop pumping. We recall that the upper bound of the head of our domain is 35. The effect of pumping flow is shown by the pollution time versus the flow. This effect is plotted in Figure 4 for different positions of the freshwater pumping well. In Figure 4(a), results are obtained for the well located to the distance d=1.5 from the saltwater injection face Γs and in Figure 4(b) for d=3. In both figures, we remark that the pollution time depends deeply on the pumping flow. It is a decreasing function of the flow. We notice also that this pollution time, for a fixed flow, increases with the distance between the well and the saltwater injection area. It can be seen by comparing Figures 4(a) and 4(b). All these observed phenomena lead us to make a parametric study in view of drawing empirical laws.Figure 3
The non-uniform dynamic of the interface and the pumping duration effect. Data: (a), (b), and (c) are a vertical section aty=−1 for pumping duration t=40,50,and60, respectively. (d), (e), and (f) are a vertical section at y=−2 for pumping duration t=40,50,and60, respectively.
(a)(b)(c)(d)(e)(f)Figure 4
Pumping flow effect on the interface dynamic shown by the pollution time versus the flow for different fixed distance,d, between the well and Γs. Data: (a) d=1.5 and (b) d=3.
(a)(b)


#### 3.2.2. Parametrical Study

We have taken different pumping flows and different positions of our pumping well. The pollution threshold remains fixed atz=34.5, this means that when the interface elevation reaches this level, we say that the pumping well is polluted. According to these flows and positions, we have different pollution times. The flow is the product between the velocity (here noted by qd, otherwise our pumping term) and the well area, so for obtaining our different flows, we fixed the velocity qd and taking different radius r. In the different figures, the values of the flow are in reality the values of the well radius. The distances well-sea are the different distances between the well and the saltwater injection part.To study the pumping well position effects, we move the well for different fixed flows. The obtained results are plotted in Figure5(a). The pollution time is given versus the distance between the pumping well and the saltwater injection area for each fixed flow. For all fixed flow, the pollution time is small if we are close to the sea, on the other hand, if one moves away from the sea, this time increases, it means that one can exploit the well longer. Thus, for sustainable use of the well in coastal aquifer, it is very important to take into account the well positions and maximize as much as possible the distance between the pumping well and the sea. It is shown too in Figure 5(a) that the pollution time depends on the using flow, which justifies the different curves in this figure. To emphasize this effect, we plot the pollution time versus the flow for different fixed distance in Figure 5(b). In Figure 5(b), the flow effect on the pollution time is studied. For each fixed position of the pumping well, we pump with different flows. We recall that the different flows are obtained by taking different well radius. The plotted results in Figure 5(b) show that the pollution time is a decreasing function of the flow for each position, which means for each distance between the pumping well and the saltwater injection part. This shows us that if we want to use a pumping well for a long time, we must manage the flow. The appropriate flow for well position given can be known if we find a relation between the pollution time and the distance for a fixed flow or the pollution time and the flow for a given distance. This leads us to look for empirical laws between pollution time and the flow but also between pollution time and the distance between the pumping well and the saltwater injection part.Figure 5
Pumping flow and distance well-sea effect on the pollution time. Data: (a) pollution time versus distance well-sea for fixed flows given by the well radius; (⋄)r=0.2, (°)r=0.4, (□)r=0.6, (⊲)r=1; (b) pollution time versus the flow for different fixed distance, d, between the well and Γs; (⋄)d=1.5, (°)d=3, (□)d=4.5, (⊲)d=6.
(a)(b)


### 3.3. Empirical Laws

The results plotted in Figures5(a) and 5(b) lead us to think about a relationship between pollution time and the distance in one hand and pollution time and the flow in another hand. For that, we fit the different curves according to a logarithmic representation. We see that for fixed flows, the pollution time can be obtained by a power law function of the distance between the pumping well and the saltwater injection part. This power law function is given by equation (31)
(31)PT=expα×distanceβ.In another hand, for fixed distance, which means fixed position of the pumping well, the pollution time is also obtained by a power law function of the flow, which is given by equation (32)
(32)PT=expγ×flowλ.Of course, the coefficients in equations (31) and (32) depend on the fixed parameters of the model. In the empirical law (31), α and β depend on the fixed flow, and we remark that β is a positive number, which shows that our increasing function is like the curves in Figure 5(a). The coefficients γ and λ in the empirical law (32) depend also on the fixed distance, and λ is a negative number. So we have here a decreasing function of the flow, which is shown in Figure 5(b). A comparison between the pollution time obtained in simulation and the one calculated with the empirical laws is done in Figure 6 (pollution time versus distance for r=0.3 and r=0.5) and in Figure 7 (pollution time versus flow for d=2 and d=3.5). The results of this comparison are satisfactory.Figure 6
Pollution time versus the distance between the well and the saltwater injection part for different fixed flows: comparison between the simulation results and the empirical law. Data: (a)r=0.3, (b) r=0.5. In both, marked line is for simulation results, and continuous bold line is for empirical law.
(a)(b)Figure 7
Pollution time versus the pumping flow for different fixed distances: comparison between the simulation results and the empirical law. Data: (a)d=2, (b) d=3.5. In both, marked line is for simulation results, and continuous bold line is for empirical law.
(a)(b)


## 3.1. Mathematical Analysis

We consider an open bounded domainΩ of ℝ3 describing the interface. The boundary of Ω, assumed C1 is denoted Γ. The time interval is 0,T, T being any non-negative real number, and we set ΩT=0,T×Ω. The space of real values functions that are square integral on Ω with respect to the Lebesgue measure dx is denoted by the relationship
(22)L2Ω=u:Ω⟶ℝsuchthat∫Ωu2xdx<∞.The spaceL2Ω is equipped with the norm
(23)uL2Ω=∫Ωux2dx1/2,and the scalar product
(24)u,vL2Ω=∫Ωuxvxdx.The Sobolev spaceH1Ω is denoted by the expression
(25)H1Ω=u:Ω⟶ℝsuchthatu∈L2Ω,∇u∈L2Ω3.The spaceH1Ω is equipped with the scalar product
(26)u,vH1Ω=∫Ωuxvx+∇ux∇vxdx,that induces the norm
(27)uH1Ω=∫Ωux2+∇ux2dx1/2.We can remark thatH1Ω is a particular case of the Sobolev space Wm,pΩ, m and p integers, where
(28)Wm,pΩ=u:Ω⟶ℝsuchthatu∈LpΩ,Dαu∈LpΩ3,∀α=α1,α2,α3∈ℕ3:∣α∣≤m,Dαu is given by Dαu=∂∣α∣u/∂x1α1∂x2α2∂x3α3,x=x1,x2,x3∈Ω,and∣α∣=α1+α2+α3.IndeedH1Ω=W1,2Ω and so H2Ω=W2,2Ω.We defineL20,T,H2Ω by the set of all function u such that ut,⋅∈H2Ω,∀t∈0,T and u⋅,x∈L20,T∀x∈Ω.The following theorem shows that our model (21), under some assumptions, has an unique solution.Theorem 1.
Forqd in L20,T×w, qsinL20,T×Γs, andzint,0inL2Ω, the system (21) admits a unique solution inL20,T,H2Ω. This solution that we denote byzintsatisfies∂zint/∂n∈L20,T,L2Γsand∂zint/∂t∈L20,T,L2Ω.For elements of proof and more documentation, see [25] or [26].


## 3.2. Numerical Simulation and Parametrical Study

In this section, using some numerical schemes, we solve the problem (21). Like the dynamic of the interface depends on the parameters of the model, a parametrical study ends this section.


### 3.2.1. Numerical Simulation

For the numerical simulation of our model, we use theP1 Lagrange finite element to deal with the spatial discretization of the problem (21). For that, it is convenient to write the variational formulation of the problem. It consists to replace the equation of the problem by an equivalent formula, said variational. This formula is obtained by multiplying the equation by a test function to integrate. The main tools for the variational formulation are the Green formula [26–28]. We obtain the following integral equation
(29)k∫Ωzint∇vdX+S∫Ω∂∂tzint⋅vdX−k∫Γs1+δqsvdσ=∫Ωδqf1wvdX,zint0,x=zint,0x.The time operator∂/∂tzint is approximated by an implicit Euler scheme ∂zint/∂t=zintm+1−zintm/dt.The variational formulation is given as follows:(30)∫ΩSzintm+1v+dt⋅k∇zintm+1∇vdX−dt⋅k∫Γs1+δqsvdσ=∫ΩSzintmv+dt⋅δqf1wvdX,zint0,x=zint,0x.The software are FreeFem++ for the resolution of the model and Python for visualization of the obtained results. The aquifer is represented by a three-dimensionalΩ of size −4,4×−4,4×−2,35. The freshwater pumping is done in a circular well ω, which is centered at 0,0,34.5, and the saltwater injection is done at the face Γs=4×−4,4×−2,20. We use the P1 finite elements with a structured mesh. The initial state of our domain is illustrated in Figure 2, where we have the interface alone in Figure 2(a) and the interface and the two fluids (freshwater and saltwater) in Figure 2(b).Figure 2
Schematic representation of the aquifer at the initial state. Data: (a) the full interface; (b) the whole aquifer with the saltwater in green, the freshwater in blue, and the interface in red.
(a)(b)To start solving our model, the following flow parametersqf=1.5, k=S=1, and qs=−1 are considered. The different results plotted in Figure 3 show a movement of the interface under the freshwater pumping effects. We notice that this movement is much more visible at the level of the pumping well, hence this kind of bump on the interface. In Figure 3, we can see that the dynamic of the interface is not uniform and depends on which position we are according to the well. In Figures 3(a), 3(b), 3(c), 3(d), 3(e), and 3(f), we plot the interface respectively in the section y=−1 and y=−2. The bump elevation due to the freshwater pumping is more remarkable at y=−1 than at y=−2. Moreover, this bump elevation depends also on the pumping duration. We notice it in Figures 3(a), 3(b), and 3(c) in the position y=−1 and in Figures 3(d), 3(e), and 3(f), for y=−2 for times t=40,50,and60, respectively. We can remark that the pumping flow plays an important role in the interface dynamic. This effect can be shown if we consider the pollution time of the well. Indeed, we saw that when we pump the freshwater, the interface moves and we have a bump on the well level; this bump can grow up to reach one certain head. If this bump level is equal to 34,5, we said there is pollution and stop pumping. We recall that the upper bound of the head of our domain is 35. The effect of pumping flow is shown by the pollution time versus the flow. This effect is plotted in Figure 4 for different positions of the freshwater pumping well. In Figure 4(a), results are obtained for the well located to the distance d=1.5 from the saltwater injection face Γs and in Figure 4(b) for d=3. In both figures, we remark that the pollution time depends deeply on the pumping flow. It is a decreasing function of the flow. We notice also that this pollution time, for a fixed flow, increases with the distance between the well and the saltwater injection area. It can be seen by comparing Figures 4(a) and 4(b). All these observed phenomena lead us to make a parametric study in view of drawing empirical laws.Figure 3
The non-uniform dynamic of the interface and the pumping duration effect. Data: (a), (b), and (c) are a vertical section aty=−1 for pumping duration t=40,50,and60, respectively. (d), (e), and (f) are a vertical section at y=−2 for pumping duration t=40,50,and60, respectively.
(a)(b)(c)(d)(e)(f)Figure 4
Pumping flow effect on the interface dynamic shown by the pollution time versus the flow for different fixed distance,d, between the well and Γs. Data: (a) d=1.5 and (b) d=3.
(a)(b)


### 3.2.2. Parametrical Study

We have taken different pumping flows and different positions of our pumping well. The pollution threshold remains fixed atz=34.5, this means that when the interface elevation reaches this level, we say that the pumping well is polluted. According to these flows and positions, we have different pollution times. The flow is the product between the velocity (here noted by qd, otherwise our pumping term) and the well area, so for obtaining our different flows, we fixed the velocity qd and taking different radius r. In the different figures, the values of the flow are in reality the values of the well radius. The distances well-sea are the different distances between the well and the saltwater injection part.To study the pumping well position effects, we move the well for different fixed flows. The obtained results are plotted in Figure5(a). The pollution time is given versus the distance between the pumping well and the saltwater injection area for each fixed flow. For all fixed flow, the pollution time is small if we are close to the sea, on the other hand, if one moves away from the sea, this time increases, it means that one can exploit the well longer. Thus, for sustainable use of the well in coastal aquifer, it is very important to take into account the well positions and maximize as much as possible the distance between the pumping well and the sea. It is shown too in Figure 5(a) that the pollution time depends on the using flow, which justifies the different curves in this figure. To emphasize this effect, we plot the pollution time versus the flow for different fixed distance in Figure 5(b). In Figure 5(b), the flow effect on the pollution time is studied. For each fixed position of the pumping well, we pump with different flows. We recall that the different flows are obtained by taking different well radius. The plotted results in Figure 5(b) show that the pollution time is a decreasing function of the flow for each position, which means for each distance between the pumping well and the saltwater injection part. This shows us that if we want to use a pumping well for a long time, we must manage the flow. The appropriate flow for well position given can be known if we find a relation between the pollution time and the distance for a fixed flow or the pollution time and the flow for a given distance. This leads us to look for empirical laws between pollution time and the flow but also between pollution time and the distance between the pumping well and the saltwater injection part.Figure 5
Pumping flow and distance well-sea effect on the pollution time. Data: (a) pollution time versus distance well-sea for fixed flows given by the well radius; (⋄)r=0.2, (°)r=0.4, (□)r=0.6, (⊲)r=1; (b) pollution time versus the flow for different fixed distance, d, between the well and Γs; (⋄)d=1.5, (°)d=3, (□)d=4.5, (⊲)d=6.
(a)(b)


## 3.2.1. Numerical Simulation

For the numerical simulation of our model, we use theP1 Lagrange finite element to deal with the spatial discretization of the problem (21). For that, it is convenient to write the variational formulation of the problem. It consists to replace the equation of the problem by an equivalent formula, said variational. This formula is obtained by multiplying the equation by a test function to integrate. The main tools for the variational formulation are the Green formula [26–28]. We obtain the following integral equation
(29)k∫Ωzint∇vdX+S∫Ω∂∂tzint⋅vdX−k∫Γs1+δqsvdσ=∫Ωδqf1wvdX,zint0,x=zint,0x.The time operator∂/∂tzint is approximated by an implicit Euler scheme ∂zint/∂t=zintm+1−zintm/dt.The variational formulation is given as follows:(30)∫ΩSzintm+1v+dt⋅k∇zintm+1∇vdX−dt⋅k∫Γs1+δqsvdσ=∫ΩSzintmv+dt⋅δqf1wvdX,zint0,x=zint,0x.The software are FreeFem++ for the resolution of the model and Python for visualization of the obtained results. The aquifer is represented by a three-dimensionalΩ of size −4,4×−4,4×−2,35. The freshwater pumping is done in a circular well ω, which is centered at 0,0,34.5, and the saltwater injection is done at the face Γs=4×−4,4×−2,20. We use the P1 finite elements with a structured mesh. The initial state of our domain is illustrated in Figure 2, where we have the interface alone in Figure 2(a) and the interface and the two fluids (freshwater and saltwater) in Figure 2(b).Figure 2
Schematic representation of the aquifer at the initial state. Data: (a) the full interface; (b) the whole aquifer with the saltwater in green, the freshwater in blue, and the interface in red.
(a)(b)To start solving our model, the following flow parametersqf=1.5, k=S=1, and qs=−1 are considered. The different results plotted in Figure 3 show a movement of the interface under the freshwater pumping effects. We notice that this movement is much more visible at the level of the pumping well, hence this kind of bump on the interface. In Figure 3, we can see that the dynamic of the interface is not uniform and depends on which position we are according to the well. In Figures 3(a), 3(b), 3(c), 3(d), 3(e), and 3(f), we plot the interface respectively in the section y=−1 and y=−2. The bump elevation due to the freshwater pumping is more remarkable at y=−1 than at y=−2. Moreover, this bump elevation depends also on the pumping duration. We notice it in Figures 3(a), 3(b), and 3(c) in the position y=−1 and in Figures 3(d), 3(e), and 3(f), for y=−2 for times t=40,50,and60, respectively. We can remark that the pumping flow plays an important role in the interface dynamic. This effect can be shown if we consider the pollution time of the well. Indeed, we saw that when we pump the freshwater, the interface moves and we have a bump on the well level; this bump can grow up to reach one certain head. If this bump level is equal to 34,5, we said there is pollution and stop pumping. We recall that the upper bound of the head of our domain is 35. The effect of pumping flow is shown by the pollution time versus the flow. This effect is plotted in Figure 4 for different positions of the freshwater pumping well. In Figure 4(a), results are obtained for the well located to the distance d=1.5 from the saltwater injection face Γs and in Figure 4(b) for d=3. In both figures, we remark that the pollution time depends deeply on the pumping flow. It is a decreasing function of the flow. We notice also that this pollution time, for a fixed flow, increases with the distance between the well and the saltwater injection area. It can be seen by comparing Figures 4(a) and 4(b). All these observed phenomena lead us to make a parametric study in view of drawing empirical laws.Figure 3
The non-uniform dynamic of the interface and the pumping duration effect. Data: (a), (b), and (c) are a vertical section aty=−1 for pumping duration t=40,50,and60, respectively. (d), (e), and (f) are a vertical section at y=−2 for pumping duration t=40,50,and60, respectively.
(a)(b)(c)(d)(e)(f)Figure 4
Pumping flow effect on the interface dynamic shown by the pollution time versus the flow for different fixed distance,d, between the well and Γs. Data: (a) d=1.5 and (b) d=3.
(a)(b)


## 3.2.2. Parametrical Study

We have taken different pumping flows and different positions of our pumping well. The pollution threshold remains fixed atz=34.5, this means that when the interface elevation reaches this level, we say that the pumping well is polluted. According to these flows and positions, we have different pollution times. The flow is the product between the velocity (here noted by qd, otherwise our pumping term) and the well area, so for obtaining our different flows, we fixed the velocity qd and taking different radius r. In the different figures, the values of the flow are in reality the values of the well radius. The distances well-sea are the different distances between the well and the saltwater injection part.To study the pumping well position effects, we move the well for different fixed flows. The obtained results are plotted in Figure5(a). The pollution time is given versus the distance between the pumping well and the saltwater injection area for each fixed flow. For all fixed flow, the pollution time is small if we are close to the sea, on the other hand, if one moves away from the sea, this time increases, it means that one can exploit the well longer. Thus, for sustainable use of the well in coastal aquifer, it is very important to take into account the well positions and maximize as much as possible the distance between the pumping well and the sea. It is shown too in Figure 5(a) that the pollution time depends on the using flow, which justifies the different curves in this figure. To emphasize this effect, we plot the pollution time versus the flow for different fixed distance in Figure 5(b). In Figure 5(b), the flow effect on the pollution time is studied. For each fixed position of the pumping well, we pump with different flows. We recall that the different flows are obtained by taking different well radius. The plotted results in Figure 5(b) show that the pollution time is a decreasing function of the flow for each position, which means for each distance between the pumping well and the saltwater injection part. This shows us that if we want to use a pumping well for a long time, we must manage the flow. The appropriate flow for well position given can be known if we find a relation between the pollution time and the distance for a fixed flow or the pollution time and the flow for a given distance. This leads us to look for empirical laws between pollution time and the flow but also between pollution time and the distance between the pumping well and the saltwater injection part.Figure 5
Pumping flow and distance well-sea effect on the pollution time. Data: (a) pollution time versus distance well-sea for fixed flows given by the well radius; (⋄)r=0.2, (°)r=0.4, (□)r=0.6, (⊲)r=1; (b) pollution time versus the flow for different fixed distance, d, between the well and Γs; (⋄)d=1.5, (°)d=3, (□)d=4.5, (⊲)d=6.
(a)(b)


## 3.3. Empirical Laws

The results plotted in Figures5(a) and 5(b) lead us to think about a relationship between pollution time and the distance in one hand and pollution time and the flow in another hand. For that, we fit the different curves according to a logarithmic representation. We see that for fixed flows, the pollution time can be obtained by a power law function of the distance between the pumping well and the saltwater injection part. This power law function is given by equation (31)
(31)PT=expα×distanceβ.In another hand, for fixed distance, which means fixed position of the pumping well, the pollution time is also obtained by a power law function of the flow, which is given by equation (32)
(32)PT=expγ×flowλ.Of course, the coefficients in equations (31) and (32) depend on the fixed parameters of the model. In the empirical law (31), α and β depend on the fixed flow, and we remark that β is a positive number, which shows that our increasing function is like the curves in Figure 5(a). The coefficients γ and λ in the empirical law (32) depend also on the fixed distance, and λ is a negative number. So we have here a decreasing function of the flow, which is shown in Figure 5(b). A comparison between the pollution time obtained in simulation and the one calculated with the empirical laws is done in Figure 6 (pollution time versus distance for r=0.3 and r=0.5) and in Figure 7 (pollution time versus flow for d=2 and d=3.5). The results of this comparison are satisfactory.Figure 6
Pollution time versus the distance between the well and the saltwater injection part for different fixed flows: comparison between the simulation results and the empirical law. Data: (a)r=0.3, (b) r=0.5. In both, marked line is for simulation results, and continuous bold line is for empirical law.
(a)(b)Figure 7
Pollution time versus the pumping flow for different fixed distances: comparison between the simulation results and the empirical law. Data: (a)d=2, (b) d=3.5. In both, marked line is for simulation results, and continuous bold line is for empirical law.
(a)(b)


## 4. Conclusion

The work done in this paper is very interesting, and the results obtained are very important. The global model, which governs the dynamic of the interface between the saltwater and the freshwater, is obtained considering the flow model in each phase. Therefore, the salt concentration of each fluid is not considered, then we have no transport equations like the previous works in this field. The dynamic of the interface is due to an injection of saltwater by the sea and a freshwater pumping through a well. In this work, only one well is considered, and the injection flow of saltwater is constant. The pollution time is studied for several well positions and pumping flows. We can conclude that this study is very important in coastal aquifer. In fact, for an efficient use of the pumping wells, it is interesting to consider the distance between the well position and the saltwater injection part because pollution time is seen to be an increasing function of this distance for any fixed pumping flow. Moreover, a good management of the pumping flow can help for a longer use of the well. The reason is for a fixed distance, the pollution time is a decreasing function of the pumping flow. Empirical laws are found for the pollution time versus the distance respectively versus the flow for fixed flow respectively fixed distance. The results obtained here for saltwater intrusion can be used in other field, where we have two phases of fluid in the similar conditions. For future works, we can consider two or several pumping wells and compare the results with the use of a single well. The saltwater injection, which is taken as constant for any time, can be considered as varying and a control of the injections, and the pumping flows can be done. To obtain our global model, we assumed the interface as a contact surface, otherwise a pressure jump is to be expected.



---

*Source: 2715152-2023-03-10.xml*