Next Article in Journal
The Use of Chemical Methods and Magnetic Field in Conditioning and Dewatering of Digested Sewage Sludge
Next Article in Special Issue
Application of Fractional Flow Theory for Analytical Modeling of Surfactant Flooding, Polymer Flooding, and Surfactant/Polymer Flooding for Chemical Enhanced Oil Recovery
Previous Article in Journal
Structure and Physicochemical Properties of Water Treated under Methane with Low-Temperature Glow Plasma of Low Frequency
Previous Article in Special Issue
B-Spline Method of Lines for Simulation of Contaminant Transport in Groundwater
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Fully Implicit Finite Volume Scheme for a Seawater Intrusion Problem in Coastal Aquifers

by
Abdelkrim Aharmouch
1,
Brahim Amaziane
2,*,
Mustapha El Ossmani
3 and
Khadija Talali
3
1
LISAC, Université Sidi Mohamed Ben Abdellah, Faculté des Sciences (Dhar El Mehraz), Fès 30003, Morocco
2
LMAP, Universite de Pau et des Pays de l’Adour, CNRS, E2S UPPA, 64000 Pau, France
3
L2M3S-ENSAM, Université Moulay Ismaïl, Meknès 50000, Morocco
*
Author to whom correspondence should be addressed.
Water 2020, 12(6), 1639; https://doi.org/10.3390/w12061639
Submission received: 21 April 2020 / Revised: 25 May 2020 / Accepted: 29 May 2020 / Published: 8 June 2020

Abstract

:
We present a numerical framework for efficiently simulating seawater flow in coastal aquifers using a finite volume method. The mathematical model consists of coupled and nonlinear partial differential equations. Difficulties arise from the nonlinear structure of the system and the complexity of natural fields, which results in complex aquifer geometries and heterogeneity in the hydraulic parameters. When numerically solving such a model, due to the mentioned feature, attempts to explicitly perform the time integration result in an excessively restricted stability condition on time step. An implicit method, which calculates the flow dynamics at each time step, is needed to overcome the stability problem of the time integration and mass conservation. A fully implicit finite volume scheme is developed to discretize the coupled system that allows the use of much longer time steps than explicit schemes. We have developed and implemented this scheme in a new module in the context of the open source platform DuMu X . The accuracy and effectiveness of this new module are demonstrated through numerical investigation for simulating the displacement of the sharp interface between saltwater and freshwater in groundwater flow. Lastly, numerical results of a realistic test case are presented to prove the efficiency and the performance of the method.
MSC:
76S05; 65N08; 81T80; 97N80; 97P40

1. Introduction

The numerical modeling and analysis of seawater intrusion in coastal aquifers have been a problem of interest for many years and many methods have been developed. There is an extensive amount of literature on this subject. We refer to the books [1,2,3,4,5,6]. Two main models are often used to predict marine intrusion into coastal aquifers. The first one is the 2D sharp interface approach which assumes that the saltwater and the freshwater are immiscible separated by an abrupt interface [2], therefore the mixing zone is not taken into consideration. Secondly, the 3D variable density flow and solute transport approach [1] which considers that the two fluids are miscible and a transition zone caused by dispersion is then considered. This approach seems more realistic to simulate the seawater intrusion problem and to track the movement and changes of the transition zone between fresh and saltwater. The vertical section is then considered and the salt concentration is expressed in 3D. However, in the numerical context, this approach is costly (CPU time) compared to the 2D model with a sharp interface. In addition, we cannot determine exactly the transition area as in the sharp interface model. In many coastal aquifers, the thickness of the aquifer is negligible compared to its horizontal surface. The flow is supposed to be horizontal and the problem is reduced to a 2D model. In this case, the sharp interface approach can be adopted to study large scale seawater intrusion in such coastal aquifers. Thus, the extent of the transition (mixing) zone can be neglected when compared to the lateral dimensions of the aquifer. The two fluids can then be supposed separated by a sharp interface.
A review of the mathematical models for 2D saltwater intrusion in coastal aquifers can be viewed in [5,7,8]. Most of the previous numerical models used finite difference and finite element methods; see, for instance [9,10,11,12,13,14]. More recently, finite volume methods were developed and analyzed for the sharp interface model [15,16,17]. This approach leads to robust schemes applicable for unstructured grids and the approximate solution has various interesting properties which correspond to the properties of the physical solution. These methods have been useful for advective flow problems because they combine element by element conservation of mass with numerical stability and minimal numerical diffusion, see for instance [18].
The purposes of this paper are to derive a robust and accurate scheme on unstructured grids by using a fully implicit finite volume method for the coupled system modeling sharp interface seawater intrusion in coastal aquifers and present numerical simulations for complex case studies. This work aims to develop and implement a code coupling approach for this model in the framework of the parallel open-source platform DuMu X [19,20], based on the distributed and unified numerics environment (DUNE) [21], allowing simulations for large-scale field applications involving seawater intrusion in coastal aquifers. The model described above is built in DuMu X framework, which handles general input/output, memory management, grid generation, parallelism, etc. The code is an object-oriented software written in C++ and equipped with efficient solvers and massively parallel computation capability. Let us mention that the variable density flow and solute transport model is already implemented in DuMu X (module 1 p n c : one phase n components) which permits us to compare it to the sharp interface one.
This paper deals with a numerical method based on a finite volume scheme for solving the sharp interface model. More precisely, the time discretization is done by an implicit Euler method and in space we use a cell-centered finite volume method. The nonlinear system is solved by the Newton method and a biconjugate gradient stabilized (BiCGSTAB) method with incomplete LU factorisation (ILU) preconditioner is used to solve the linear systems. Numerical differentiation techniques are used to approximate the derivatives in the calculation of the Jacobian matrix. The control of the time-step is based on the number of iterations required by the Newton method to achieve convergence for the last time iteration. The time-step is reduced, if the number of iterations exceeds a specified threshold, whereas it is increased if the method converges within less iterations. This numerical scheme has been implemented and integrated in the platform DuMu X by the creation of a new module named two-phase saltwater intrusion ( 2 p S W I ). Our numerical model is verified with the field-scale free aquifer presented in [16]. We also apply the method to perform numerical simulations of the sharp interface approach in the Souss–Chtouka aquifer field case [22].
The outline of the paper is as follows. In Section 2, we briefly describe the governing equations for the sharp interface seawater intrusion model, while Section 3 describes the fully implicit finite volume scheme for solving the nonlinear coupled system. A description of the implementation of our scheme in the DuMu X framework is given in Section 4 and numerical results which illustrate the robustness of our approach are discussed for a field-scale free aquifer test and the Souss—Chtouka case study. Section 5 contains concluding remarks.

