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.

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 [13]. 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 [916]. 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 and the height of the freshwater under the sea level (). 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 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 . 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 composed of , the saltwater phase, and , 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 be the saltwater content, be the density, and be the velocity. The continuity equation is given by where (respectively ) is the provided mass flow (respectively the taken mass flow).

The effective velocity of the flow is thus related to the pressure P through the so-called Darcy law where and are respectively the density and the viscosity of the fluid, is the permeability of the soil, and g is the gravitational acceleration constant. Introducing the piezometric head defined by where is the elevation of the considering particle, we write equation (2) as follows: where is the hydraulic conductivity, which expresses the ability of the underground to conduct the fluid. The saltwater content, given by , where is the volume of and is the volume of , verifies

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 is then given by where and are respectively the pressure of saltwater and the specific coefficient of compressibility of the media. The partial derivative of the pressure according the time is given by the following equation: from equation (3), we have , therefore

Combining the state equation of the saltwater given as follows: and equation (7), we deduce this following relation: with the compressibility coefficient of the saltwater. And since , we obtain

Developing the continuity equation (1) of the saltwater, we obtain:

Since 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 equation (12) becomes

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 movement where is given like in equation (13) replacing the “s” indice by “f”. Being in the same homogeneous media , we take . The specific storage coefficient of the saltwater and the freshwater, and , respectively, depends on and , respectively. To close the models (15) and (16), we consider the case and are constants.

2.2. The Global Model

In the previous section, we have local systems that govern the flow in each compartment and 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: like

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, and in the freshwater, are given respectively by .

Let a particle at the interface with elevation , means , with the continuity of the pressure at the interface, , we have .

We obtain the following relation

Only pumping freshwater throughout the pumping well is considered, and let , the pumping term. We assume that there are not providing freshwater and not taking saltwater, i.e., and . Let , the part of the boundary of , where the saltwater injection is done and let , 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 given in equation (19), we have the following equation where is given in equation (18). Like in the local models, we study only the case . 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.

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 describing the interface. The boundary of , assumed C1 is denoted . The time interval is , being any non-negative real number, and we set . The space of real values functions that are square integral on with respect to the Lebesgue measure is denoted by the relationship

The space is equipped with the norm and the scalar product

The Sobolev space is denoted by the expression

The space is equipped with the scalar product that induces the norm

We can remark that is a particular case of the Sobolev space , and integers, where is given by .

Indeed and so .

We define by the set of all function such that and .

The following theorem shows that our model (21), under some assumptions, has an unique solution.

Theorem 1. For in , in, andin, the system (21) admits a unique solution in. This solution that we denote bysatisfiesand.

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 the 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

The time operator is approximated by an implicit Euler scheme .

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 . The freshwater pumping is done in a circular well , which is centered at , and the saltwater injection is done at the face . We use the 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 , , and 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 and . The bump elevation due to the freshwater pumping is more remarkable at than at . 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 and in Figures 3(d), 3(e), and 3(f), for for times , 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 from the saltwater injection face and in Figure 4(b) for . 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.

3.2.2. Parametrical Study

We have taken different pumping flows and different positions of our pumping well. The pollution threshold remains fixed at , 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 , otherwise our pumping term) and the well area, so for obtaining our different flows, we fixed the velocity and taking different radius . 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 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.

3.3. 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 and ) and in Figure 7 (pollution time versus flow for and ). The results of this comparison are satisfactory.

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.

Data Availability

Data supporting this research article are available from the corresponding author or first author on reasonable request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research work was supported by UFR SAT of Gaston BERGER University. We appreciate the several discussions with colleagues who helped to improve the quality of this work.