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

: We present a numerical framework for efﬁciently simulating seawater ﬂow in coastal aquifers using a ﬁnite volume method. The mathematical model consists of coupled and nonlinear partial differential equations. Difﬁculties arise from the nonlinear structure of the system and the complexity of natural ﬁelds, 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 ﬂow dynamics at each time step, is needed to overcome the stability problem of the time integration and mass conservation. A fully implicit ﬁnite 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 ﬂow. Lastly, numerical results of a realistic test case are presented to prove the efﬁciency and the performance of the method.


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

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: Saltwater flow equation: 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: 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: Saltwater flow equation: This system gives a 2D 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.

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 γ kl between two adjacent elements V k and V l . We denote by f k = 1 |V k | V k f dV 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: Finite volume discretization of saltwater equation: where ∆t n is the time step, n kl denotes the unit outer normal to γ kl , V(k) is the set of adjacent elements of V k . The gradient operators ∇u n+1 kl and ∇v n+1 kl on the interfaces γ kl 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 γ kl . 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.

Numerical Simulations
In this section, we present a brief description of DuMu X and our module 2p − SW 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.

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 2p − SW 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.

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.

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.

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.

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.

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.

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.

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.

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.

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 Figures 15-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.

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].