2. Mathematical Model

In this section, we recall the equations modeling a subsurface freshwater–saltwater flow in the case where saltwater and freshwater are immiscible separated by a sharp interface. Thus the coastal aquifer is divided into two regions. The freshwater flows in the upper part of a vertical section and the saltwater in the lower part, separated by a freshwater/saltwater interface, as shown in Figure 1.
The freshwater–saltwater model for groundwater flow with free aquifer consists of two vertically integrated governing equations, one for freshwater flow and the other for saltwater flow [2]. Let Ω be a bounded domain in R 2 representing the aquifer and [ 0 , T ] the time interval of interest. Invoking the continuity of flux and pressure on the interface and combining Darcy’s law with the continuity equation, the governing equations of the sharp interface freshwater-seawater model in coastal aquifers are expressed as follows (see, e.g., [1,2]):
Freshwater flow equation:
S f b f u , v + ϕ u t ϕ Z t div b f u , v K f u = Q f in Ω × 0 , T ,
Saltwater flow equation:
S s b s u , v v t + ϕ Z t div b s u , v K s v = Q s in Ω × 0 , T ,
where u and v are the hydraulic heads [m] of the freshwater and saltwater respectively and are the main unknowns of system (1) and (2). Z [ m ] is the saltwater front elevation given by:
Z = ( 1 + δ ) v δ u .
Q f and Q s indicate flows (fresh and salt) pumped or injected per unit surface of aquifer. The other different notations used in Equations (1) and (2) are exhibited in Table 1.
By replacing the functions Z, b f and b s by their formulas, system (1) and (2) becomes:
Freshwater flow equation:
( 1 + δ ) S f ( u v ) + ϕ u t ϕ ( 1 + δ ) v t div ( 1 + δ ) u v K f u = Q f ,
Saltwater flow equation:
S s ( 1 + δ ) v δ u Z B + ( 1 + δ ) ϕ v t δ ϕ u t div ( 1 + δ ) v δ u Z B K s v = Q s .
This system gives a 2 D description for tracking a sharp interface problem in a free aquifer. It represents two coupled parabolic partial differential equations that should be solved simultaneously for the freshwater head u and saltwater head v.

3. Numerical Scheme

The spatial discretization of the coupled system (3) and (4), subject to boundary and initial conditions, employs a conservative Finite Volume (FV) method. The time discretization is done by an implicit Euler method. For the sake of simplicity of exposition, here we present the scheme for a regular mesh. The extension to unstructured grids is straightforward. Here, we choose a cell-centered FV method. It consists in integrating Equations (3) and (4) on a control volume V k (Figure 2) and evaluating the fluxes at the interface γ k l between two adjacent elements V k and V l .
We denote by f k = 1 | V k | V k f d V the average of a function f on each element V k . By using the implicit Euler scheme for the time discretization and due to the fact that the approximation of the primary unknowns (u, v) and the physical parameters are constant on each element V k , the cell-centered FV schemes corresponding to the discretization of the freshwater and saltwater equations are given by:
Finite volume discretization of freshwater equation:
S f ( u k n + 1 v k n + 1 ) + ϕ k ( u k n + 1 u k n ) ϕ k ( v k n + 1 v k n ) Δ t n | V k | l V ( k ) | γ k l | ( u n + 1 v n + 1 ) K f k l u n + 1 k l · n k l = Δ t n ( 1 + δ ) Q f n + 1 k ;
Finite volume discretization of saltwater equation:
S s ( 1 + δ ) v k n + 1 δ u k n + 1 Z B + ( 1 + δ ) ϕ k ( v k n + 1 v k n ) δ ϕ k ( u k n + 1 u k n ) Δ t n | V k | l V ( k ) | γ k l | ( 1 + δ ) v n + 1 δ u n + 1 Z B K s k l v n + 1 k l · n k l = Δ t n Q s n + 1 k ;
where Δ t n is the time step, n k l denotes the unit outer normal to γ k l , V ( k ) is the set of adjacent elements of V k . The gradient operators u n + 1 k l and v n + 1 k l on the interfaces γ k l are calculated by a Two Point Flux Approximation (TPFA), see for instance [19,23]. A harmonic average of the values between two adjacent elements is used to calculate the diffusion coefficients at the interface γ k l .
By integrating the boundary and initial conditions into the FV discretization of system (5) and (6), we have two sets of nonlinear equations which are solved implicitly using Newthon’s method with variable time stepping. The control of the time-step is based on the number of iterations required by the Newton method to achieve convergence for the last time iteration. The time-step is reduced, if the number of iterations exceeds a specified threshold, whereas it is increased if the method converges within less iterations. Numerical differentiation techniques are used to approximate the derivatives in the calculation of the Jacobian matrix. This allows to transform the nonlinear system of equations for each iteration step into a linear system of equations. For solving the occurring linearized systems of equations, an iterative linear solver is used, namely, Bi-Conjugate Gradient STABilized (BiCGSTAB) method with ILU (Incomplete LU factorisation) preconditioner. This solver is integrated in the ISTL-Library of DUNE.

