Parametrical Study of Freshwater – Saltwater Interface Dynamic

,


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][2][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 densitydependent problem [9][10][11][12][13][14][15][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 twophase flow.The two phases are separated by a freshwatersaltwater 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.

Model Description
In the present study, conceptual, unconfined coastal aquifer is considered as shown by a schematic section in Figure 1.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 z int .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.
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 U s be the velocity.The continuity equation is given by where q p s (respectively q t s ) is the provided mass flow (respectively the taken mass flow).
The effective velocity U s of the flow is thus related to the pressure P through the so-called Darcy law 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 h s defined by where z s is the elevation of the considering particle, we write equation (2) as follows: where K = kρ s g/μ is the hydraulic conductivity, which expresses the ability of the underground to conduct the fluid.
The saltwater content, given by θ s = V s /V, where V s is the volume of Ω s and V is the volume of Ω, verifies We assume that the solid matrix of the aquifer is nondeformable.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  Modelling and Simulation in Engineering the strain.The state equation of the subdomain Ω s is then given by where p s and α are respectively the pressure of saltwater and the specific coefficient of compressibility of the media.The partial derivative of the pressure p s according the time is given by the following equation: from equation (3), we have p s = ρ s gðh s − z s Þ, therefore Combining the state equation of the saltwater given as follows: and equation ( 7), we deduce this following relation: with β s the compressibility coefficient of the saltwater.And since ðz s − h s Þρ s β s g ≪ 1, we obtain Developing the continuity equation ( 1) of the saltwater, we obtain: Since ρ s ≠ 0 and does not depend on the space variable, with relations ( 6) and (8), equation ( 11) becomes Considering equation (10) and setting the specific storage coefficient of With Darcy's equation in each phase, we obtain the following system in the saltwater phase and by similarity the governing system for the freshwater where S f is given like in equation ( 13) replacing the "s" indice by "f".Being in the same homogeneous media Ω, we take K s = K f = k.The specific storage coefficient of the saltwater and the freshwater, S s and S f , respectively, depends on θ s and θ f , respectively.
To close the models ( 15) and ( 16), we consider the case S s and S f 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: 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, h s and in the freshwater, h f are given respectively by h s = ðp s /ρ s gÞ + z s and h f = ðp f /ρ f gÞ + z f .
Let a particle at the interface with elevation z int , means z s = z f = z int , with the continuity of the pressure at the    Modelling and Simulation in Engineering interface, p s = p f , we have We obtain the following relation Only pumping freshwater throughout the pumping well ω is considered, and let q f , the pumping term.We assume that there are not providing freshwater and not taking saltwater, i.e., q p f = 0 and q t s = 0. Let Γ s , the part of the boundary of Ω s , where the saltwater injection is done and let q s , 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 z int given in equation ( 19), we have the following equation where S is given in equation (18).Like in the local models, we study only the case S s = S f = S = constant.Under the pumping of freshwater and the injection of saltwater, we obtain the following system governing the movement of the interface The next step is devoted to the resolution of the system (21) for studying the movement of the interface.

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.

Mathematical Analysis.
We consider an open bounded domain Ω of ℝ 3 describing the interface.The boundary of Ω, assumed C 1 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

Modelling and Simulation in Engineering
The space L 2 ðΩÞ is equipped with the norm and the scalar product The Sobolev space H 1 ðΩÞ is denoted by the expression The space H 1 ðΩÞ is equipped with the scalar product that induces the norm We can remark that H 1 ðΩÞ is a particular case of the Sobolev space W m,p ðΩÞ, m and p integers, where Indeed H 1 ðΩÞ = W 1,2 ðΩÞ and so H 2 ðΩÞ = W 2,2 ðΩÞ.We define L 2 ð0, T, H 2 ðΩÞÞ by the set of all function u such that uðt, ⋅ Þ ∈ H 2 ðΩÞ, ∀t ∈ ½0, T and uð⋅ ,xÞ ∈ L 2 ð½0, TÞ ∀x ∈ Ω.
The following theorem shows that our model (21), under some assumptions, has an unique solution.
For elements of proof and more documentation, see [25] or [26].

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.

Numerical Simulation.
For the numerical simulation of our model, we use the P 1 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][27][28].We obtain the following integral equation The time operator ð∂/∂tÞz int is approximated by an implicit Euler scheme ∂z int /∂t = z m+1 int − z m int /dt.The variational formulation is given as follows: 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 = f4g × ½−4, 4 × ½−2, 20.We use the P 1 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).
To start solving our model, the following flow parameters q f = 1:5, k = S = 1, and q s = −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, and 60, 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 Modelling and Simulation in Engineering 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.

Parametrical Study.
We have taken different pumping flows and different positions of our pumping well.The pollution threshold remains fixed at z = 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 q d , otherwise our pumping term) and the well area, so for obtaining our different flows, we fixed the velocity q d 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 wellsea 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 Figure 5(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.

Empirical Laws.
The results plotted in Figures 5(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) 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) 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

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 9 Modelling and Simulation in Engineering 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.

Figure 1 :
Figure 1: Schematic section of an unconfined coastal aquifer.

Figure 2 :
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.

Figure 6 :
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.

Figure 7 :
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.
(pollution time versus flow for d = 2 and d = 3:5).The results of this comparison are satisfactory.