4. Numerical Simulations

In this section, we present a brief description of DuMu X and our module 2 p S W I developed and integrated in this plateform. The new algorithm is first validated on some classical tests. In cases we obtain a very good level of accuracy, showing also numerical convergence results; we furthermore confirm mass conservation up to machine precision. Then, to validate our approach, we consider two test cases. The first one is proposed in [16] for an unconfined coastal aquifer subject to pumping while the second one considers a realistic case with a highly complex geometry with different types of rocks [22].
Let us mention that throughout all numerical experiments, we observed that in no instance more than a maximum of 20 iterations was needed for the convergence of Newton’s method. Consequently, for this study the adopted strategy for the management of the time step is sufficient. However, various types of local time-stepping strategies have been proposed in the literature; see, for instance [24,25,26].
In explicit and semi-explicit schemes, the time step is restricted by the famous Courant–Friedrichs–Lewy (CFL) condition to ensure stability. This condition is usually very restrictive in reservoir-scale models, and it is therefore more common to use implicit schemes. Furthermore, sequential approaches introduce operator splitting errors and restrictions on the time step are mandatory to ensure mass conservation for instance. Implicit discretizations capable of taking large(r) time steps are therefore often preferred in practical computations. Although implicit schemes are more diffusive than their explicit counterparts, they yield better stability and mass conservation. Although this guarantees numerical stability in the solution, this does not guarantee nonlinear convergence. The computational performance of the code depends greatly on the convergence of Newton’s method and linear systems in an optimal number of iterations which is highly correlated to the choice of time step. The fully coupled fully implicit finite volume scheme, developed in this study, combined with time-accurate local time stepping allow Newton’s method and linear solver to converge within a reasonable number of iterations which saves computing time.
Finally a remarkable property of the scheme is that the discrete maximum principles (nonnegativity of the thickness of fresh and saltwater in the aquifer) is satisfied wich is crucial to obtain physically meaningful approximate solutions. This has been verified in all our simulations. A proof of this result for a simplified model could be find in [27]. The numerical analysis of the scheme including this property is out of the scope of this paper and will be considered in a future work.

4.1. DuMu X : Numerical Simulator

All our numerical developments have been implemented in DuMu X . It is a parallel free and open-source simulator for flow and transport processes in porous media. It is based on the Distributed and Unified Numerics Environment DUNE. It provides many tools to solve numerically partial differential equations and allowing, among other things, the management of mesh, discretization or linear and nonlinear solvers. The code is an object-oriented software written in C++ and has massively parallel computation capability. The modular concept of DuMu X makes it easy to integrate new modules adapted to our numerical scheme. For this, we have developed a new module named 2 p S W I in DuMu X . This module allows to numerically solve the coupled system (3) and (4) with a fully implicit scheme in time and a cell-centered finite volume method in space.

4.2. Numerical Tests

Our approach has been validated by solving several tests to approximate solutions of seawater flows with sharp interface in coastal aquifers. The first one performed is the rotating interface of Keulegan [28] who proposed an analytical solution that consists of describing the movement of the freshwater–saltwater interface without any external forcing in a confined aquifer. The numerical results are satisfactory and replicated to those in [29,30]. The results of these simulations are omitted since nothing startling was found. Instead, we concentrate on the results obtained in realistic two case studies in the next section. More precisely, we present two numerical studies, one proposed in [16] dealing with a hypothetical conceptual homogeneous unconfined aquifer subject to eleven scenarios with different pumping rates and different well locations, and the other focusing on a real field test case dealing with the contamination of the Souss–Chtouka aquifer. This second test is a coastal aquifer, with complex geometry, and geologic material heterogeneities and subject to a large stress period along with high pumping withdrawals.
All computations were performed on a laptop with Intel Xeon(R) CPU E3-1505M Processor (3.00 GHz) and 8 GB RAM. One of the objectives of this paper is to deliver computational performance also suitable for limited computational resources. Let us mention that in view of the CPU times required for the examples treated in this paper, all the simulations were performed sequentially. However, the new module developed can be used on multicore/multinode systems. The parallelization in DuMu X is carried out using the DUNE parallel library package based on MPI providing high parallel efficiency and allowing simulations with several tens of millions of degrees of freedom to be carried out, ideal for large-scale field applications. DuMu X has the ability to run on anything from single processor systems to highly parallel supercomputers with specialized hardware architectures.

4.3. Test 1: A Field-Scale Free Aquifer

The purpose of this test, considered in [16], is to see the influence of the pumping rate and well position on the displacement of the interface. As reported in [16], it aims to assess the validity of the sharp-interface approach for an unconfined coastal aquifer subjected to pumping by comparison with dispersive approach results. For this, numerical simulation for several scenarios have been achieved and compared to the results obtained in [16] where good agreement is observed.
The test considers a free aquifer of thickness 30 m and length 500 m. A Dirichlet condition for the saltwater head v = 30 m was imposed in x = 500 m which corresponds to the seaside. The constant freshwater flux 0.1 m 3 /day at land boundary ( x = 0 ) was considered and homogeneous Neumann boundary conditions are imposed on the rest of the border.
The total extraction rate in the pumping wells is assumed to be constant irrespective of the proportions of saltwater and freshwater: Q t = Q f + Q s in a similar manner described by [31]. The hydraulic head of freshwater and saltwater under steady-state conditions (before pumping) are used as initial conditions for the transient regime (after pumping) simulations. The simulation duration of the transient model is approximately 30 years. The parameters and properties of the aquifer are presented in Table 2.
Table 3 gives the different scenarios of the pumping rate and the wells position. Q t is the quantity of the water pumped, X w represents the distance between the well and the seaside ( x = 500 m) and Z w is the depth of the well. The well screen position is presented by a solid rectangular on Figure 3.
For all simulations, a uniform rectangular mesh of 250 × 10 cells is used for control volumes. Simulations were achieved with an initial time step of 0.01 s and a maximum time-step equal to 1 day have been considered. The tolerances for the Newton method and the BICGSTAB method are respectively 10 8 and 10 6 . In this case, Newton’s method converges rapidly in less than 5 iterations. The CPU time consumed for such simulation is less than 1 min.

4.3.1. Numerical Results by Varying the Pumping Rate

Here, we test the effect of the amount of water pumped on the evolution of the saltwater interface. We will simulate the latter by varying the pumping rate. We assume that the well is located at the point ( X w = 150 m, y = 0 m) in a depth Z w = 15 m. The different pumping scenarios correspond respectively to 0.1 m 3 /day (Sc-1), 0.05 m 3 /day (Sc-2), 0.07 m 3 /day (Sc-3) and 0.15 m 3 /day (Sc-4). The numerical results obtained in this case are presented in Figure 3. We observe that the salt front rises and the cone of the polluted water tends towards the pumping well. The progression of this local “upconing” is dependent on the amount of water pumped. In the second scenario, which corresponds to the smallest amount of water pumped, the interface is too far from the well base. However, we see that the interface touches the well bottom in the fourth scenario. The latter is salinized by the seawater, which can be dangerous for the operator. It can, therefore, be deduced that the salinization of the well is strongly related to the quantity of water pumped.

4.3.2. Numerical Results by Varying the Depth of the Well

Here, we focus on the effect of the well position combined with the pumping rate on the evolution of the sharp interface. The well is placed in the middle, bottom and top of the aquifer. Two different pumping rates are considered. Figure 4 shows the evolution of the sharp interface for Q t = 0.1 m 3 /day in the case of full (Sc-5) and partial (Sc-1, Sc-6) penetration of the pumping well.
In Figure 4 (Sc-5), the salt cone is gone and the pumping well is filled by the seawater. The risk of salinization decreases gradually as the well moves vertically up the aquifer.
We can make the same remarks in the case where we reduce the pumping rate to Q t = 0.07 m 3 /day (Figure 5). The effects are less comparable to those observed with Q t = 0.01 m 3 /day.

4.3.3. Numerical Results by Varying the Longitudinal Position of the Well

In order to visualize the impact of the well screen position on the evolution of the salt front, four longitudinal places are considered. The pumping rate is assumed to be constant ( Q t = 0.1 m 3 /day). The numerical results obtained in this case are presented in Figure 6.
Remark 1.
We have presented a comparison of our results versus those obtained in [16] for a test case dealing with a homogeneous unconfined aquifer subject to eleven scenarios with different pumping rates and different well locations. For this test case, in [16], it was also studied the effect of pumping rates and well screen location on the evolution of the salt front. They deduce that if the pumping wells are deep and close to the coast, the risk of salinization seawater is higher. It is also the case when the pumping rate is high. As shown above, the similarity between the two results validates our approach.

4.4. Test 2: Souss–Chtouka Aquifer Field Case

4.4.1. Geographic Location and Geologic Settings

The Souss—Chtouka aquifer is located at the southwestern part of Morocco in the Souss plain and its extension, the Chtouka plain. It is a hydrogeologic structure embraced between the alpine Haut-Atlas chain at the north, the eruptive massif of Siroua at the east, the Anti-Atlas massif at the south and is in contact with the Atlantic ocean at the west. Deep aquifers are identified in this hydrogeologic structure, in addition of the generalized phreatic aquifer, subject of this study. The phreatic aquifer of the Souss and Chtouka plains consists of heterogeneous fitting material of the valley. According to its geology, the deposits correspond to the quaternary alluvial sands and gravels of the Oued-Souss river, to the Moghrebian sandstones and coastal marine sands, to the Pliocene limestone with marly and conglomeratic intercalations of the downland areas of the Souss plain and to fluvial-lacustrine deposits of the Souss unit extending to the Anti-Atlas chain border.

4.4.2. Studied Domain and Discretization

The studied domain, of approximately 24 Km 2 , corresponds to the downland areas of the Souss plain and its extension, the Chtouka plain. It is located between the Haut-Atlas chain at the north, the 100 m contour of the 1968s piezometry at the east, the Oued-Massa river at the south and the Atlantic ocean at the west (Figure 7).
The geometry and boundaries of the aquifer are given in the left of Figure 8. For the mesh, we used an unstructured triangulation of 19,520 elements and 9871 vertexes as shown in the right of Figure 8.
In passing, we remark that the horizontal surface area of the aquifer is about 24 Km 2 , while the thickness of the aquifer is less than 1 Km (about 600 m in the north and 100 m in the south). In spite of the presence, locally, of relatively high depth regions in the northern part of the Souss–Chtouka coastal aquifer, the thickness of the latter remains small when compared to its lateral extent, which allows the Dupuit assumption to be valid. The fluid movement is therefore assumed to be horizontal.

4.4.3. Parameters and Boundary Conditions

Hydraulic conductivity of different geologic units that constitute the Souss–Chtouka aquifer, are determined following trial-and-error calibration operations in [22]. Ten zones have been recognized, for which the hydraulic conductivity values vary between 1.21 m/day and 40 m/day. The specific storage coefficient values varying between 10 5 m 1 and 4 × 10 5 m 1 and taken to be the same for fresh and salty phases are also specified in the right of Figure 8. Porosity values varying between 0.1 and 0.25 have been used, depending on the geologic material lithology. The specific storage coefficient distribution and the values of hydraulic conductivity and porosity are specified in Figure 9. Figure 10 shows the topography of the aquifer bottom. The Souss—Chtouka aquifer is fed, mainly, by the precipitation, the irrigation returns, the vertical leakance of the underlain Turonian limestone, the infiltration from the Oued-Souss river and the recharge from the Haut-Atlas chain at the north. In [22], the author delimited 8 zones in the study area with fluxes values varying between 1.05 × 10 6 m/day and 8.06 × 10 5 m/day.
Owing to the intensive exploitation of the Souss–Chtouka aquifer and the difficulty to have appropriate data, a general lowering of the water table is assumed. Many pumping wells given by the Moroccan authority of water (ONEP) are presented in Table 4. However, the numerical simulations obtained with these data show that the interface does not move for almost 80 years. In order to predict a significant displacement of the interface in the long term, we have multiplied the pumping rates given by the ONEP by a factor of 10.
To close the problem, boundary conditions have to be specified. Fixed head ( u = 100 m) is used on the upstream of the domain, at the east. At the northwest, at the contact with the Haut-Atlas chain and the south, fixed heads are also imposed representing the measured head [22]. On the western boundary, at the contact with the ocean a zero head is imposed ( u = v = 0 m). To obtain the initial conditions, the present model was run, under steady conditions, with the parameters resulted from calibration conducted in [22].
For all simulations, an unstructured mesh of 19,520 cells is used for control volumes. Simulations were achieved with an initial time step of 10 s. The tolerances for the Newton method and the BICGSTAB method are respectively 10 8 and 10 6 . In this case, Newton’s method converges in less than 20 iterations. As expected, the time-step is increased when Newton’s algorithm converges within less iterations. A remarkable attribute of the algorithm is that the total CPU time required for a 80 years simulation is less than 15 min on a laptop.
Let us end this section with the following remark. A second simulation for the Souss—Chtouka aquifer was performed with a refined mesh (37,320 cells and 18,861 vertex). The obtained results are very close to those of the previous coarse mesh. However, the CPU time is 15 min. In the following we will present results corresponding to the coarse mesh.

4.4.4. Numerical Results

We show the piezometry of the freshwater of the Souss–Chtouka aquifer between 1968 and 2048 in Figure 11.
For better visualization, all the curves are represented also on the bottom of the aquifer. To give a more realistic vision of the freshwater potential and to visualize the impact of the dramatic exploitation of the Souss–Chtouka Plain, we illustrate the piezometric contours of the freshwater head on the bottom of the aquifer (Figure 12).
Up to 1968, the aquifer was not subjected to any external force. We can see that the equipotentials are regular and vary from 0 m at the coast to 100 m at the landside. For the 2048 predictions, however, following a drastic exploitation scenario, consisting of a general lowering of the freshwater level at the upstream of the aquifer ( u = 100 0.625 t), along with over pumping rates, the water table decreased considerably, reflected in a generalized decrease from the coast to the entire basin. The level of freshwater decreases by 50 m in the upper part of the aquifer and an exceptional depression is located around the pumping wells position. Initially, the hydraulic head of saltwater is zero and the piezometric lines after 80 years are given in Figure 13.
Figure 14 illustrates the extension of the salt wedge, which is more prominent in the north towards the east (the depth of the reservoir reaches 650 m) and remains practically parallel to the western limit when going south. After 80 years of activity, the salt bevel has experienced a significant displacement in the north caused by intensive freshwater pumping, especially in regions where pumping wells are placed. However, in the southwestern region, the figure shows that the salt bevel is stable and does not advance in the continent (regions with low permeabilities).
To show the extent of the salt intrusion according to the location, we chose five sections perpendicular to the Atlantic coast and oriented from the East to the West. The profiles of u, v and Z according to different sections are presented in Figure 15, Figure 16 and Figure 17.
The first remark we can make is the progress of the freshwater/saltwater interface for the great depths to the east in the Agadir region, over a distance estimated at 7000 m on the substratum at a depth of 600 m. Indeed, in the area concerned, the high value of the permeability of the formations contributes to the advance of the bevel. Significant vertical advancement of the interface is noticed in the first section and an “upconing” is developed under a large flow well. The pumping well is close to the coast, which accelerates the advancement of the salt bevel during the 80 years of activity (Figure 15).
On the second section, (Figure 16 (left)), we can see a significant displacement of the salt bevel laterally around 4 km and especially vertically, so that an "upconing" developed below a large flow of three pumping wells. A slight cone of depression can also be noted at the location of these wells.
Along the third section, (Figure 16 (right)), the salt bevel is displaced laterally over a distance of at least 6 km after 80 years of activity and an “upconing” is developed below the pumping well. The pumping effects are less visible compared to the first two sections since the corresponding flow rate is low. On the other hand, a strong cone of depression developed at the well site. The regular shape of the bedrock and the high permeability values in this region contributed to the lateral advancement of the salt bevel.
On the fourth section, (Figure 17 (left)), there is always a sustained movement of the salt bevel over the years to reach its maximum value after 80 years of service. Finally, the fifth section, (Figure 17 (right)) shows almost no movement of the salt bevel despite a linear decrease in the free surface area, which is generalized to the entire slick. It should be noted that the region interested in logging is free of any pumping zones.
Remark 2.
The numerical results of the second test showed that the new developed module is capable of providing a stable solution and can predict the location, the shape and the extent of the water table and the freshwater–saltwater interface in coastal aquifers under various stress conditions, demonstrating its robustness with satisfactory rate of convergence.
Let us end this section by the following remark. For numerous tests, the obtained results are satisfactory and the numerical computations for the coupled system have demonstrated that this approach yields physically realistic flow fields in highly heterogeneous fields. Furthermore, the fully coupled fully implicit scheme greatly reduce the runtime of the simulations.

5. Conclusions

In this paper, we have presented a fully implicit finite volume method for solving a seawater intrusion problem in coastal aquifers. We have shown that this method is efficient, and accurate, and is capable of solving seawater flow problems into complex heterogeneous aquifers. The proposed approach was implemented in the framework of the parallel open-source platform DuMu X . Numerical results on a realistic heterogeneous aquifer were presented. The use of an existing and well-established framework such as DuMu X to implement the simulations has presented several advantages. First, the framework provided most of the basic numerical tools for implementing the new methods. Then, the structure of the framework means that extensions to quite varied geometrical and physical situations will be reasonably straightforward. Its open architecture facilitates further development for specific needs. Future works will focus on the numerical analysis of the scheme and the extension of the methodology to more advanced models, like the mixte sharp diffuse interface model recently introduced in [32,33,34].

Author Contributions

Data curation, A.A.; Formal analysis, B.A., and M.E.O.; Investigation, A.A., B.A., M.E.O. and K.T; Software, M.E.O. and K.T.; Supervision, B.A. and M.E.O.; Writing—review & editing, A.A., B.A., M.E.O. and K.T. All authors have read and agreed to the published version of the manuscript.

Funding

The publication of this article was funded by ADERA, Aquitaine Region, France.

Acknowledgments

This work was partially supported by the International Research Group of CNRS GDRI LEM2I (Laboratoire Euro-Maghrébin de Mathématiques et de leurs Interactions). The financial supports are gratefully acknowledged. The authors gratefully thank the anonymous referees for their insightful comments and suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bear, J.; Cheng, A.H.D. Modeling Groundwater Flow and Contaminant Transport; Springer: Dordrecht, The Netherlands, 2000. [Google Scholar]
  2. Bear, J.; Cheng, A.H.D.; Sorek, S.; Ouazar, D.; Herrera, I. Seawater Intrusion in Coastal Aquifers, Concepts, Methods and Practices; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1999. [Google Scholar]
  3. Bear, J.; Verruijt, A. Modelling Groundwater Flow and Pollution; D. Reidel Publishing Company: Dordrecht, The Netherlands, 1987. [Google Scholar]
  4. Cheng, A.H.D. Multilayered Aquifier Systems. In Fundamentals and Applications; CRC Press Published: New York, NY, USA, 2000. [Google Scholar]
  5. Cheng, A.H.D.; Ouazar, D. Coastal Aquifer Management-Monitoring, Modeling, and Case Studies; CRC Press: New York, NY, USA; Washington, DC, USA, 2003. [Google Scholar]
  6. Diersch, H.-J.G.G. FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media; Springer: Berlin/Heidelberg, Germany, 2014. [Google Scholar]
  7. Bobba, G.H. Mathematical models for saltwater intrusion in coastal aquifers. Water Resour. Manag. 1993, 7, 3–37. [Google Scholar] [CrossRef]
  8. Werner, A.D.; Bakke, M.; Post, V.E.A.; Vandenbohede, A.; Barry, D.A. Seawater intrusion processes, investigation and management: Recent advances and future challenges. Adv. Water Resour. 2013, 51, 3–26. [Google Scholar] [CrossRef]
  9. Cobaner, M.; Yurtal, R.; Dogan, A.; Motz, L.H. Three dimensional simulation of seawater intrusion in coastal aquifers: A case study in the Goksu Deltaic Plain. J. Hydrol. 2012, 464–465, 262–280. [Google Scholar] [CrossRef]
  10. Essaid, H.L. A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: Model developpement and application. Water Resour. Res. 1990, 26, 1431–1454. [Google Scholar] [CrossRef]
  11. Fahs, M.; Koohbor, B.; Belfort, B.; Ataie-Ashtiani, B.; Simmons, C.T.; Younes, A.; Ackerer, P. A generalized semi-analytical solution for the dispersive Henry problem: Effect of stratification and anisotropy on seawater intrusion. Water 2018, 10, 230. [Google Scholar] [CrossRef] [Green Version]
  12. Houben, G.H.; Stoeckl, L.; Mariner, K.E.; Choudhury, A.S. The influence of heterogeneity on coastal groundwater flow-physical and numerical modeling of fringing reefs, dykes and structured conductivity fields. Adv. Water Resour. 2017, 113, 155–166. [Google Scholar] [CrossRef]
  13. Huyakorn, P.S.; Wu, Y.S.; Park, N.S. Multiphase approach to the numerical solution of a sharp-interface saltwater intrusion problem. Water Resour. Res. 1996, 32, 93–102. [Google Scholar] [CrossRef]
  14. Povich, T.J.; Dawson, C.N.; Farthing, M.W.; Kees, C.E. Finite element methods for variable density flow and solute transport. Comput. Geosci. 2013, 17, 529–549. [Google Scholar] [CrossRef]
  15. Bouzouf, B.; Chen, Z. A comparison of finite volume method and sharp model for two dimensional saltwater intrusion modeling. Can. J. Civ. Eng. 2014, 41, 191–196. [Google Scholar] [CrossRef]
  16. Mehdizadeh, S.; Vafaie, F.; Abolghasemi, H. Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping. Environ. Earth Sci. 2015, 73, 8345–8355. [Google Scholar] [CrossRef]
  17. Vafaie, F.; Mehdizadeh, S. Investigation of sea level rise effect on saltwater intrusion in an unconfined coastal aquifer using sharp-interface approach. Int. J. Glob. Warm. 2015, 8, 501–515. [Google Scholar] [CrossRef]
  18. Eymard, R.; Gallouët, T.; Herbin, R. Finite volume methods. In Handbook of Numerical Analysis; Ciarlet, P.G., Lions, J.L., Eds.; North-Holland: Amsterdam, The Netherlands, 2000; pp. 713–1022. [Google Scholar]
  19. DuMuX DUNE for Multi-{Phase, Component, Scale, Physics, ...} Flow and Transport in Porous Media. Available online: http://www.dumux.org (accessed on 1 June 2020).
  20. Flemisch, B.; Darcis, M.; Erbertseder, K.; Faigle, B.; Lauser, A.; Mosthaf, V.K.; Muthing, S.; Nuske, P.; Tatomir, A.; Wolf, M.; et al. DuMuX: Dune for multi-{Phase, Component, Scale, Physics, …} flow and transport in porous media. Adv. Water Resour. 2011, 34, 1102–1112. [Google Scholar] [CrossRef]
  21. DUNE, the Distributed and Unified Numerics Environment. Available online: https://www.dune-project.org (accessed on 1 June 2020).
  22. Abouelmahassine, N. Modélisation Préliminaire du Biseau salé par Éléments Finis dans la Zone Côtière de Souss–Chtouka. Master’s Thesis, Ecole Mohammadia d’Ingénieurs, Rabat, Morocco, 2002. [Google Scholar]
  23. Eymard, R.; Gallouët, T.; Guichard, C.; Herbin, R.; Masson, R. TP or not TP, that is the question. Comput. Geosci. 2014, 18, 285–296. [Google Scholar] [CrossRef] [Green Version]
  24. Linga, G.; Møyner, O.; Nilsen, H.M.; Moncorgé, A.; Lie, K.A. An implicit local time-stepping method based on cell reordering for multiphase flow in porous media. J. Comput. Phys. 2020, 100051. [Google Scholar] [CrossRef]
  25. Vohralik, M.; Wheeler, M.F. A posteriori error estimates, stopping criteria, and adaptivity for two-phase flows. Comput. Geosci. 2013, 17, 789–812. [Google Scholar] [CrossRef] [Green Version]
  26. Younes, A.; Ackerer, P. Empirical versus time stepping with embedded error control for density-driven flow in porous media. Water Resour. Res. 2010, 46. [Google Scholar] [CrossRef]
  27. Ait Hammou Oulhaj, A. Numerical analysis of a finite volume scheme for a seawater intrusion model with cross-diffusion in an unconfined aquifer. Numer. Methods Partial Differ. Equ. 2018, 34, 857–880. [Google Scholar] [CrossRef] [Green Version]
  28. Keulegan, H.G. An Example Report on Model Laws for Density Current; U.S. National Bureau of Standards: Gaithersburg, MD, USA, 1954.
  29. Abudawia, A.; Rosier, C. Numerical analysis for a seawater intrusion problem in a confined aquifer. Math. Comput. Simul. 2015, 118, 2–16. [Google Scholar] [CrossRef]
  30. Marion, P.; Najib, K.; Rosier, C. Numerical simulations for a seawater intrusion problem in a free aquifer. Appl. Numer. Math. 2014, 75, 48–60. [Google Scholar] [CrossRef]
  31. Shi, L.; Cui, L.; Park, N.; Huyakorn, P.S. Applicability of a sharp-interface model for estimating steady-state salinity at pumping wells-validation against sand-tank experiments. Contam. Hydrol. 2011, 124, 35–42. [Google Scholar] [CrossRef]
  32. Choquet, C.; Diédhiou, M.; Rosier, C. Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer, Numerical simulations. SIAM J. Appl. Math. 2016, 76, 138–158. [Google Scholar] [CrossRef] [Green Version]
  33. Cherfils, L.; Choquet, C.; Diédhiou, M. Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows. Math. Comput. Simul. 2017, 137, 246–265. [Google Scholar] [CrossRef]
  34. Abudawia, A.; Mourad, A.; Rodrigues, J.H.; Rosier, C. A finite element method for a seawater intrusion problem in unconfined aquifers. Appl. Numer. Math. 2018, 127, 349–369. [Google Scholar] [CrossRef]
Figure 1. Interface and other related elevation from [16].
Figure 1. Interface and other related elevation from [16].
Water 12 01639 g001
Figure 2. Discretization by the cell-centered finite volume method.
Figure 2. Discretization by the cell-centered finite volume method.
Water 12 01639 g002
Figure 3. Evolution of the salt front by varying the pumping rate and fixing the position of the well.
Figure 3. Evolution of the salt front by varying the pumping rate and fixing the position of the well.
Water 12 01639 g003
Figure 4. Evolution of the salt front by varying the depth of the well and fixing the pumping rate at 0.1 m 3 /day.
Figure 4. Evolution of the salt front by varying the depth of the well and fixing the pumping rate at 0.1 m 3 /day.
Water 12 01639 g004
Figure 5. Evolution of the salt front by varying the depth of the well and fixing the pumping rate at 0.07 m 3 /day.
Figure 5. Evolution of the salt front by varying the depth of the well and fixing the pumping rate at 0.07 m 3 /day.
Water 12 01639 g005
Figure 6. Evolution of the salt front for the Sc-1, Sc-7, Sc-8 and Sc-9 scenarios.
Figure 6. Evolution of the salt front for the Sc-1, Sc-7, Sc-8 and Sc-9 scenarios.
Water 12 01639 g006
Figure 7. Study area and geolocation of the aquifer.
Figure 7. Study area and geolocation of the aquifer.
Water 12 01639 g007
Figure 8. Left: Geometry and boundaries of the aquifer and well locations. Right: mesh for the aquifer.
Figure 8. Left: Geometry and boundaries of the aquifer and well locations. Right: mesh for the aquifer.
Water 12 01639 g008
Figure 9. Distribution maps of specific storage coefficient, hydraulic conductivity and porosity, respectively, of the aquifer.
Figure 9. Distribution maps of specific storage coefficient, hydraulic conductivity and porosity, respectively, of the aquifer.
Water 12 01639 g009
Figure 10. Bottom elevation contours of the aquifer.
Figure 10. Bottom elevation contours of the aquifer.
Water 12 01639 g010
Figure 11. Plan view of the piezometry of the plain Souss–Chtouka before (left) and after (right) solicitation.
Figure 11. Plan view of the piezometry of the plain Souss–Chtouka before (left) and after (right) solicitation.
Water 12 01639 g011
Figure 12. Contours representing the freshwater potential of the Souss–Chtouka Plain illustrated on the bottom of the aquifer in 1968 (left) and 2048 (right).
Figure 12. Contours representing the freshwater potential of the Souss–Chtouka Plain illustrated on the bottom of the aquifer in 1968 (left) and 2048 (right).
Water 12 01639 g012
Figure 13. Piezometric contours for saltwater hydraulic heads in 2048.
Figure 13. Piezometric contours for saltwater hydraulic heads in 2048.
Water 12 01639 g013
Figure 14. Position of the saltwater/freshwater interface in 1968 (black) and in 2048 (red).
Figure 14. Position of the saltwater/freshwater interface in 1968 (black) and in 2048 (red).
Water 12 01639 g014
Figure 15. Left: Plan view of level positions. Right: vertical cross section showing the interface position (red line), the free surface position (blue line) and the variation of the saltwater hydraulic head (green line) for different times: cut on the level 1.
Figure 15. Left: Plan view of level positions. Right: vertical cross section showing the interface position (red line), the free surface position (blue line) and the variation of the saltwater hydraulic head (green line) for different times: cut on the level 1.
Water 12 01639 g015
Figure 16. Vertical cross section showing the interface position (red line), the free surface position (blue line) and the variation of the saltwater hydraulic head (green line), for different times: cut on the levels 2 (left) and 3 (right).
Figure 16. Vertical cross section showing the interface position (red line), the free surface position (blue line) and the variation of the saltwater hydraulic head (green line), for different times: cut on the levels 2 (left) and 3 (right).
Water 12 01639 g016
Figure 17. Vertical cross section showing the interface position (red line), the free surface position (blue line) and the variation of the saltwater hydraulic head (green line), for different times: cut on the levels 4 (left) and 5 (right).
Figure 17. Vertical cross section showing the interface position (red line), the free surface position (blue line) and the variation of the saltwater hydraulic head (green line), for different times: cut on the levels 4 (left) and 5 (right).
Water 12 01639 g017
Table 1. Notations of parameters and functions used in system (1) and (2).
Table 1. Notations of parameters and functions used in system (1) and (2).
Parameters
Z T : water table elevation [m] Z B : bottom of the aquifer [m]
K f , K s : hydraulic conductivities [m/day] S f , S s : specific storage coefficients [1/m]
ρ f , ρ s : densities of fresh and saltwater [kg/m 3 ] Q f , Q s : flows [m/day]
ϕ : porosity of medium [%] δ = ρ f ρ s ρ f
Functions
b f = u Z : freshwater thickness [m] b s = Z Z B : saltwater thickness [m]
Table 2. Parameters and properties of the aquifer.
Table 2. Parameters and properties of the aquifer.
Parameters δ K f K s ϕ ρ f ρ s S f S s
Values40 1.0 41 0.35 10001025 0.0 0.0
Table 3. Parameters for different scenarios.
Table 3. Parameters for different scenarios.
ScenariosSc-1Sc-2Sc-3Sc-4Sc-5Sc-6Sc-7Sc-8Sc-9Sc-10Sc-11
Q t [m 3 /d]0.10.050.070.150.10.10.10.10.10.070.07
X w [m]150150150150150150200300100150150
Z w [m]15151515025151515025
Table 4. Well position and rate pumping.
Table 4. Well position and rate pumping.
Wells pumping P 1 P 2 P 3 P 4 P 5
x [m]99,811100,277101,250100,750102,000
y [m]384,285383,966375,920374,280375,400
Rate [m 3 /day]−2918.84−2686.99−1011.10−673.31−710.22
Wells pumping P 6 P 7 P 8 P 9
x [m]99,507.7101,710100,0621,000,000
y [m]375,246374,431372,311374,996
Rate [m 3 /day]−1113.03−1533.89−1402.35−1488.54

Share and Cite

MDPI and ACS Style

Aharmouch, A.; Amaziane, B.; El Ossmani, M.; Talali, K. A Fully Implicit Finite Volume Scheme for a Seawater Intrusion Problem in Coastal Aquifers. Water 2020, 12, 1639. https://doi.org/10.3390/w12061639

AMA Style

Aharmouch A, Amaziane B, El Ossmani M, Talali K. A Fully Implicit Finite Volume Scheme for a Seawater Intrusion Problem in Coastal Aquifers. Water. 2020; 12(6):1639. https://doi.org/10.3390/w12061639

Chicago/Turabian Style

Aharmouch, Abdelkrim, Brahim Amaziane, Mustapha El Ossmani, and Khadija Talali. 2020. "A Fully Implicit Finite Volume Scheme for a Seawater Intrusion Problem in Coastal Aquifers" Water 12, no. 6: 1639. https://doi.org/10.3390/w12061639

APA Style

Aharmouch, A., Amaziane, B., El Ossmani, M., & Talali, K. (2020). A Fully Implicit Finite Volume Scheme for a Seawater Intrusion Problem in Coastal Aquifers. Water, 12(6), 1639. https://doi.org/10.3390/w12061639

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop