M Id Moulay

S Hamidi

S Eddamiri

Y Zehraoui

S Raaidi

Mathematical modeling of seawater intrusion in coastal aquifers

Keywords: Saltwater intrusion, finite volume, TPFA, MPFA, DuMu X, numerical simulation, numerical analysis, convergence, heterogeneity, porous media. Résumé Intrusion marine, volumes finis, TPFA, MPFA, DuMu X, simulation numérique, analyse numérique, convergence, hétérogénéité, milieux poreux. Chapter 1. General introduction

Je remercie sincèrement les membres du Jury : Mesdames Z. Mghazli et C. Rosier et Monsieur A. Taakili d'avoir accepté d'être mes rapporteurs. Je leur exprime toute ma gratitude pour leurs remarques pertinentes et constructives qui m'ont permis d'améliorer ce document. J'exprime ma gratitude à Z. Mghazli pour m'avoir fait l'honneur de présider ce Jury. Des remerciements particuliers sont adressés à Messieurs M. Afif et A. Ezziani, qui ont gentiment accepté d'être membres de mon Jury. J'exprime toute ma reconnaissance à Monsieur A. Aharmouch pour ses discussion fructueuses, sa collaboration et pour sa participation au Jury.

Je tiens à remercier tout particulièrement mes directeurs de thèse, Messieurs M. El Ossmani et B. Amaziane, pour m'avoir proposé ce sujet et pour m'avoir accompagné dans ma formation de Mathématiques Appliquées. Tout au long de ce travail, j'ai bénéficié de leurs compétences scientifiques en mathématiques appliquées, notamment en analyse numérique, en programmation C++ et en simulation numérique. Leur profonde humanité m'a beaucoup encouragée et m'a permis de préparer cette thèse dans les meilleures conditions possibles. Je leur suis très reconnaissante de m'avoir permis de participer à des conférences internationales. Je n'oublierai pas non plus la gentillesse et le soutien de Monsieur E. Ahusborde pour son aide pour l'initiation et durant le développement des codes dans la plateforme DuMuX.

Introduction

In this thesis, we are interested in the mathematical modeling and numerical simulation of seawater intrusion into coastal aquifers. Thus, this chapter presents a general introduction of this manuscript while outlining the objective, the structure and the different developments and improvements made in the context of this work.

We introduce the motivation and the problematic of this research in section 1.2. Due to accelerated urbanization and the population explosion and industrial activities, the exploitation of water resources becomes more and more considerable, particularly in coastal regions which are pole of attraction. Therefore, the potential of freshwater decreases and the seawater advances towards the land-side which degrade the quality of freshwater. A good management of these areas requires the use of appropriate strategies and efficient tools for the simulation and prediction of seawater intrusion in theses zones. Different approaches used in the modeling of marine intrusion are discussed and presented. After that, we present in section 1.3 the main results achieved in this thesis. In this work, two marine intrusion modules, based on recent and efficient tools, are developed and integrated in the platform DuMu X . A brief description of theses modules is presented and the different developments done in the implementation of the finite volume numerical methods is explained. The structure of this manuscript is explained in section 1.4.

Motivation and problem setting

Groundwater is an essential source of water that represents only 3% of the world total water. Most groundwater is in contact with saltwater. These valuable resources suffer from drought and contamination by seawater. Freshwater extraction is increasing with population growth and industrial activity. The water potential of the aquifers is threatened, both qualitatively and quantitatively, by the significant decline in peizometric levels and ocean intrusion due to high pumping rates. A good management of theses resources is becoming an important issue and a major challenge, see for instance [92,[START_REF] Werner | Seawater intrusion processes, investigation and management: recent advances and future challenges[END_REF].

The coastal regions are often, poles of attraction for industrial, agricultural and tourist activities. The large increase in demand for drinking water due to population growth and the expansion of economic activities in these areas presents a significant challenge to the hydraulic potential of coastal aquifers. Over-exploitation and excessive pumping of these resources lead to the invasion of aquifers by salt water. As a consequence, seawater moves towards the landside under freshwater because of its important density compared to freshwater. Indeed, over-pumping from wells located in sensitive areas of the water table can cause the advancement of the saltwater wedge and thus cause the degradation of the quality of the water in the coastal water table. The prediction of the hydrodynamics of freshwater and saltwater and the implementation of efficient tools is essential for a good management of these resources in the coastal zones, because once contaminated, the freshwater hardly recovers its state before the pollution, hence the importance of prediction.

The numerical modeling and simulation of marine intrusion in coastal aquifers has been a significant research topic and different approaches have been developed to predict seawater intrusion. Two classical approaches are often employed to tackle marine intrusion problems. We refer to the books [START_REF] Bear | Modelling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Bear | Seawater Intrusion in Coastal Aquifers, Concepts, Methods and Practices[END_REF][START_REF] Bear | Modelling Groundwater Flow and Pollution[END_REF][START_REF] Cheng | Multilayered Aquifier Systems[END_REF][START_REF] Cheng | Coastal Aquifer Management-Monitoring, Modeling, and Case Studies[END_REF][START_REF] Diersch | FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media[END_REF].

The first one is the 2D sharp interface approach [START_REF] Bear | Seawater Intrusion in Coastal Aquifers, Concepts, Methods and Practices[END_REF] which assumed that the freshwater 1.3. The objective of the thesis and saltwater fluids do not interact with each other and are separated by an abrupt interface. In this case, flow at the saltwater and freshwater zones is considered and no mass transfer occurs at the transition zones. In natural conditions where external forces are neglected, saltwater is supposed to be immobile; and only the freshwater movement is then considered. In such case, we are simulating only freshwater flow.

Secondly, the 3D variable density flow and solute transport approach [START_REF] Bear | Modelling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Diersch | FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media[END_REF] which considers that the two fluids are miscible. Two processes then govern the transport of the dissolved salt, namely convection and dispersion. This approach seems more credible to track movement and changes in the freshwater/saltwater transition zone, but no interface is defined to determine the area of saltwater wedge occupation or that of the water table which may be moved by pumping. In this case, the vertical variation of the fluids is considered. Furthermore, in terms of CPU time, this approach (3D model) is expensive compared to the sharp interface approach and produces significant numerical diffusion in the transition zones between freshwater and saltwater, and saturated and unsaturated zones.

In the coastal areas, the thickness of the aquifer is insignificant compared to its horizontal thickness. The problem is simplified and reduced to a 2D model since the fluids flow is assumed to be horizontal. In such case, the sharp interface approach can be adopted and extended to deal with large scale aquifers and to provide physically efficient results. In this context, a new approach (2D model) coupling the physical sides of the 3D variable density model and the simplicity of the sharp interface approch is introduced in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. The dynamic of the diffusive zones is tacking into account.

The objective of the thesis

The purpose of this research is the development and implementation of finite volume numerical schemes for seawater intrusion models into coastal aquifers. This study is based mainly on the of an efficient and robust numerical framework simulating seawater intrusion models and able to take into account various hydraulic and hydrogeological phenomena. Furthermore, we are interested in the numerical analysis of a conservative finite volume method for the flow of two immiscible fluids in an anisotropic and heterogeneous porous medium.

In this section, we describe the different development done in the implementation of our finite volume seawater intrusion numerical schemes. We have adopted a fully coupled, fully implicit approach for the resolution of freshwater and saltwater flows into coastal aquifers. We employ an Euler implicit scheme for the discretisation in time and a Cell Centered finite volume method to discretize in space. We linearise the coupled systems using Newton's method in each time step wherein each iteration, a linear system is solved using iterative Krylov methods. We use mainly BiCGSTAB (Biconjugate gradient stabilized method) solvers combined with a preconditioner (ILU) and GMRES (Generalized minimal residual method) solver coupled with a preconditioned (AMG Backend). Numerical differentiation techniques are used to approximate the derivatives in the calculation of the Jacobian matrix.

We have developed and implemented two modules in the the framework of the parallel open source DuMu X , based on the Distributed and Unified Numerics Environment (DUNE), allowing simulations for large-scale field applications involving seawater intrusion in coastal aquifers. The two modules called respectively "2p-SWI" (for the sharp interface approach) and "2pdiff-SWI" (for the sharp-diffuse interfaces approach). The object-oriented code provides recent and efficient numerical tools and allows us to predict marine intrusion while taking into account various hydraulic and hydrogeological phenomena. All our recent developments have been implemented in DuMu X versions 2.12 and 3.0.

The sharp interface module

The 2p-SWI module is based on a TPFA (See [START_REF] Droniou | Finite volume schemes for diffusion equations: Introduction to and review of modern methods[END_REF][START_REF] Eymard | Finite volume methods[END_REF]) numerical finite volume scheme by means of the sharp interface approach. The validation of the integrated module is performed on different 1D, 2D test cases, in particular the Souss-Chtouka plain in Morocco. In the case of a confined aquifer, we have compared our numerical results with an analytical solution proposed in [START_REF] Keulegan | An example report on model laws for density current[END_REF]. The numerical results feet well with those presented in [START_REF] Abudawia | Numerical analysis for a seawater intrusion problem in a confined aquifer[END_REF][START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF]. The efficiency and accuracy of our module in the case of a free aquifer are investigated by performing the numerical convergence of the scheme. Afterwards, we have extent our module to deal with heterogenous and realistic test cases. We have applied our method to two real test cases: "the field scale free aquifer" presented in [START_REF] Mehdizadeh | Vertical leakage in sharp-interface seawater intrusion models of layered coastal aquifers[END_REF], the Souss-Chtouka plain [5] and the Tripoli aquifer test case [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF]. Numerical simulations for the long-term demonstrate the applicability of our developed approach in highly heterogeneous coastal aquifers, and proved that our module is efficient and able to solve numerically 2D seawater intrusion problem using the sharp interface approach. Chapter 3 summarizes all the numerical illustrations of these test cases.

The extension of the methodology to more advanced models, like the mixte sharp-diffuse interfaces model recently introduced in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] will be the subject of the next section.

The sharp-diffuse interfaces module

Firstly, we have extended the sharp interface module to more general module which take into account the width of the transition zones. It is based on a TPFA finite volume numerical method. This module provides more informations about the position of the diffusion zones and yields physically realist results. To ensure the validity of our module, we proceed to the study of the numerical convergence of the scheme and compared our numerical results with an analytical solution presented in [START_REF] Ferris | Cyclic fluctuations of water level as a basis for determining aquifer transmissibility[END_REF]. The ability of our 2D module to provide an efficient and accurate results is checked by comparing the numerical results with the classical 3D model for miscible displacements.

After that, we have developed the "2pdiff-SWI" module to deal with the anisotropy and the heterogeneity of the coastal aquifer and to provide accurate reults. The developed module uses an MPFA [START_REF] Droniou | Finite volume schemes for diffusion equations: Introduction to and review of modern methods[END_REF] finite volume numerical method, where the subfluxes are approximated on the half-edges of the control volume. In fact, this method tackles the non-linearity of the problem, the anisotropy and the heterogeneity of the aquifer on unstructured mesh, which improves the efficiency and accuracy of the numerical computations of the system of equations. Numerical results for a test case showed the impact of various hydraulic and hydrogeological phenomena on the displacement of the freshwater and saltwater flows.

1.4. Thesis structure 5

Thesis structure

The purpose of this section is to present the general structure of the thesis, which is divided into six chapters.

Chapter 2 is devoted to the mathematical modeling of seawater intrusion in coastal aquifers. First, we present a state of the art on the existing approaches modeling marine intrusion. After that, we describe the derivation of the governing equations modeling saltwater and freshwater flows using the sharp interface model. The set of physical assumptions used in the formulation of the mathematical model is presented.

In chapter 3, we develop a fully coupled, fully implicit finite volume method for the numerical resolution of the sharp interface model. To discretize in time, we use an implicit Euler method and a TPFA finite volume method (See [START_REF] Droniou | Finite volume schemes for diffusion equations: Introduction to and review of modern methods[END_REF][START_REF] Eymard | Finite volume methods[END_REF]) in space. The coupled system is linearized using Newton's method, wherein each iteration a linear system is solved. In this context, we have integrated "2p-SWI" module in the framework of DuMu X . A brief description of the platform is presented. Afterwards, we illustrate the numerical results obtained in the validation of our module. The module is first tested on a reference case concerning the evolution of the interface under density contrast in a confined aquifer. The efficiency of the module is examined by conducting the numerical scheme study in the case of a free aquifer. We applied the method to different test cases presented in [START_REF] Abudawia | Numerical analysis for a seawater intrusion problem in a confined aquifer[END_REF][START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF][START_REF] Mehdizadeh | Vertical leakage in sharp-interface seawater intrusion models of layered coastal aquifers[END_REF]. The results obtained are in good agreement, which guarantees the validity of our implemented module. Moreover, we test our module in a two realistic cases corresponding to the Souss-Chtouka plain [5] in Morocco and Tripoli aquifer [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF] in Lebanon. The numerical model is based on available hydrogeological data in real scale by simulating the progress of seawater intrusion for several years. Numerical results have proved the ability of our developed module to predict water table and salt front elevations under various stress conditions and its ability to yield physically accurate and efficient solutions. Part of the results of this chapter have been published in [4].

Chapter 4 aims to develop a finite volume method for the numerical simulation of the sharp-diffuse interfaces model. First, we briefly outline the derivation of the mathematical model. Then, we present different numerical results used for the validation of our module. First, we proceed to the study of the numerical convergence of the scheme. After that, we consider a test case which concerns a free aquifer subjected to tidal fluctuations where the analytical solution is known. The numerical results are compared to the analytical solution and proved a good agreement with those presented in [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF][START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. Both test cases show that our fully implicit approach is suitable to simulate a seawater intrusion problem taking into account the dynamics of transition zones. In addition, we propose to compare the numerical results of the 2pdiff-SWI module (2D model) with those provided by the classical 3D model for miscible displacements. The numerical results have shown that this approach yields physically realistic and performant results.

Chapter 5 is devoted to the study of the convergence of an MPFA finite volume scheme for the sharp-diffuse interfaces model in an anisotropic and heterogeneous free aquifer. We develop a fully coupled, fully implicit approach finite volume method to discretize the coupled system on an unstructured mesh. The numerical scheme ensures a discrete maximum principle taking into account the heterogeneity and anisotropy of the aquifer. Using a fixed point theorem, we have established the existence of the discrete solution.

Based on a priori estimates and compactness arguments, we prove the convergence of the numerical approximation to the weak solution of the continuous problem. Furthermore, we have integrated the numerical scheme in the DuMu X framework. Numerical results proving the robustness and the ability of our module to deal with heterogeneity of the aquifer are presented.

In chapter 6, we conclude and summarize the work carried out in this thesis and give some research perspectives in the context of modeling seawater intrusion in coastal aquifers.

-α = ρ f ρs-ρ f : Density contrast [-], -γ = ρs-ρ f ρ f , -Z B : Bottom of the aquifer [m], -Z T : Water table elevation [m], -δ: Width of the transition zones [m]. Functions -b f = Z T -Z B : Thickness of freshwater [m], -b s = Z -Z B : Thickness of saltwater [m], -ϕ(u) = 1 2 (u -Z B ) 2 , -ϕ(Z) = 1 2 (Z -Z B ) 2 .

Introduction

In this chapter, we are interested in the mathematical modeling of seawater intrusion in coastal aquifers. In section 2.2, we present an overview of marine intrusion models. We consider two classical approaches. The sharp interface approach (2D model), see for instance [START_REF] Bear | Modelling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Bear | Seawater Intrusion in Coastal Aquifers, Concepts, Methods and Practices[END_REF][START_REF] Bear | Modelling Groundwater Flow and Pollution[END_REF][START_REF] Cheng | Multilayered Aquifier Systems[END_REF][START_REF] Cheng | Coastal Aquifer Management-Monitoring, Modeling, and Case Studies[END_REF][START_REF] Diersch | FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media[END_REF], which assumes that freshwater and saltwater are immiscible and separated by a sharp interface. Furthermore, we have the variable density approach (3D model) which supposes that the two fluids are miscible and a mixing zone caused by dispersion is considered. In this context, another approach modeling saltwater intrusion in coastal aquifers (2D model) has recently been introduced in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. It is based on the assumption that the two fluids are immiscible and takes into account the dynamics of transition zones.

Chapter 2. Mathematical modeling of seawater intrusion in coastal aquifers

In section 2.3, we present the mathematical formulation of the equations modeling marine intrusion using the sharp interface model. The hydraulic flow of freshwater and saltwater is governed by Darcy's law and the mass conservation equation of each fluid. We define the main unknowns and physical parameters of the aquifer and subsequently introduce the laws that contribute in the construction of the mathematical system. Furthermore, we consider two hypotheses, corresponding to Dupuit's approximation [START_REF] Dupuit | Traité théorique et pratique de la conduite et de la distribution des eaux[END_REF], which suppose that the flow is horizontal and Bear's hypothesis [START_REF] Bear | Dynamics of fluids in porous media[END_REF], where the vertical variation of the fluids is neglected. Thereafter, we present the sharp interface seawater intrusion model in a coastal aquifer. Such flows are governed by a coupled non-linear system of two parabolic equations simulating freshwater and saltwater heads respectively. Finally, section 2.4 contains concluding remarks.

A state of the art on marine intrusion models

Groundwater is the main source of drinking water in coastal areas, and better management of these resources is becoming a major challenge, see e.g. [START_REF] Werner | Seawater intrusion processes, investigation and management: recent advances and future challenges[END_REF]. With the population explosion and increased economic activities in these areas, the high demand for drinking water is threatening the water potential of coastal aquifers, both in terms of quantity and quality. Indeed, over-exploitation can lead to the invasion of aquifers by seawater, thus degrading groundwater quality. Marine intrusion is caused by excessive pumping from the water table, thus reducing the load of freshwater and destabilizing the prevailing balance between freshwater and saltwater. Therefore, seawater progresses towards the land side as a saltwater bevel under freshwater because of its high density compared to freshwater. The implementation of appropriate strategies for the management of coastal aquifers is a requirement. Much work is also being done to develop technical tools to determine the amount of water that can be pumped without causing such a negative impact.

There are many reference books, see for instance [START_REF] Bear | Modelling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Bear | Seawater Intrusion in Coastal Aquifers, Concepts, Methods and Practices[END_REF][START_REF] Bear | Modelling Groundwater Flow and Pollution[END_REF][START_REF] Cheng | Multilayered Aquifier Systems[END_REF], and different theses, such as [START_REF] Abudawia | Analyse numérique d'une approximation élément fini pour un modèle d'intrusion saline dans les aquifères côtiers[END_REF][START_REF] Diédhiou | Approche mixte interface nette/diffuse pour les problèmes d'intrusion saline en sous-sol: Modélisation, analyse mathématique et illustrations numériques[END_REF][START_REF] Kalaoun | Modélisation de l'intrusion saline dans l'aquifère de Tripoli-Liban: impact des changements globaux[END_REF][START_REF] Mourad | Identification de la conductivité hydraulique pour un problème d'intrusion saline: Comparaison entre l'approche déterministe et l'approche stochastique[END_REF][START_REF] Oulhaj | Conception et analyse de schémas non-linéaires pour la résolution de problèmes paraboliques: application aux écoulements en milieux poreux[END_REF][START_REF] Zidane | Risk of subsidence and aquifer contamination due to evaporite dissolution: Modelization of flow and mass transport in porous and free flow domains[END_REF] on the modeling of marine intrusion in coastal aquifers. Two approaches has been adopted. The first one is the variable density flow and solute transport approach, which assumes that freshwater and saltwater are miscible and a mixing zone caused by dispersion is tacking into account. We can mention in this case [START_REF] Cobaner | Three dimensional simulation of seawater intrusion in coastal aquifers: A case study in the goksu deltaic plain[END_REF][START_REF] Fahs | A generalized semi-analytical solution for the dispersive Henry problem: effect of stratification and anisotropy on seawater intrusion[END_REF][START_REF] Houben | The influence of heterogeneity on coastal groundwater flow-physical and numerical modeling of fringing reefs, dykes and structured conductivity fields[END_REF][START_REF] Povich | Finite element methods for variable density flow and solute transport[END_REF]. The second method is based on an abrupt approximation, see e.g. [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Huyakorn | Multiphase approach to the numerical solution of a sharp interface saltwater intrusion problem[END_REF]. This approach assumes that the two fluids are immiscible, therefore separated by a sharp interface.

Models with a sharp interface can be divided into two main categories. Models considering only freshwater flow and saltwater are assumed to be stationary. In this case, we are talking about the one fluid dynamics approach of freshwater [START_REF] Aharmouch | Development and evaluation of a numerical model for steady state interface and/or free surface groundwater flow[END_REF][START_REF] Essaid | A comparison of the coupled fresh water-salt water flow and the ghybenherzberg sharp interface approaches to modeling of transient behavior in coastal aquifer systems[END_REF][START_REF] Essaid | A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: model development and application[END_REF]. This approach remains valid for long-term predictions under steady-state conditions but suffers from an insufficient description of the behavior of the freshwater-saltwater interface in the short term. Secondly, we have the two-fluid dynamics approach [START_REF] Bouzouf | A comparison of finite volume method and sharp model for two dimensional saltwater intrusion modeling[END_REF][START_REF] Huyakorn | Multiphase approach to the numerical solution of a sharp interface saltwater intrusion problem[END_REF] considering the dynamics of the two fluids (freshwater and saltwater). Different analytical models are derived from the sharp interface approach for the prediction of marine intrusion [START_REF] Ferris | Cyclic fluctuations of water level as a basis for determining aquifer transmissibility[END_REF][START_REF] Glover | The pattern of fresh-water flow in a coastal aquifer[END_REF][START_REF] Keulegan | An example report on model laws for density current[END_REF][START_REF] Strack | A single-potential solution for regional interface problems in coastal aquifers[END_REF][START_REF] Vappicha | An approximate solution for the transient interface in a coastal aquifer[END_REF] but they are limited to very simple geometries and to situations where the Ghyben-Herzberg hypothesis is satisfied. Under simplified conditions, these analytical solutions are used as a test case to validate numerical codes.

Recently, several studies [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF][START_REF] Abudawia | Numerical analysis for a seawater intrusion problem in a confined aquifer[END_REF][START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF][START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF][START_REF] Choquet | Mathematical analysis of a seawater intrusion model including storativity[END_REF] have focused on saline intrusion modeling with a new method called sharp-diffuse interfaces approach. The last one combines the simplicity of the sharp interface model with the physical accuracy of models with a diffuse interface. The modeling of the transition zone is based on a phase-field theory developed in [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening[END_REF] and used for the description of the exchange between freshwater and saltwater. A detailed description of the derivation of the mathematical model is given in [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF].

Sharp interface approach

The sharp interface approach assumes that saltwater and freshwater are immiscible separated by an abrupt interface. This approximation is all the more legitimate as the vertical dimension of the aquifer is very small in front of the horizontal dimensions of the aquifer, see for instance [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Essaid | A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: model development and application[END_REF]. For the same reasons, the interface between the saturated and unsaturated zone is assumed to be abrupt in the case of a free aquifer. The presentation of the sharp interface approach is shown in Figure 2.1. In steady state conditions, according to a hypothesis proposed in [START_REF] Dagan | Solving the problem of local interface upcoming in a coastal aquifer by the method of small perturbations[END_REF], there is an approximate relationship between the thickness of freshwater above sea level and that of freshwater below sea level. In hydrostatic conditions, they assumed that saltwater is static and freshwater is the only fluid in dynamic. According to this relation, if the water table in an unconfined coastal aquifer is lowered by 1 m, the saltwater interface will rise 40 m. In [START_REF] Werner | Correction factor to account for dispersion in sharp-interface models of terrestrial freshwater lenses and active seawater intrusion[END_REF], a density difference factor was modified in the relationship introduced in [START_REF] Dagan | Solving the problem of local interface upcoming in a coastal aquifer by the method of small perturbations[END_REF]. The developed method is used as an analytical solution to predict seawater intrusion in steady-state conditions in [92].

Unlike previous two-fluid flow models that were established using freshwater and saltwater hydraulic heads as variables, the model introduced in [START_REF] Huyakorn | Multiphase approach to the numerical solution of a sharp interface saltwater intrusion problem[END_REF] used a mixed formulation with one fluid potential and a Pseudosaturation as double variables. The numerical solution solved by the finite element method was verified against the analytical solution proposed in [START_REF] Strack | A single-potential solution for regional interface problems in coastal aquifers[END_REF]. In the same context, another mathematical formulation of the sharp interface model is proposed in [START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF]. It is based on the mathematical model presented in [START_REF] Essaid | A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: model development and application[END_REF] which simulates freshwater and saltwater heads, therefore the position of the interface is calculated using the continuity of the flux on the interface. The variables of the proposed model are the depth of the salt front and the capillary fringe.

In the case of confined aquifers, the validity of the sharp interface model was assessed in [START_REF] Sakr | Validity of a sharp-interface model in a confined coastal aquifer[END_REF] according to four factors: infiltration factor, dispersion/advection ratio, geometric ratio, and time scale factor. The results have shown that the abrupt interface is applicable in steady and unsteady state conditions provided that advection is dominant. Similarly, its validity was discussed in [START_REF] Llopis-Albert | Discussion about the validity of sharpinterface models to deal with seawater intrusion in coastal aquifers[END_REF] by comparison with the diffusive approach in a free aquifer subjected to pumping. Considering the mathematical formulation presented in [START_REF] Essaid | A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: model development and application[END_REF], a numerical model based on a finite volume method is performed to predict seawater intrusion in Llobregat Delta [START_REF] Bouzouf | A comparison of finite volume method and sharp model for two dimensional saltwater intrusion modeling[END_REF][START_REF] Bouzouf | Integrating hydrogeochemical and geophysical data for testing a finite volume based numerical model for saltwater intrusion[END_REF].

Recent work has examined other factors that increased seawater intrusion in a coastal aquifer. The impact of Sea-Level Rise (SLR) on seawater intrusion has been investigated in [START_REF] Ketabchi | Sea-level rise impacts on seawater intrusion in coastal aquifers: Review and integration[END_REF][START_REF] Vafaie | Investigation of sea level rise effect on saltwater intrusion in an unconfined coastal aquifer using sharp-interface approach[END_REF] for a single aquifer and in [START_REF] Mehdizadeh | Sea level rise effect on seawater intrusion into layered coastal aquifers (simulation using dispersive and sharpinterface approaches)[END_REF] for multi-layer aquifer. The location of the toe (intersection of the interface with the base of the aquifer) is assessed in [START_REF] Ketabchi | Sea-level rise impacts on seawater intrusion in coastal aquifers: Review and integration[END_REF] under the impact of SLR combined with the Land Surface inundation and the Aquifer Bed Slope. The pumping impact was studied in [START_REF] Goswami | Laboratory-scale investigation of saltwater intrusion dynamics[END_REF][START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF][START_REF] Shi | Applicability of a sharp-interface model for estimating steady-state salinity at pumping wells-validation against sand tank experiments[END_REF] according to two parameters: the pumping quantity and the position of the wells. Several scenarios are presented to visualize the motion of the salt front while the effect of the random heterogeneity of the aquifer structure on saline intrusion was investigated in [START_REF] Dagan | Seawater-freshwater interface in a stratified aquifer of random permeability distribution[END_REF]. The resolution of the system in a multilayer aquifer is considered as a kind of heterogeneity of the medium [START_REF] Huyakorn | Multiphase approach to the numerical solution of a sharp interface saltwater intrusion problem[END_REF][START_REF] Mehdizadeh | Sea level rise effect on seawater intrusion into layered coastal aquifers (simulation using dispersive and sharpinterface approaches)[END_REF][START_REF] Mehdizadeh | Vertical leakage in sharp-interface seawater intrusion models of layered coastal aquifers[END_REF]. The impact of the population increase on the exploitation of of freshwater potential in the Tripoli aquifer was examined in [START_REF] Kalaoun | Assessing the contribution of demographic growth, climate change, and the refugee crisis on seawater intrusion in the Tripoli aquifer[END_REF]. Results for different pumping scenarios and their impact on the evolution of the freshwater/saltwater interface are presented.

Many works are presented for the application of the sharp groundwater model in complex and realistic examples. We can cite the case of Tongatapu Island [START_REF] Babu | Sharp interface approach for regional and well scale modeling of small island freshwater lens: Tongatapu island[END_REF], the Tripoli Aquifer in Lebanon [START_REF] Kalaoun | Assessing the contribution of demographic growth, climate change, and the refugee crisis on seawater intrusion in the Tripoli aquifer[END_REF], the Qom aquifer in Iran [START_REF] Ranjbar | Development of an efficient surrogate model based on aquifer dimensions to prevent seawater intrusion in anisotropic coastal aquifers, case study: the Qom aquifer in Iran[END_REF], the Greek Island Kalymnos [START_REF] Kopsiaftis | Pumping optimization in coastal aquifers: comparison of sharp interface and density dependent models[END_REF], the Llobregat aquifer in Spain [START_REF] Bouzouf | A comparison of finite volume method and sharp model for two dimensional saltwater intrusion modeling[END_REF], the Soquel-Aptos basin -Santa Cruz Country in California [START_REF] Essaid | A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: model development and application[END_REF], the Walawe River basin in Sri Lanka [START_REF] Ranjan | Numerical modelling of saltwaterfreshwater interaction in the Walawe river basin, Sri Lanka[END_REF], Yermasogia aquifer in Cyprus and the Crau aquifer in France [START_REF] Gemitzi | Hydra model: Simulation of salt intrusion in coastal aquifers using visual basic and gis[END_REF].

Variable density approach

In the variable density approach, the two volumes of water are considered to interact with each other, thus forming a transition zone characterized by variations in salt concentration. This interface is, in fact, a diffuse zone in which freshwater and saltwater mix. It is called a dispersion zone (or transition zone) (Figure 2.2). This approach seems more realistic and represents physical reality since fresh and saltwater are miscible. The vertical section is then considered and the salt concentration is expressed in 3D. The density is considered as a function of salt concentration. During the last decades, many simulations models were developed with finite elements and finite difference techniques using this approach [START_REF] Chun | Assessing impacts of climate change and sea-level rise on seawater intrusion in a coastal aquifer[END_REF][START_REF] Cobaner | Three dimensional simulation of seawater intrusion in coastal aquifers: A case study in the goksu deltaic plain[END_REF][START_REF] Hamidi | Modeling of 2D density-dependent flow and transport in porous media using finite volume method[END_REF][START_REF] Karahanoglu | Finite element simulation of seawater intrusion into a quarry-site coastal aquifer, Kocaeli-Darıca, Turkey[END_REF][START_REF] Ketabchi | Sea-level rise impacts on seawater intrusion in coastal aquifers: Review and integration[END_REF], also the finite volume methods [START_REF] Badaruddin | Water table salinization due to seawater intrusion[END_REF][START_REF] Hamidi | Numerical modeling of seawater intrusion in coastal aquifer using finite volume unstructured mesh method[END_REF][START_REF] Liu | A finite volume simulation model for saturated-unsaturated flow and application to Gooburrum, Bundaberg, Queensland, Australia[END_REF]. In the same context, different works are devoted to the development of software for the problem with variable density such as SEAWAT [START_REF] Guo | User's guide to SEAWAT; a computer program for simulation of three-dimensional variable-density ground-water flow[END_REF], SUTRA [START_REF] Voss | A finite-element simulation model for saturated-unsaturated, fluid-densitydependent ground-water flow with energy transport or chemically-reactive singlespecies solute transport[END_REF] based on finite element methods, FEMWATER [START_REF] Lin | A three-dimensional finite-element computer model for simulating density-dependent flow and transport in variable saturated media: version 3[END_REF], MODFLOW [START_REF] Mcdonald | A modular three-dimensional finite-difference ground-water flow model[END_REF], MOCDENS-3D [START_REF] Sanford | A two-constituent solute-transport model for ground water having variable density[END_REF] and FEFLOW [START_REF] Diersch | Variable-density flow and transport in porous media: approaches and challenges[END_REF].

Under steady-state conditions, a semi-analytical solution proposed in [START_REF] Henry | Effects of dispersion on salt encroachment in coastal aquifers, in " seawater in coastal aquifers[END_REF] for the variable density diffusion problem in a confined aquifer. The classical problem of saltwater intrusion in a coastal aquifer using the variable density approach under simplified conditions, introduced in [START_REF] Henry | Effects of dispersion on salt encroachment in coastal aquifers, in " seawater in coastal aquifers[END_REF], is widely used for the validation of numerical models. All semi-analytical methods proposed in previous works [START_REF] Fahs | A new benchmark semi-analytical solution for density-driven flow in porous media[END_REF][START_REF] Koohbor | Semianalytical solutions for contaminant transport under variable velocity field in a coastal aquifer[END_REF][START_REF] Younes | A semi-analytical solution for the reactive henry saltwater intrusion problem[END_REF][START_REF] Younes | Extension of the Henry semi-analytical solution for saltwater intrusion in stratified domains[END_REF][START_REF] Zidane | Risk of subsidence and aquifer contamination due to evaporite dissolution: Modelization of flow and mass transport in porous and free flow domains[END_REF] for the resolution of Henry's problem are limited to homogeneous and isotropic cases that have recently extended to a heterogeneous and anisotropic coastal aquifer [START_REF] Fahs | A generalized semi-analytical solution for the dispersive Henry problem: effect of stratification and anisotropy on seawater intrusion[END_REF]. This model has many applications in randomly heterogeneous aquifers [START_REF] Siena | Groundwater withdrawal in randomly heterogeneous coastal aquifers[END_REF] and fractured media [START_REF] Sebben | Seawater intrusion in fractured coastal aquifers: A preliminary numerical investigation using a fractured henry problem[END_REF]. For the application of the method in realistic case, we can mention, the Kocaeli-Darıca plain in Turkey [START_REF] Karahanoglu | Finite element simulation of seawater intrusion into a quarry-site coastal aquifer, Kocaeli-Darıca, Turkey[END_REF], the Wadi Ham aquifer in the Fujairah Emirate of the UAE [START_REF] Hussain | Three dimensional simulation of seawater intrusion in a regional coastal aquifer in UAE[END_REF], the Gooburrum-Bundaberg-Queensland aquifers in Australia [START_REF] Liu | A finite volume simulation model for saturated-unsaturated flow and application to Gooburrum, Bundaberg, Queensland, Australia[END_REF].

Sharp-diffuse interfaces approach

A new model combined between the sharp interface approach and the diffusive approach has recently been introduced in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. The freshwater and saltwater zones are modeled using the Dupuit approach which assumes that all fluid movements are essentially horizontal.

The modeling of the diffusive interface (resp. between the saturated and unsaturated zone) is derived from a model based on the phase field theory. However, the latter was initially developed in [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening[END_REF] for the description of the transition phenomena between freshwater and seawater. The reader can consult [START_REF] Abudawia | Analyse numérique d'une approximation élément fini pour un modèle d'intrusion saline dans les aquifères côtiers[END_REF][START_REF] Diédhiou | Approche mixte interface nette/diffuse pour les problèmes d'intrusion saline en sous-sol: Modélisation, analyse mathématique et illustrations numériques[END_REF][START_REF] Mourad | Identification de la conductivité hydraulique pour un problème d'intrusion saline: Comparaison entre l'approche déterministe et l'approche stochastique[END_REF] for documentation including an overview of the sharp-diffuse modeling approach and the mathematical formulation of the proposed model. It combines the simplicity of the abrupt interface with the realism of diffuse interface models. Furthermore, the model takes into account the dynamics of the diffuse interface. In [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF], a numerical comparison between the results given by this model with those obtained from the variable-density flow and solute transport model, is proposed. Several scenarios, including pumping and injection processes, were considered to demonstrate that the diffuse interface provided by this model passes through the middle of the transition zone. The salt concentration in the transition zone is assumed to be half that of the saltwater zone.

Another comparison between this model and the classical sharp interface approach is proposed in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] to assess the validity of the mixte approach. The diffusive interface given by the mixed approach coincides with that provided by the sharp interface model when the width of the transition zone is small. Tidal effects on the saltwater front were studied in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] using this model. The effect of the width of the transition zone was also investigated in the same reference. In [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF], the authors approximate the sharp-diffuse model using finite volume and finite element methods in unconfined aquifers. The numerical results were compared and showed good agreement.

Mathematical modeling of the sharp interface model

In this section, we present the basic equations used in the derivation of the mathematical model, including the mass conservation equations and Darcy's law. The governing flow equations are provided firstly in 3D. Based on different assumptions, presented in subsection 2.3.3, the equations are vertically integrated and the problem is reduced to a 2D model. Afterwards, the system of equations is closed by evaluating the flow and the pressure continuity on the sharp interface.

Continuity equation

The continuity equation of a fluid in movement on a Representative Elementary Volume (REV) is given by the following equation, see e.g. [START_REF] Bear | Modelling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Bear | Modelling Groundwater Flow and Pollution[END_REF]]

∂(φρ) ∂t + div(ρ q) = ρQ, (2.3.1)
with:

φ: porosity [-],
ρ: fluid density kg • m -3 , q: Darcy's velocity m.s -1 , Q: quantity of water injected per unit of VER s -1 .

We define the coefficient of compressibility of the soil and the fluid respectively by Therefore, the mass conservation equation (2.3.1) is written as follows

β P = 1 (1 -
ρ φα P + β P (1 -φ) ∂P ∂t + ρ div( q) + ∇ρ • q = ρQ. (2.3.3)

Darcy's law

Darcy's law describes the flow of a fluid through a porous medium in a macroscopic context. It expresses the flow velocity as a function of the gradient of the hydraulic head.

It is given by the following equation: We define the macroscopic hydraulic head by

q = - κ µ (∇P + ρg∇x 3 ) , (2.3.4)
Φ[m] = P ρg + x 3 . (2.3.5)
The freshwater and the saltwater are incompressibles with a constant density. In this case, the hydraulic head is identified at the piezometric level. Equation (2.3.4) is, therefore, written as a function of the hydraulic head:

q = -D∇Φ, (2.3.6)
where D m.s -1 = κρg µ is the tensor of the hydraulic conductivity.

Assumptions

In this subsection, we present different assumptions used in the derivation of the sharp interface model. Some hypotheses are related to the nature of the soil and the fluids and others are used for the up-scaling of the 3D problem to a 2D model. The aquifer is supposed weakly deformable and the elastic behavior of the soil is very weak. Moreover, Chapter 2. Mathematical modeling of seawater intrusion in coastal aquifers freshwater and saltwater are both weakly compressible fluids. Therefore, the compressibility coefficient of the fluid and the soil, defined in (2.3.2), are neglected, i.e. (β p , α P 1).

Two main assumptions are also applied: Dupuit's approximation where the movement of the fluids is horizontal since the thickness of many coastal aquifers is neglected compared to its horizontal surface. This has been used and discussed in [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Dupuit | Traité théorique et pratique de la conduite et de la distribution des eaux[END_REF]. In this concept, the transition zone between freshwater and saltwater (resp. saturated and unsaturated zone in the case of a free aquifer) is presumed an abrupt interface. The problem is thereby reduced to a 2D model describing the freshwater and saltwater flows separated by an abrupt interface. In a second step, we consider Bear's hypothesis [START_REF] Bear | Dynamics of fluids in porous media[END_REF], which consists to ignore the density variation in the flow direction.

Integration of the equations

This subsection aims to present different details of the derivation of the mathematical model. The mass conservation equation (2.3.3) and Darcy's equation (2.3.4) for each fluid can be rewritten as follows

       ρ ρg φα P + β P (1 -φ) S ∂Φ ∂t + ρ div( q) + ∇ρ • q = ρQ, q = -D∇Φ, (2.3.7)
with S = ρg φα P + β P (1 -φ) is the storage coefficient representing the capacity of the aquifer to contain groundwater. Usually, its value is very small due to the low compressibilty of the soil and the fluid. i.e, ( β P 1 and α P 1).

Using Bear's hypothesis, the term ∇ρ • q in equation (2.3.7) is negligible. One has ρS ∂Φ ∂t + ρ div( q) + ∇ρ • q =0 = ρQ, with q = -D∇Φ.

(2.3.8)

We introduce specific indices for freshwater (f) and saltwater (s) in (2.3.8). The governing freshwater and saltwater flows are as follows

S f ∂Φ f ∂t + div( q f ) = Q f , S f = ρ f g φα P + β P (1 -φ) , q f = -D f ∇Φ f . (2.3.9) S s ∂Φ s ∂t + div( q s ) = Q s , S s = ρ s g φα P + β P (1 -φ) , q s = -D s ∇Φ s . (2.3.10)
We consider a coastal aquifer represented by a domain Ω × (Z B , Z T ) with Ω ⊂ R 2 is a part of its horizontal surface (see ) and (2.3.10) are expressed in 3D. However, the vertical dimension of a coastal aquifer is very small compared to its horizontal dimension. We, therefore, use Dupuit's approximation. In the following, we will present all the steps of upscaling the 3D problem to a 2D model. In this part, we refer to the work performed in [START_REF] Diédhiou | Approche mixte interface nette/diffuse pour les problèmes d'intrusion saline en sous-sol: Modélisation, analyse mathématique et illustrations numériques[END_REF].

We integrate equation (2.3.9) between Z and Z T . We obtain

Z T Z S f ∂Φ f ∂t + div( q f ) dx 3 = Z T Z Q f dx 3 . (2.3.11)
According to the Leibniz formula, we have

Z T Z S f ∂Φ f ∂t dx 3 = S f Z T Z ∂Φ f ∂t dx 3 , = S f ∂ ∂t Z T Z Φ f dx 3 -S f Φ f (Z T ) ∂Z T ∂t + S f Φ f (Z) ∂Z ∂t .
(2.3.12)

Let Φ f be the vertical average of Φ f :

Φ f (x 1 , x 2 ) = 1 b f Z T Z Φ f (x 1 , x 2 , x 3 ) dx 3 , ∀(x 1 , x 2 ) ∈ Ω. Since S f ∂ ∂t Z T Z Φ f dx 3 = S f ∂ ∂t b f Φ f , = S f b f ∂Φ f ∂t + S f Φ f ∂b f ∂t , equation (2.3.12) becomes Z T Z S f ∂Φ f ∂t dx 3 = S f b f ∂Φ f ∂t + S f Φ f -Φ f (Z T ) ∂Z T ∂t -S f Φ f -Φ f (Z) ∂Z ∂t .

Chapter 2. Mathematical modeling of seawater intrusion in coastal aquifers

Similarly, applying the Leibniz formula to the second term on the right-hand side of equation (2.3.11), we obtain

Z T Z div ( q f ) = Z T Z ∂q fx 1 ∂x 1 + ∂q fx 2 ∂x 2 dx 3 + Z T Z ∂q fx 3 ∂x 3 dx 3 , = ∇ • Z T Z q f dx 3 -q f (Z T ) • ∇ (Z T ) + q f (Z) • ∇ (Z) + q fx 3 (Z T ) -q fx 3 (Z), = ∇ • b f q f -q f (Z) • ∇ (Z T ) -∂ x 3 x 3 + q f (Z) • ∇ (Z) -∂ x 3 x 3 ,
where:

q f = q fx 1 , q fx 2 , q fx 3 q f = (q fx 1 , q fx 2 ) ∇ • () = ∂() ∂x 1 + ∂() ∂x 2 + ∂() ∂x 3 , ∇ • () = ∂() ∂x 1 + ∂() ∂x 2 .
The vertical average of freshwater flux is given by

q f 1 b f Z T Z q f dx 3 -D f ∇ Φ f , with D f (x 1 , x 2 ) = 1 b f Z T Z D f (x 1 , x 2 , x 3 ) dx 3 , ∀(x 1 , x 2 )
∈ Ω, and x 3 ∈ (Z, Z T ) Using Dupuit's assumption, the following approximations hold up Φ

f (x 1 , x 2 , x 3 ) Φ f (x 1 , x 2 ) and Q f (x 1 , x 2 , x 3 ) Q f = (x 1 , x 2 ) with (x 1 , x 2 )
∈ Ω and x 3 ∈ (Z, Z T ). One obtains the equation modeling freshwater flow process [START_REF] Essaid | A comparison of the coupled fresh water-salt water flow and the ghybenherzberg sharp interface approaches to modeling of transient behavior in coastal aquifer systems[END_REF][START_REF] Essaid | A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: model development and application[END_REF]

S f b f ∂Φ f ∂t -div b f D f ∇ Φ f + q f (Z T ) • ∇(x 3 -Z T ) -q f (Z) • ∇(x 3 -Z) = Q f . (2.3.13) Likewise, using Q s (x 1 , x 2 , x 3 ) Q s = (x 1 , x 2
) and integrating between Z B and Z the equation (2.3.10), leads to the following equation describing the saltwater flow

S s b s ∂Φ s ∂t -div b s D s ∇ Φ s + q s (Z) • ∇(x 3 -Z) -q s (Z B ) • ∇(x 3 -Z B ) = Q s . (2.3.14)
The vertical averages of Φ s and D s are given by

Φ s (x 1 , x 2 ) = 1 b s Z Z B Φ s (x 1 , x 2 , x 3 ) dx 3 , and D s (x 1 , x 2 ) = 1 b s Z Z B D s (x 1 , x 2 , x 3 ) dx 3 .
The derived system (2.3.13)-(2.3.14) is composed of two equations with four unknowns ( Φ f , Φ s , Z and Z T ). To reduce the number of unknowns, we close the system by approaching the fluxes at the interfaces. More precisely, we use the continuity of the fluxes and the pressure at the interfaces.

Pressure continuity at the interface

We assume the continuity of the pressure at the interface. The salt front elevation Z can be written as a function of the hydraulic head (2.3.5)

P f |x 3 =Z = P s|x 3 =Z ⇐⇒ ρ f (Φ f -Z) = ρ s (Φ s -Z), ⇐⇒ Z = 1 + ρ f ρs-ρ f Φ s - ρ f ρs-ρ f Φ f , ⇐⇒ Z = (1 + α)Φ s -αΦ f , (2.3.15) with α = ρ f ρs-ρ f is the density contrast. Approximation of the flux q s|x 3 =Z B • ∇(x 3 -Z B )
The base layer of the aquifer is impermeable. There is, therefore, no flow between the salt zone and the bottom x 3 = Z B . One has

q s|z=Z B • ∇(x 3 -Z B ) = 0.
The base of the aquifer Z B is assumed spatially variable.

Approximation of the flux

q i|x 3 =Z • ∇(x 3 -Z), i = f, s
In the case of the sharp interface approach, the normal components of the velocity at the interface are equal in each zone. This corresponds to no transfer of mass across the interface. Thus, we obtain

q f |x 3 =Z φ -v • n = q s|x 3 =Z φ -v • n = 0,
where n is the normal unit vector with respect to the abrupt interface, v m.s -1 is the interface velocity. We get

q f |x 3 =Z • ∇(x 3 -Z) = q s|x 3 =Z • ∇(x 3 -Z) = φ ∂Z ∂t .
Approximation of the flux

q f |x 3 =Z T • ∇(x 3 -Z T )
The modeling of the upper surface depends on the nature of the aquifer, we will distinguish between free and confined aquifers.

In the case of a confined aquifer, the upper surface of the aquifer is impermeable, so there is no flow between the roof and the freshwater area. The flux through the upper interface of the aquifer is given by:

q f |x 3 =Z T • ∇(x 3 -Z T ) = 0.
In the case of a free aquifer, the upper surface consists of a permeable layer, which can be crossed by water at reasonable speeds. This limit corresponds to the piezometric surface where there are exchanges with the outside world such as injections and rainwater supply. The flow is given by

q f |x 3 =Z T • ∇(x 3 -Z T ) = φ ∂Φ f ∂t .
We introduce a parameter β to combine the two cases of the confined and unconfined aquifer,

β = 1 if the aquifer is free, 0 if the aquifer is confined.
By replacing the flux in the equations (2.3.13) and (2.3.14) and Z by its formula, we obtain a 2D system modeling freshwater and saltwater flows

S f b f + (β + α)φ(x) ∂Φ f ∂t -(1 + α)φ(x) ∂Φ s ∂t -div b f D f ∇Φ f = Q f , S s b s + (1 + α)φ(x) ∂Φ s ∂t -αφ(x) ∂u ∂t -div b s D s ∇Φ s = Q s .
The thicknesses of freshwater and saltwater are respectively

b f = Z T -Z and b s = Z -Z B .
The salt front is obtained by using the continuity of pressure on the interface Z = (1 + α)Φ s -αΦ f . The upper surface Z T is replaced by Φ f in the case of a free aquifer and by a constant in the confined case.

Mathematical model of the sharp interface problem

In this subsection, we present the seawater intrusion sharp interface model, where freshwater and saltwater heads are the unknowns. Let Ω be an open polygonal domain,

∂Ω = Γ = Γ D ∪ Γ N , x = (x 1 , x 2 ) ∈ R 2 and ]0, T [ a time interval. Let Ω T =]0, T [×Ω. We note u = Φ f , v = Φ s and D i = D i , i = f, s.
The mathematical system is given by:

                   β f (u, v) ∂u ∂t -ω(.) ∂v ∂t -div (b f (u, v)D f ∇u) = Q f in Ω T , β s (u, v) ∂v ∂t -λ(.) ∂u ∂t -div (b s (u, v)D s ∇v) = Q s in Ω T , u(0, .) = u 0 (.), v(0, .) = v 0 (.) in Ω, v = v D , u = u D on ]0, T [ × Γ D , b f (u, v)D f ∇u • n = 0, b s (u, v)D s ∇v • n = 0 on ]0, T [×Γ N . (2.3.16) with β f (u, v) = S f b f (u, v) + φ(x)(β + α), ω(x) = φ(x)(1 + α), β s (u, v) = S s b s (u, v) + φ(x)(1 + α), and λ(x) = φ(x)α.
The thickness of freshwater is given by

b f (u, v) = (1 + α)(u -v) free aquifer, Z T -(1 + α)v + αu confined aquifer, while that of saltwater is as follows b s (u, v) = αu -(1 + α)v -Z B .
The initial conditions u 0 (x) and v 0 (x) are given for u and v at (t=0) to describe the initial state of the aquifer.

The problem is a coupled system of two nonlinear parabolic partial differential equations simulating two immiscible fluids. The unknowns are u the hydraulic head of freshwater and v the hydraulic head of saltwater. The depth of the freshwater/saltwater surface can be calculated from the following relationship

Z = (1 + α)v -αu.

Conclusion

In this section, we conclude the work carried out in the context of marine intrusion modeling. We first presented a state of the art of saline intrusion models. Three approaches were outlined. The classical 2D sharp interface approach, which assumes that the two fluids are immiscible. This approach simulates respectively the location of the salt front and the depth of the free surface. The second is the 3D variable density approach, which assumes that the two masses of water interact with each other. It is characterized by the variation of the salt concentration and the hydraulic head of the freshwater. Furthermore, we presented a new approach [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] corresponding to the sharp-diffuse interfaces model. It is a combination of the two previous approaches. The two fluids are assumed to be immiscible and the dynamics of the transition zones are taken into account.

The other part of this chapter is devoted to the formulation of the mathematical model using the sharp interface approach. We first introduced the physical laws of the continuity equation and Darcy's law. Moreover, we consider Bear's hypothesis and Dupuit's approximation to upscale the initial 3D problem. We presented the mathematical model of the marine intrusion using the sharp interface approach in a free and confined aquifer. The system of equations simulating freshwater and saltwater heads is further completed with Neumann and Dirichlet boundary conditions. The elevation of the salt front is predicted as a function of the head of freshwater and saltwater. The numerical resolution of this model will be the purpose of the next chapter.

Chapter 3

A finite volume method for numerical simulation of the sharp interface model

Introduction

This chapter deals with a numerical method for performing numerical simulations of the sharp interface model in a coastal aquifer. In this work, we develop a fully implicit, cellcentered finite volume method to predict seawater intrusion by means of this approach. A complete presentation of the derivation of the numerical scheme is presented in section 3.2. To approximate the diffusive terms, a Two-Point Flux Approximation (TPFA) method [START_REF] Eymard | Finite volume methods[END_REF] on structured mesh is performed. An upwind scheme is used to evaluate the convective fluxes. The nonlinear system is solved by a Newton's method, wherein each iteration a linear system is solved by an iterative Krylov method.

After that, we have implemented and integrated this scheme in a new module in the DuMu X framework [START_REF]DuMu X DUNE for multi-{ Phase, Component, Scale, Physics[END_REF]. A short introduction of the platform is presented in section 3.

model

The developed module enables us to solve a coupled system of two parabolic equations modeling saline intrusion in coastal aquifers. Our numerical scheme has been tested and validated by different 1D, 2D test cases, including real-scale Benchmarks. First of all, we validate our model in the case of a confined aquifer with the rotating interface problem of Keulegan [START_REF] Keulegan | An example report on model laws for density current[END_REF]. Numerical results, presented in subsection 3.4.1, are compared to the analytical solution and showed a good agreement with those presented in [START_REF] Abudawia | Numerical analysis for a seawater intrusion problem in a confined aquifer[END_REF][START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF]. After that, we validate our model in the case of an unconfined aquifer. We proceed to the study of the numerical convergence of the scheme on different grid resolutions. Then, we applied our approach to two real test cases: the first test "field-scale free aquifer" is presented in [START_REF] Mehdizadeh | Vertical leakage in sharp-interface seawater intrusion models of layered coastal aquifers[END_REF]. Numerical simulations over 30 years of exploitation of the aquifer are presented in subsection 3.4.2. The comparison between our results and those presented in [START_REF] Mehdizadeh | Vertical leakage in sharp-interface seawater intrusion models of layered coastal aquifers[END_REF] for eleven scenarios illustrates the robustness and efficiency of our approach. The second test case, presented in [5], corresponds to the Souss-Chtouka plain located in the southwest of Morocco. Hydrological data and physical parameters are used to predict seawater intrusion in the Souss-Chtouka aquifer. The numerical results showing the evolution of different quantities during 80 years of over-exploitation of the aquifer are presented in subsection 3.4.3. Afterwards, we validated our implemented module on a real test case corresponding to the Tripoli aquifer in the northeast of Lebanon. Numerical simulations over 20 years of exploitation are presented in section 3.4.4. The comparison of our numerical results and those presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF] showed a good agreement, which ensures the validity of our module to provide efficient and accurate results in real scale test cases. Section 3.5 contains concluding remarks.

A fully implicit finite volume numerical scheme

This section aims to develop a finite volume numerical scheme for the sharp interface approach. We consider the system of equations, presented in chapter 2, which simulates freshwater head (u) and saltwater head (v) into confined and unconfined aquifers:

                   β f (u, v) ∂u ∂t -ω(.) ∂v ∂t -div (b f (u, v)D f ∇u) = Q f in Ω T , β s (u, v) ∂v ∂t -λ(.) ∂u ∂t -div (b s (u, v)D s ∇v) = Q s in Ω T , u(0, .) = u 0 (.), v(0, .) = v 0 (.) in Ω, v = v D , u = u D on ]0, T [ × Γ D , b f (u, v)D f ∇u • n = 0, b s (u, v)D s ∇v • n = 0 on ]0, T [×Γ N . (3.2.1)
which involve the following functions:

β f (u, v) = S f b f (u, v) + φ(x)(β + α), ω(x) = φ(x)(1 + α), β s (u, v) = S s b s (u, v) + φ(x)(1 + α), λ(x) = φ(x)α.
The thickness of freshwater is given by

b f (u, v) = (1 + α)(u -v) free aquifer, Z T -(1 + α)v + αu confined aquifer, while that of saltwater is as follows b s (u, v) = αu -(1 + α)v -Z B .
We generalize our model in both cases of free and confined aquifers by introducing a parameter β = 1, if the aquifer is free, 0, if the aquifer is confined.

The source terms Q f and Q s represent the amount of freshwater and saltwater injected ( or pumped) per unit surface, respectively. The domain Ω is illustrated in Figure 3.1

Discretization and basic notations

In this subsection, we give the parameters and the notations which used in the derivation of the numerical scheme. Here, we choose a cell-centered finite volume methods for the discretization in space. For the sake of simplicity, we present the scheme on a structured mesh which suitable with TPFA method. The main advantage of this scheme is to preserve the conservation property of the mass and the positivity of the freshwater and saltwater thicknesses. The extension to unstructured grids will be presented in the next chapter for the so called the Multi Point Flux Approximation (MPFA) [3].

We introduce the following notations

Let {t 0 , • • • , t N } be a partition of [0, T ], ∆t n = t n+1 -t n the time step size and ∆t = max n ∆t n . let (T ) be a rectangular partition of Ω, such that Ω = ∪ V k ∈T V k . ∂V k = V k \ V k the boundary of V k .
The physical parameters are given as constant values on each cell. One has

φ k = 1 | V k | V k φ(x) dx, w k = 1 | V k | V k w(x) dx and λ k = 1 | V k | V k λ(x) dx
The freshwater thickness can be written as b

f (u, v) = χ(u) -(1 + α)v with χ(u) =
(1 + α)u, in the case of free aquifer, Z T + αu, in the case of confined aquifer.

We define the nonlinear function η ∈ C 1 (Ω) by:

η(u) = u 0 χ(s) ds.
Also, we define the function

ϕ(v) = (1 + α) 2 v 2 . model
For an initial condition u 0 ∈ L ∞ (Ω) (resp. v 0 ), we put :

u 0 k = 1 |V k | V k u 0 (x)dx (resp. v 0 k ).
Other notations are presented in Figure 3.2. We present the derivation of the numerical scheme using a cell-centered finite volume method in space. For the time discretization, we use an Euler implicit method. The approximation of the convective subfluxes is done with an upwind scheme. In the diffusion term, we approximate the gradient on the interface by a TPFA method. The control volume V k where the numerical scheme will be written is illustrated in Figure 3.2.

We integrate the freshwater and the saltwater equations in the system (3.2.1) on ]t n , t n+1 [× V k , we obtain:

The freshwater equation:

t n+1 tn V k β f (u, v) ∂u ∂t dx dt - t n+1 tn V k ω(x) ∂v ∂t dx dt - t n+1 tn V k div (b f (x, u, v)D f ∇u) dx dt = t n+1 tn V k Q f dx dt.
The saltwater equation:

t n+1 tn V k β s (u, v) ∂v ∂t dx dt - t n+1 tn V k λ(x) ∂u ∂t dx dt - t n+1 tn V k div (b s (x, u, v)D s ∇v) dx dt = t n+1 tn V k Q s dx dt. Let u n k = 1 |V k | V k u(x, t n ) dx (resp, v n k = 1 |V k | V k v(x, t n ) dx
) be the approximate value of u (resp. v) in the control volume V k at time t n . Using an implicit approximation in time, the numerical scheme can be written as follows

|V k | ∆t n β f (u n+1 k , v n+1 k ) u n+1 k -u n k -w k v n+1 k -v n k 3.2. A fully implicit finite volume numerical scheme 27 - l∈∂V (k) γ kl D f ∇η u n+1 -(1 + α)D f v n+1 ∇u n+1 kl • n kl dγ = Q n+1 f,k , |V k | ∆t n β s (u n+1 k , v n+1 k ) v n+1 k -v n k -λ k u n+1 k -u n k - l∈∂V (k) γ kl D s ∇ϕ v n+1 -D s αu n+1 + Z B ∇v n+1 kl • n kl dγ = Q n+1 s,k ,
where n kl is the unit outer normal to γ kl and V (k) is the set of adjacent elements of V k . We denote by

Q n+1 f,k (resp. Q n+1 s,k ) the average value of Q f (resp. Q s ) on V k at time t n+1 .
The exact diffusion fluxes are given by

F kl (u n+1 ) = - γ kl {D f } kl {∇η(u n+1 )} kl • n kl dγ kl , F kl (v n+1 ) = - γ kl {D s } kl {∇ϕ(v n+1 )} kl • n kl dγ kl ,
while the exact convectif fluxes are given by

G kl (u n+1 ) = (1 + α) γ kl {D f } kl {v n+1 } kl ∇u n+1 • n kl dγ kl , G kl (v n+1 ) = γ kl {D s } kl (α{u n+1 } kl + Z B )∇v n+1 • n kl dγ kl .
The next step consist to approximate the exact fluxes.

Numerical diffusion fluxes:

The gradient of the freshwater and saltwater heads on the interface γ kl are approximated by a TPFA method. For the computation of the diffusion coefficients D f (resp. D s ), we use an harmonic average between two adjacent elements. Here is the expression of the numerical fluxes defined as a function of the main unknowns:

F * kl (u n+1 k ) = -{D f } har kl η(u n+1 l ) -η(u n+1 k ) d kl |γ kl |, F * kl (v n+1 k ) = -{D s } har kl ϕ(v n+1 l ) -ϕ(v n+1 k ) d kl |γ kl |,
where {D f } har kl is the harmonic average of {D f } k and {D f } l , while d kl is the distance between the centers of the control volume V k and V l .

Numerical convective fluxes:

For the approximation of the convective terms, we use a fully implicit upwind scheme. The freshwater and saltwater heads are therefore evaluated implicitly and dependly on the sign of the freshwater velocity. The numerical convective fluxes at the interface γ kl are given by:

G * kl (u n+1 k ) = {v n+1 } up kl {D f } kl {∇u n+1 } kl • n kl |γ kl |(1 + α), G * kl (v n+1 k ) = α{u n+1 } up kl + Z B {D s } kl {∇v n+1 } kl • n kl |γ kl |.
with

{v n+1 } up kl = v n+1 k if {D f } kl {∇u n+1 } kl • n kl > 0 v n+1 l else (3.2.2)
Chapter 3. A finite volume method for numerical simulation of the sharp interface model

In an analogue way, the discrete value {u n+1 } up kl is evaluated on the interface γ kl . Finally, the fully implicit cell-centered finite volume numerical scheme issued from the sharp interface model (3.2.1) yields to the nonlinear coupled equations:

|V k | ∆t n β f (u n+1 k , v n+1 k ) u n+1 k -u n k -w k v n+1 k -v n k - l∈∂V (k) {D f } har kl |γ kl | η(u n+1 l ) -η(u n+1 k ) d kl -(1 + α){v n+1 } up kl u n+1 l -u n+1 k d kl = Q n+1 f,k , (3.2.3
)

|V k | ∆t n β s (u n+1 k , v n+1 k ) v n+1 k -v n k -λ k u n+1 k -u n k - l∈∂V (k) {D s } har kl |γ kl | ϕ(v n+1 l ) -ϕ(v n+1 k ) d kl -α{u n+1 } up kl + Z B v n+1 l -v n+1 k d kl = Q n+1 s,k . (3.2.4)
The initial conditions are given by

u 0 k = 1 |V k | V k u 0 (x) dx, and v 0 k = 1 |V k | V k v 0 (x) dx. (3.2.5)
For the discretization of the boundary conditions, we distinguish between two cases Homogeneous Neumann BCs on Γ N , one has

F * N (u n+1 k ) = F * N (v n+1 k ) = G * N (u n+1 k ) = G * N (v n+1 k ) = 0. (3.2.6)
Dirichlet BCs on Γ D , we have

F * D (u n+1 k ) = -{D f } k η(u D ) -η(u n+1 k ) d kl |γ k |, (3.2.7)
F * k (v n+1 k ) = -{D f } k η(v D ) -η(u n+1 k ) d kl |γ k |, G * D (u n+1 k ) = {v n+1 } up D {D f } k {∇u n+1 } D • n D |γ k |, (3.2.8) G * D (v n+1 k ) = {u n+1 } up D {D f } k {∇v n+1 } D • n D |γ k |, (3.2.9)
with

{v n+1 } up D = v n+1 k if {D f } k {∇u n+1 } kl • n kl > 0 v n+1 D else

Resolution of the nonlinear system

The numerical scheme (3.2.3)-(3.2.8) derived from the TPFA method is formulated by two coupled nonlinear equations. Solving the governing equations is not obvious. The difficulty arises from the nonlinearity and the coupling of the system. To linearise the coupled system at each time step ∆t n , we use Newton's method wherein each iteration, a linear system is solved using iterative linear solvers provided by the external library of DUNE [2]. The residual of the system of equations (3.2.3)-(3.2.4) which is a function of the discrete solution (u n+1 k , v n+1 k ) can be written by the two equations: ∀k ∈ T , and n ∈ {0, ..., N -1}

{R f } n+1 k u n+1 k , v n+1 k = 0, {R s } n+1 k u n+1 k , v n+1 k = 0. (3.2.10)
We note that the local residual function R f (resp. R s ) of the discrete freshwater (resp. saltwater) equation is given at each element V k and any time t n+1 . In our implemented module, we use some numerical differentiation techniques to approximate the derivatives of the residual in the calculation of the Jacobian matrix. Adding or subtracting epsilon from the solution in the numerical differential for each coefficient of the local Jacobian matrix requests more time within the numerical context (CPU time). However, in practice, we do not recalculate the whole matrix at each Newton iteration, but only a part of it according to the local residual error. Thereafter, using some numerical strategies, we assemble the global residual and then the global matrix. To show and illustrate the matrix form of the linear system solved at each Newton's iteration, we propose to compute explicitly the derivatives of the local residue from our numerical scheme (3.2.3)-(3.2.4). To do so, we consider a control volume V k sharing an edge l ∈ ∂V (k) with its neighbor V l . We begin by computing the local residual of freshwater and saltwater with respect to the discrete solution u n+1 k and v n+1 k . We introduce the following notations

α v kl = 1 if {D f } kl {∇u n+1 } kl • n kl > 0 0 else α u kl = 1 if {D f } kl {∇v n+1 } kl • n kl > 0 0 else (3.
∂{R f } n+1 k ∂u n+1 k = |V k | ∆t n ∂β f (u n+1 k , v n+1 k ) ∂u n+1 k u n+1 k -u n k + β f (u n+1 k , v n+1 k ) - l∈∂V (k) {D f } har kl ∂ ∂u n+1 k η(u n+1 l ) -η(u n+1 k ) d kl |γ kl | + (1 + α) l∈∂V (k) {D f } har kl α v kl v n+1 k ∂ ∂u n+1 k u n+1 l -u n+1 k d kl |γ kl |, (3.2
.12)

∂{R f } n+1 k ∂v n+1 k = |V k | ∆t n ∂β f (u n+1 k , v n+1 k ) ∂u n+1 k u n+1 k -u n k -ω k + (1 + α) l∈∂V (k) {D f } har kl α v kl u n+1 l -u n+1 k d kl |γ kl |. (3.2.13)
Now, we calculate the derivative of the local residual {R f } n+1 k with respect to the approximate solutions on the element V l .

∂{R f } n+1 k ∂u n+1 l = - l∈∂V (k) {D f } har Dl ∂ ∂u n+1 l η(u n+1 l ) -η(u n+1 D ) d kl |γ kl |
Chapter 3. A finite volume method for numerical simulation of the sharp interface model

+ (1 + α) l∈∂V (k) {D f } har kl α v kl v n+1 l ∂ ∂u n+1 l u n+1 l -u n+1 k d kl |γ kl |, (3.2.14) ∂{R f } n+1 k ∂v n+1 l = (1 + α) l∈∂V (k) {D f } har kl α v kl u n+1 l -u n+1 k d kl |γ kl |. (3.2.15)
The derivative of the local saltwater residual {R s } n+1 k on the element V k at time t n+1 with respect to the discrete unknowns v n+1 k and u n+1 k are as follows respectively:

∂{R s } n+1 k ∂v n+1 k = |V k | ∆t n ∂β s (u n+1 k , v n+1 k ) ∂v n+1 k v n+1 k -v n k + β s (u n+1 k , v n+1 k ) - l∈∂V (k) {D s } har kl ∂ ∂v n+1 k ϕ(v n+1 l ) -ϕ(v n+1 k ) d kl |γ kl | + l∈∂V (k) {D s } har kl α u kl u n+1 k ∂ ∂v n+1 k v n+1 l -v n+1 k d kl |γ kl |, (3.2.16
)

∂{R s } n+1 k ∂u n+1 k = |V k | ∆t n ∂β s (u n+1 k , v n+1 k ) ∂u n+1 k v n+1 k -v n k -λ k + l∈∂V (k) {D s } har kl (αα v kl + Z B ) v n+1 l -v n+1 k d kl |γ kl |. (3.2.17)
For the derivative of the local residue of the saltwater equation with respect to the discrete unknowns u n+1 l and v n+1 l , we have

∂{R s } n+1 k ∂v n+1 l = - l∈∂V (k) {D s } har kl ∂ ∂v n+1 l ϕ(v n+1 l ) -ϕ(u n+1 k ) d kl |γ kl | + l∈∂V (k) {D s } har kl α u kl u n+1 l + Z B ∂ ∂v n+1 l v n+1 l -v n+1 k d kl |γ kl |, (3.2.18) ∂{R s } n+1 k ∂u n+1 l = l∈∂V (k) {D s } har kl α u kl v n+1 l -v n+1 k d kl |γ kl |. (3.2.19)
Now we consider two adjacent elements V k and V l . We illustrate the blocks of the matrix on these two cells without considering the other cells. Therefore, the structure of the Jacobian matrix is given by

J kk J kl J lk J ll (3.2.20)
The block J kk (resp. J ll ) represents the derivative of the residual issued from the system (3.2.3)-(3.2.4) on the control volume V k (resp. V l ) with respect to the discrete unknowns of freshwater and the saltwater heads on the element V k (resp. V l ). The explicit from of the blocks is given by

J kk =      ∂{R f } n+1 k ∂u n+1 k ∂{R f } n+1 k ∂v n+1 k ∂{R s } n+1 k ∂u n+1 k ∂{R s } n+1 k ∂v n+1 k      (3.2.21) J ll =      ∂{R f } n+1 l ∂u n+1 l ∂{R f } n+1 l ∂v n+1 l ∂{R s } n+1 l ∂u n+1 l ∂{R s } n+1 l ∂v n+1 l      (3.2.22)
The blocks J kl (resp. J lk ) corresponds to the derivative of the residual on the control volume V k (resp. V l ) with respect to the discrete unknowns u n+1 l , v n+1 l on the element V l (resp. u n+1 k , v n+1 k on the element V k for the flux J lk ). The explicit form of this terms is as follows

J kl =      ∂{R f } n+1 k ∂u n+1 l ∂{R f } n+1 k ∂v n+1 l ∂{R s } n+1 k ∂u n+1 l ∂{R s } n+1 k ∂v n+1 l      (3.2.23) J lk =      ∂{R f } n+1 l ∂u n+1 k ∂{R f } n+1 l ∂v n+1 k ∂{R s } n+1 l ∂u n+1 k ∂{R s } n+1 l ∂v n+1 k      (3.2.24)

Implementation of the numerical scheme

The implementation of the code is done in the platform DuMu X . It is based on the DUNE library (Distributed and Unified Numerical Environment). It offers a variety of tools to numerically solve PDEs by the finite volume methods, finite element methods, finite difference methods and discontinuous Galerkin method. It also provides other features such as mesh management, discretization, linear and nonlinear solvers...etc.

The main feature of DuMu X is the modularity of the concepts. It can be easily integrated and combine different properties depending on the problem to be solved. For each DuMu X module, different elements interact with each other via common interfaces. The user can select different properties according to the nature of the problem.

DuMu X has two types of numerical schemes: the fully implicit approach and the sequential approach. The global implict strategy employs an implicit Euler discretization in time , combined with the "Cell-Centered " or "Vertex-Centered" method in space. The derived equations are assembled into a single nonlinear system, solved at each time step. The sequential approach deals with conservation equations where the flow equation is approximated implicitly, while the transport equations are solved using an explicit scheme. Flow and transport equations are solved sequentially at each time step. Besides, there is Chapter 3. A finite volume method for numerical simulation of the sharp interface model a general module for solving flow and transport problems, called "NpMc", where N is the number of phases and M is the number of components. Moreover, the choice of linear and nonlinear solvers, the type of grids as provided by the DUNE library, is maintained by the user in the problem file. Finally, we have the parameters, properties and physical law for each problem solved in the source files.

In this work, we have developped a new module, called 2p-SWI, to predict seawater intrusion problem using the sharp interface approach by means of the fully implicit approach. This method enable us to decrease in a significant way the CPU time while maintaining a good accuracy of the results. The nonlinear system is linearized by Newton's method, wherein each iteration a linear system is solved by a BiConjugate Gradient STABilized (BiCGSTAB) combined to the precondinor (ILU). Other iterative linear solvers and preconditioners are available in the DUNE library. For the approximation of the derivatives in the Jacobian matrix, we use numerical differentiation techniques. This is mainly to ensure the occurrence and efficiency of the implemented model. The management 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 simulations

This section aims to present numerical simulations of seawater intrusion problem using the sharp interface approach. Several 1D, 2D test cases are used to validate our method. Numerical results for each test case are presented. In the case of a confined aquifer, we have tested our module 2p-SWI against a linear analytical solution of Keulegan. The obtained results correspond well to the analytical solution and are in agreement with those proposed in [START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF]. After that, we have tested and validated our strategy in the case of an unconfined aquifer by proceeding to the study of the numerical convergence of the scheme. Furthermore, we examined three test cases: the first is proposed in [START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF] and the two others are realistic test cases with real data and complex geometry. The complete similarity of our results with those presented in [START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF][START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF] validates our approach and proves that our module is applicable and able to simulate seawater intrusion in coastal aquifers.

We note that all computations were performed on a laptop with Intel Xeon(R) CPU E3-1505M Processor (3.00 GHz) and 8 GB RAM.

Finally, a remarkable property of the TPFA scheme is that the discrete maximum principles (nonnegativity of the thickness of freshwater 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 [START_REF] Oulhaj | Numerical analysis of a finite volume scheme for a seawater intrusion model with cross-diffusion in an unconfined aquifer[END_REF].

Test 1: The rotating interface problem

To verify and validate the numerical model obtained from the TPFA method, the numerical results are compared with an analytical solution proposed in [START_REF] Keulegan | An example report on model laws for density current[END_REF]. This test consists to observe the movement of the interface without any external forces. In this case, we consider a confined aquifer of uniform thickness D with saltwater on the x < 0 part and freshwater on the x > 0 part. Homogeneous Neumann conditions are imposed at the borders. At the time t = 0, the interface starts to move due to the density difference α.

The elevation of the interface is therefore described by a linear profile passing through a fixed point (0, -D/2) in each time t ≥ 0. It is given by the following equation

Z(x, t) = - D 2 1 + x L(t) . (3.4.1)
The intersection of the interface with the base of the aquifer (toe the interface) is given for Z = -D by

L(t) = DD f t αφ . (3.4.2)
The initial position of the interface is given for L(t 0 ) = 20 m, position corresponding to t 0 = 12.28 days. The aquifer is approximated by a square ]-50, 50[×]0, 100[. The physical and geometric parameters used are those used in [START_REF] Essaid | A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: model development and application[END_REF][START_REF] Shamir | Motion of the seawater interface in coastal aquifers: a numerical solution[END_REF] and they are summarized in Table 3.1. To visualize the state of equilibrium of the salt front, a second simulation was performed with the same initial condition as the previous test. The interface between these two liquids rotates counterclockwise until it becomes horizontal as shown in Figure 3.5. The simulations described below permit to follow the rotation of the interface. We refer to Figure 3.6 for better visualization of the rotary character of the interface.

Parameters D [m] α φ D f [m/day] D s [m/day] S f , S s [m -1 ] Values 10 40
Here, we plot its evolution along the section x 2 = 50. The equilibrium state of the interface achieved at t = 800 days.

Injection and pumping processes

We consider the case of a free aquifer with soil compressibility. We associate homogeneous Neumann boundary conditions to our problem, which allows the solution to evolve freely. We examine two cases:

Injection of a significant amount of freshwater for a short period of 0.8 days through a well with a radius of 1.0 m centered at [START_REF] Bear | Seawater Intrusion in Coastal Aquifers, Concepts, Methods and Practices[END_REF][START_REF] Guo | User's guide to SEAWAT; a computer program for simulation of three-dimensional variable-density ground-water flow[END_REF].

Pumping freshwater through a well with a radius of 1.0 m centered at [START_REF] Bear | Seawater Intrusion in Coastal Aquifers, Concepts, Methods and Practices[END_REF][START_REF] Guo | User's guide to SEAWAT; a computer program for simulation of three-dimensional variable-density ground-water flow[END_REF]. We observe a local depression of the water table and the appearance of a salt cone " upconing" in the pumping area. The salt interface rises globally towards the horizon, which can be dangerous for the users. As to the injection scenario, the effects of freshwater filling with a decrease in the depth of salt water are observed, especially in the injection zone. model

The complete similarity of our results with those presented in [START_REF] Abudawia | Numerical analysis for a seawater intrusion problem in a confined aquifer[END_REF][START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF] validate our model.

The soil Compressibility

The purpose of this subsection is to visualize the effect of soil compressibility on the evolution of the groundwater surface Z T and the salt front Z. The higher the storage coefficient S f , the more water the soil can contain. Z when S f is zero. Usually, the injection process is used to decrease the salt cone and limit the marine intrusion problem. However, we remark the opposite phenomenon during a pumping process. Therefore, we see here that neglecting the compressibility of the soil can produce estimation errors on the evolution of the Z T and Z interfaces, particularly in the injection area.

Test 2: a field-scale free aquifer problem

We consider a test case presented in [START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF] which aims to ensure the applicability of the sharp-interface approach for an unconfined coastal aquifer subjected to pumping by comparison with dispersive approach results. The validity of the sharp interface approach is based on two parameters: the pumping rate and the position of the well. Numerical simulations for several scenarios have been considered. We briefly review the main results of this comparison. The sharp interface gives better results compared to the dispersive approach for higher pumping rates where the saltwater has reached the well screen. Moreover, it's produce good results when the pumping wells are deep and close to the coast. As a conclusion, for such problems, the sharp interface approach models well the saline intrusion.

To prove the performance of our developed module in the case of a free aquifer, we proceed to the study of the numerical convergence of the TPFA scheme. The efficiency and the accuracy are investigated through 2D simulations with different grid resolutions. After that, we compared our numerical results with those presented in [START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF] for several senarios. A good agreement between both results is observed.

Test description

We consider a hypothetical free aquifer of thickness 30 m and length 500 m. A Dirichlet condition for the saltwater head v = 30 m is imposed at the seaside boundary (x = 500 m). The constant freshwater flux 0.1 m 3 /day at land boundary (x = 0) is considered. Homogeneous Neumann boundary conditions are imposed on the rest of the boundaries. 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 .
The last quantities are calculated in the same manner proposed in [START_REF] Shi | Applicability of a sharp-interface model for estimating steady-state salinity at pumping wells-validation against sand tank experiments[END_REF]. One has,

Q f = D f l f D f l f + D s l s Q t , Q s = D s l s D f l f + D s l s Q t .
The thickness of the freshwater l f and the saltwater l s in the well screen respectively are calculated as follow

l f = max{min(Z wt , u) -max(Z wb , Z), 0}, l s = min{Z wt , u} -max{Z wb , Z B } -l f .
Z wt and Z wb are respectively the elevation of the top and the bottom of the well screen.

The freshwater head is noted by u and Z B is the elevation of the bottom of the aquifer. The physical well screen length is about 2 m while the diameter is set at 1 m. Other notations are presented in Figure 3.9.

Figure 3.9: Schematic of the well screen in an unconfined aquifer [START_REF] Shi | Applicability of a sharp-interface model for estimating steady-state salinity at pumping wells-validation against sand tank experiments[END_REF].

First of all, we have run the model under steady state conditions. In this case, the aquifer is assumed to be static (no pumping). The obtained results are used as initial conditions for the transient model (after pumping). The time required for the simulation is approximately 30 years. The properties and parameters of the aquifer are summarized in

Numerical convergence

In order to validate our implemented module, we proceed to the study of the numerical convergence. To do so, we select a scenario in the test case described in [START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF], for example Scenario 2. In this case, the pumping well is placed at a depth of 15 m and at a distance of x w = 350 m from the shoreline. We pump a constant amount of freshwater with Q t = 0.05 m 3 /day. We calculate the solution of the problem, corresponding to the saltwater and freshwater heads, for different grids. We have refined the initial mesh several times by the factor 2. The grids considered for the study of the numerical convergence are summarized in Table 3.4. The reference solution is calculated on a refined mesh, corresponding to 409600 cells with a resolution 3200 × 128 in the x 1 /x 2 direction with a very small time step equal to dt = 12 hours.

Grid levels

To visualize the convergence of the solution calculated on different grids to the reference solution, we make several sections in space and time. In Figure 3.10, we plot different sections in space from the point (0, 10), which represents the seashore, to the point (500, 10) at the landside. This section allows us to visualize the evolution of the solution, corresponding to the freshwater and saltwater heads and salt front, during the 30 years of activity. In Figure 3.10, we observe the convergence of the solution calculated on different grids towards the reference solution.

Section for freshwater head u.

Section for saltwater head v.

Section for the elevation Z. In the sequel, we visualize the convergence of the solutions in time. To do so, we select a point in the zone of pumping of freshwater corresponding to the point (350, 10). In Figure 3.11, we plot the curves of the freshwater head u, the saltwater head v and the salt front elevation Z on 30 years (9.4 × 10 8 s). The evolution of the solution calculated on different grids is compared to the reference solution computed on a fine mesh. We can see a good convergence of these quantities to the reference solution during 30 years of activity. model

Evolution in time for freshwater head u.

Evolution in time for saltwater head v.

Evolution in time for the elevation Z. As we are interested in the accuracy of the scheme, we are going to evaluate the error of the numerical scheme. We denote by er l 2 the discrete l 2 -relative norm. Moreover, we study the error between the reference solution and the numerical one in l ∞ and we note by er l ∞ the discrete l ∞ -relative norm. Let

er l 2 = i=N Ref i=1 |u h (i) -u ref (i)| 2 i=N Ref i=1 |u ref (i)| 2 and er l ∞ = sup i=1,...,N Ref | u h (i) -u ref (i) | sup i=1,...,N Ref |u ref (i)| . u ref (i)
is the reference solution and u h (i) is the numerical solution calculated on a certain grid at element i while N Ref represents the number of elements in the fine mesh. We present in Table 3.5 the discrete l 2 -relative norm and in Table 3.6 the discrete l ∞ -relative norm.

Number of cells 400 1600 6400 25600 102400 Freshwater head u 4.37e -3 2.05e -3 6.29e -4 3.01e -4 9.15e -5 Saltwater head v 2.00e -3 1.08e -3 3.77e -4 1.76e -4 6.28e Let p be the order of convergence of the scheme in space. Then, there exists a constant C such that: er l 2 (h) Ch p . To estimate the order of convergence numerically, we have refined the initial mesh several times by a factor of 2 for example. For two successive mesh sizes h and h/2, one has

er l 2 ( h 2 ) ≈ C h 2 p ≈ er l 2 (h) 2 p .
We obtain

p ≈ log er l 2 (h) er l 2 ( h 2 ) log(2) . (3.4.3)
The next step consists to calculate the order of convergence using equation (3.4.3). Tables 3.7-3.8 represent the order of convergence of the scheme calculated for u, v and Z in the l 2 and l ∞ relative norms. It is of order 1 for the primary variables, which seems logical considering the complexity and nonlinearity of the problem.

Number of cells

Numerical results

For the simulations of this test, we use a uniform rectangular mesh of 250 × 10 cells in the x 1 /x 2 direction for the control volumes. We performed the simulations with an initial time step of 0.01 s and a maximum time step of 1 day. The time step is incremented or decremented according to the number of iterations required by Newton's method in the previous time iteration. If Newton's method does not converge to a fixed maximum number of iterations, the time step is then divided by 2. In the opposite case, the time step is multiplied by two. This process runs until the maximum time step is reached.

The tolerances for the Newton's method and the linear solver BICGSTAB are respectively 10 -8 and 10 -6 . For this simulation, Newton's method converges quickly with less than 5 iterations while the CPU time required is inferior to 1 min. In the following, we will compare our developed model 2p-SWI implemented in DuMu X (Curves: SWI-DuMu X ) with the one presented in [START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF] (Curves: Mehdizadeh and al) while changing the pumping rate and well position. The well position is represented by a solid rectangular on Figure 3.12. The plots below illustrate that our numerical solutions lead to results similar to those obtained in [START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF].

a) Case 1: varying the pumping rates

The objective here is to visualize the effect of the quantity of water pumped on the variation of the saline front. In this case, the well screen is fixed and its central point is located at (x w = 350, x 2 = 0) with a diameter of 0.5 m. The depth of the well with respect to the horizontal surface is z w = 15 m. We simulate the saltwater interface with different pumping rates, which correspond respectively to 0.1 m 3 • d -1 (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 corresponding to each scenario are shown in Figure 3.12. throughout the pumping time, the salt interface gradually rises and the cone of the polluted water tends towards the pumping well. The progression of the cone is related to the quantity of the water pumped. The interface is too far from the base of the well in (Sc-2) which corresponds to the smallest amount of freshwater pumped. In this case, pumping do not present any risk of pollution. However, saltwater reaches the base of the well in (Sc-4). The latter is salinized by the seawater, which can be dangerous for the users. The interface is little far from the well bottom in (Sc-1, Sc-3) but the risk is still present. It can, therefore, be deduced that the salinization of the well increases with the increase in the quantity of water pumped.

b) Case 2: changing the depth of the well

In this test, we mainly concentrate on the effect of the well position together with the pumping rate on the variation of the sharp interface. The well is placed in the center, at the bottom and at the upper part of the aquifer. We consider two different pumping rates Q t = 0.1 m 3 /day, correspond to (Sc-1, Sc-5, Sc-6), and Q t = 0.07 m 3 /day for (Sc-3, Sc-10, Sc-11). We plot in Figure (3.13) the salt front for different position of the well screen. model As shown in Figure 3.13 (Sc-5) and Figure 3.14 (Sc-10), the well screen is placed at the bottom of the aquifer. We are talking about the full penetration of the pumping well where the salt cone has disappeared and the well is completely filled with saltwater. Steady-state conditions are reached quickly with Q t = 0.1 m 3 /day, where more time is needed for Q t = 0.07 m 3 /day.

For the other simulations, we have a partial penetration of the strainer. We observe the evolution of the salt interface and the appearance of the "upconing". The progress of the latter depends in particular on the location of the well . In (Sc-6), the polluted cone is too far from the base, which does not present any risk, while it is too close in (Sc-1). We can conclude that the deeper the well, the faster the saltwater reaches the base.

The same remarks can be made when we decrease the pumping rate to Q t = 0.07 m 3 /day but the effects are less than those obtained with Q t = 0.1 m 3 /day. -8) corresponds to the greatest distance from the coast compared to the other simulations, it can be seen that the saltwater is too far from the well position, whereas it is completely polluted in the nearest one ( Sc-9). The saltwater cone is too close in model (Sc-1, Sc-7) which means that the risk is higher in these cases. We conclude that if the pumping well is close to the coast, the probability of its pollution by seawater is higher.

In addition, the risk decreases as you move away from the coast.

Remark 3.4.1. We have presented a comparison of our results versus those obtained in [START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF] for a test case dealing with a homogeneous unconfined aquifer subject to eleven scenarios with different pumping rates and different well locations. We can observe a very good agreement between both calculations, ensures the validity of our developed model. From these series of simulations, it can be seen that if the well position is deep and close to the coast, the risk of pollution is higher, the same remarks can be made with a large quantity of water pumped.

We note that these simulations are carried out until steady-state conditions are obtained, which means that the aquifer will be exploited for many years. It can be noted, therefore, that these simulations can be useful to control the use of water resources and avoid pollution of the aquifer.

Test 3: Souss-Chtouka case study

In this subsection, we present numerical results obtained using the 2p-SWI module in the case of the Souss-Chtouka test case presented in [5]. We start by giving a short description of the Geographical location of the Souss-Massa basin. Moreover, we present different parameters corresponding to the geometry, boundary conditions, geological data and pumping rates, used for the prediction of saline intrusion in the Souss and Chtouka aquifer. The plain is under coastal groundwater conditions and assumed to be static before 1986. Therefore, the model is first run under steady state conditions. The comparison between the approximate and the measured piezometric map of 1968 showed a good agreement. The obtained results are used as initial conditions for the transition model. Different sections in space and time are used to visualize the saltwater intrusion under 80 years of over-exploitation of the plain. Numerical simulations have shown that the salt front elevation has moved to the landside and that the groundwater flow decreases in the whole plain. This means that our module is efficient and able to simulate the saltwater intrusion in a realistic test case.

Geographical location

The Souss-Massa basin is located in the south-western of Morocco (see Figure 3.16), which contains an important hydrological catchments with an area of 27000 Km 2 . The natural limits of the Souss-Massa river basin are: the Anti-Atlas mountains in the south, the High Atlas mountains in the north, the Siroua massif in the east, and the Atlantic ocean in the west (see Figure 3.16). In this region, significant groundwater tables have been identified occupying the Souss plain (4500 Km 2 ), the Chtouka plain (1260 Km 2 ) and Tiznit plain (1200 Km 2 ). 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 marl and conglomeratic intercalations of the down-land areas of the Souss plain and to fluvial-lacustrine deposits of the Souss unit extending to the Anti-Atlas chain border [5]. 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.

Physical 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 [5]. Ten zones have been recognized, for which the spatial hydraulic conductivity values vary between 1.21 m/day and 40 m/day. The value is higher in the lower part of the valley, especially in the north-west, where agricultural activity is very important, than in the upper part. Porosity values varying between 0.1 and 0.25 have been used, depending on the geologic material lithology. The porosity and hydraulic conductivity distributions are specified in Figure 3.19. For the initial conditions, the plain is supposed to be in hydrostatic equilibrium. Saltwater is static and freshwater is in dynamic. The problem is, therefore, run under steady-state Chapter 3. A finite volume method for numerical simulation of the sharp interface model conditions, using the parameters obtained from the calibrations done in [5]. The obtained results are compared with the measuring head of 1968 piezometric (see Figure 3.20 right).

The results are satisfactory and used as initial conditions for the transition regime. 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 3.10. 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 [5]. On the western boundary, at the contact with the ocean a zero head is imposed (u = v = 0 m).

Points

For all simulations, a triangulation mesh of 19, 520 cells is used for control volumes. We have ensured that the orthogonality condition is satisfied while constructing the mesh. This condition is necessary for consistent and efficient approximations of the TPFA method. Simulations were achieved with an initial time step of 10 s. The tolerances for the Newton method and the BICGSTAB solver are respectively 10 -8 and 10 -6 . In this case, Newton's method converges in less than 5 iterations. As expected, the timestep 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 5 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 sequel, we will present results corresponding to the coarse mesh.

Numerical results

In this subsection, we will present numerical results obtained using our 2p-SWI module for modeling seawater intrusion into the Souss-Chtouka aquifer. We consider the hydrogeological data and physical parameters described above. The model is, thereafter, run under transition condition over an 80-year period. Figure 3.22 illustrates the piezometry of freshwater in 1968 and after 80 years of exploitation. Up to 1968, the Souss-Chtouka aquifer was not subject to any exploitation. It follows the principle of the functioning of coastal aquifers. It is supposed to be in hydrostatic equilibrium. The aquifer is inflated and the equipotentials are regular and vary from 0 downstream to 100 upstream as shown in Figure 3.22 at left. However, for the 2048 predictions, equipotentials are no longer regular. Indeed, after 80 years of drastic exploitation, consisting of a general decline of the water table from the coast to the upstream ( u(t) = 100 -0.625t), the level of freshwater is decreased by 50 m in the upper part of the aquifer and a local depression of the water table in pumping wells position illustrated by negative values of the freshwater head.

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 (see Figure 3 Until 1968, the plain follows the principle of coastal aquifers in stable conditions. The only fluid in movement is freshwater. The saltwater level is supposed to be identified with the sea level, which is equal to zero. After 80 years of drastic exploitation, the aquifer has experienced a general decrease in seawater level from the coast to the whole basin illustrated by negative values.

Figure 3.25 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, Figure 3.25 shows that the salt bevel is stable and does not advance in the continent (in regions with low permeability). On the second section, (Figure 3.27 (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 3.27 (right)), the salt bevel is displaced laterally over a distance of at least 6K m 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 3.28 (right)), 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 3.28 (left)), shows almost no movement of the salt bevel despite a linear decrease in the free surface area, which is generalized to the entire basin. It should be noted that the region interested in logging is free of any pumping zones.

To understand the functioning of the Souss-Chtouka aquifer during the period of its activity, we propose to visualize the evolution of several quantities over time. To do so, we have selected five points on the surface of the aquifer. Figures 3.29 show the evolution of u, v and Z over time.

Points position

Elevation of salt front.

Freshwater head. Saltwater head. The results show that the salt front elevation is stable and did not move during the period of its exploitation at the second and third points. This seems physically logical since the corresponding area is free of any external exploitation. At the third and fifth points, the salt interface has moved about 150 m towards the ground during the 80 years of operation. The area concerned is under immense stress. In addition, the hydrogeological parameters (high permeability and porosity) improve the progression of the salt front. The third point is located at the position of the pumping well with a rate of 2686.99 m 3 /day. As a result, the interface has moved rapidly about 250 m to the surface, particularly over the past 50 years. This is the result of a combination of over-exploitation and the nature of the soil.

A general decline in freshwater levels is observed in Figure 3.29 (right) over the 80-years period. At the fifth point, the freshwater level decreases by 13 m during the pumping period. The evolution of the saltwater level is shown in Figure 3.29 (left). The saltwater level in the area subjected to stress decreases in the first few years and then stabilizes again.

Remark 3.4.2. The numerical results of the "Souss-Chtouka" test case proved that the new module developed is capable of providing an accurate solution and can predict the location, shape and extent of the water table and of the the freshwater/saltwater interface in coastal aquifers under a variety of different stress conditions, thereby illustrating its robustness with a sufficient convergence rate.

Chapter 3. A finite volume method for numerical simulation of the sharp interface model

To end this section, we have 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 heterogeneous fields. Furthermore, the fully coupled fully implicit scheme greatly reduce the CPU time of the simulations.

Test 4: Tripoli aquifer case study

In this section, we present the numerical results obtained in the case of the Tripoli test case proposed in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF] using the sharp interface module. It deals with the impact of demographic evolution on seawater intrusion in the lower Tripoli aquifer in Lebanon. We present a brief description of the geographical location of the Tripoli aquifer and the different physical parameters employed in the simulations. To ensure the validity of our developed module in a real case, we propose to compare our numerical results with those presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF]. Different scenarios are considered and presented based on the available hydrological data. Numerical simulations have shown a reduction of the water potential and an advance of the freshwater/saltwater interface.

Geographical location

The city of Tripoli is the second pole of attraction located in northern Lebanon (see Figure 3.30), wich covers a total area of 400 Km 2 . It is located between the Mediterranean Sea (about 15 Km) in the west, the Mountain of "Bcharré" in the east and the Mountains of Qalhat and Torbol from the south and the north, respectively. The groundwater is mainly supplied by precipitation and snow-melt from the "Bcharré" region located at the east of the Tripoli aquifer (see Figure 3.30).

Figure 3.30: Position of the Tripoli aquifer, study zone and wells position (cf. [START_REF] Kalaoun | Modélisation de l'intrusion saline dans l'aquifère de Tripoli-Liban: impact des changements globaux[END_REF]).

In this work, we consider only the lower Tripoli aquifer. The bedrock of the aquifer is the impermeable layer "C6" formed by marl and marl limestone (See Figure 3.31). The thickness of the aquifer above the impermeable layer is about 200 m.

Figure 3.31: Geological map and stratigraphy of Tripoli (cf. [START_REF] Kalaoun | Modélisation de l'intrusion saline dans l'aquifère de Tripoli-Liban: impact des changements globaux[END_REF]).

The studied domain and physical parameters

The study area is approximately 3 Km × 6 Km. The Tripoli aquifer is an homogeneous porous medium with a thickness of 200 m. The geometry and boundaries of the aquifers are given in Figure 3.32. For the mesh, we use a triangulation mesh of 37568 elements an 10285 vertexes.

The boundary conditions are divided into two parts: Γ 1 represents the shoreline of the aquifer and Γ 2 coresponds to a fracture separating the lower and the upper zone of the city. A Dirichlet boundary condition (u = v 0 = 200 m) is used on Γ 1 while a freshwater inflow Q inf low is imposed on Γ 2 . This flux has been calibrated in [START_REF] Kalaoun | Modélisation de l'intrusion saline dans l'aquifère de Tripoli-Liban: impact des changements globaux[END_REF] for the year 2008. model The physical parameters used in the simulation are listed in [START_REF] Kalaoun | Modélisation de l'intrusion saline dans l'aquifère de Tripoli-Liban: impact des changements globaux[END_REF]) .

The direction of the flow in the Tripoli plain can be assumed mainly horizontal since its thickness (200 m) is negligible compared to its horizontal surface (18 Km). This ensures the validity of the Dupuit approximation and allows us to vertically integrate the freshwater and saltwater equations in 3D.

Due to the increase in population over the past decades in the city of Tripoli, the high demand for drinking water threatens the freshwater potential. Figure 3.33 represents the pumping areas as described in [START_REF] Kalaoun | Modélisation de l'intrusion saline dans l'aquifère de Tripoli-Liban: impact des changements globaux[END_REF][START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF]. Private wells are located in the center of Tripoli city and represented by a blue rectangle in Figure 3.33. In addition, six active public wells maintained by the Tripoli Water Authority are considered and presented by yellow markings in Figure 3.33. Four wells are located in the village of "Bahsas" and two others in the neighborhoods of "El-Jisr" and "Malouleh". 3.12: Wells position and associated rates.

Four scenarios are used to predict the extent and progression of the salt wedge in the lower Tripoli aquifer:

"NOP": represents the no pumping scenario and corresponds to the results of the year 2008. In this case, the plain is assumed to be in hydrostatic equilibrium and the model is run under steady state conditions.

"REC": refers to a scenario with 150 L/day as required by the water authorities after 10 years of exploitation.

"10 Y": refers to a scenario with the current pumping rate of 250 L/day for 10 years of exploitation.

"20 Y": refers to a scenario with the current pumping rate of 250 L/day after 20 years of exploitation.

Numerical results

In this subsection, we will present numerical reults obtained using our developed module to predict seawater intrusion into Tripoli case study. We performed the simulations with an initial time step of 0.01 s and a maximum time step of 20 days. The tolerances model for the Newton's method and the linear solver BICGSTAB are respectively 10 -8 and 10 -6 . For this simulation, Newton's method converges quickly with less than 4 iterations while the CPU time required is less than 13 min. In the following, we will compare the result obtained with our developed model 2p-SWI implemented in DuMu X with the one presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF] in differents situations. The comparison between our numerical results and those presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF] showed good agreement between both calculations.

We consider the hydrogeological data and the physical parameters described above. The model is run under steady state conditions. The obtained results are labeled as "NOP" pumping scenario and used as initial conditions for the other scenarios. Figure 3.34 presents the evolution of the freshwater/saltwater interface in 2008 in the natural conditions.

Presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF]. Implemented in DuMu X . Presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF]. Implemented in DuMu X . Presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF]. Implemented in DuMu X . We illustrate in 3.38 the piezometric head of freshwater in natural conditions. We can see that the equipotentials are regular and vary from 200 m at the coast to 206 m at the landside. In Figure 3.39, we illustrate the position of the freshwater head u in different situations. In Figure 3.39a), the pumping rate in private wells as recommended by the water authorities is about 150 L/day. We can observe that the water potential is stable and has not changed during 10 years of activity while it decreases with a rate of 250 L/day for the same period. During 20 years of over-exploitation with a rate of 250 L/day, the freshwater level is decreased in the whole plain and an exceptional depression is located around the position of the pumping wells.

Presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF] Implemented in DuMu X . In this section, we presented a comparison between our results and those presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF] for a homogeneous test case dealing with the impact of demographic change on seawater intrusion in the lower Tripoli aquifer in Lebanon. We can see good agreements between the two calculations, which ensures the validity of our developed model. The numerical results of the test case "Tripoli aquifer" showed the ability and capacity of our model to predict the seawater intrusion in real scale test cases.

Conclusion

In this chapter, we have developed a fully coupled, fully implicit approach of the sharp interface seawater intrusion problem. The first part deals with the numerical scheme of the mathematical model. The coupled system is discretized by a cell-centered finite volume method in space. For the time discretization, we used an implicit Euler method, which allows us to take large time steps. The CPU time is, therefore, significantly decreasing while maintaining a good accuracy of the scheme. The nonlinear system is solved by Newton's method, wherein each iteration a linear system is solved by a preconditioned gradient BiConjugate STABilized (BiCGSTAB). The time step is chosen according to the number of iterations required by the Newton method to achieve convergence in the last iteration. Thereafter, We have developed and implemented a new module, called 2p-SWI, in the context of the parallel open-source platform DuMu X , based on DUNE. A short description of the platform is also provided. Afterwards, the 2p-SWI module has been validated in several examples.

The first test case is described in [START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF]. Numerical simulations are compared to an analytical solution proposed in [START_REF] Keulegan | An example report on model laws for density current[END_REF]. The obtained results are satisfactory and in good agreement with those in [START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF]. Next, we have considered the field-scale free aquifer test case [START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF]. Our numerical simulations are compared to those in [START_REF] Mehdizadeh | Assessment of sharp-interface approach for saltwater intrusion prediction in an unconfined coastal aquifer exposed to pumping[END_REF]. The numerical results demonstrate the accuracy and robustness of our approach.

After that, we have applied our model to the Souss-Chtouka field case located in the southwestern of Morocco. The numerical model is based on available hydrogeological data in real scale by simulating the progress of seawater intrusion during 80 years. Until 1968, the aquifer is supposed to be in hydrostatic equilibrium. To do so, the model is first performed under steady-state conditions. The obtained results constitute, therefore, the initial values of the transition model. After 80 years of exploitation, the salt front has moved to the landside. Also, the level of the freshwater head decreases in the entire basin.

The numerical results show that marine intrusion is dependent on several factors, such as the amount of extraction carried out and the location of wells close to the shoreline. It should also be noted that marine intrusion is favored in areas of high permeability, even in the absence of pumping. Afterwards, we validated our implemented module on a real test case corresponding to the Tripoli aquifer in the northeast of Lebanon. Numerical simulations over 20 years of exploitation are presented. The comparison of our numerical results and those presented in [START_REF] Kalaoun | Impact of demographic growth on seawater intrusion: Case of the Tripoli aquifer, Lebanon[END_REF] showed a good agreement, which ensures the validity of our module to provide efficient and accurate results in real scale test cases.

Introduction

In this chapter, we present a finite volume scheme and numerical simulation of the sharpdiffuse interfaces model into coastal aquifers, which was recently introduced in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF][START_REF] Choquet | Mathematical analysis of a seawater intrusion model including storativity[END_REF]. We present the derivation of the mathematical model in section 4.2. Such flows are governed by a coupled system of two parabolic partial differential equations describing a two-phase immiscible seawater/freshwater flow and integrating the dynamics of transition zones.

The non-degeneracy of the equations comes from the presence of the term modeling the transition zones.

In section 4.3, we developed a fully coupled, fully implicit finite volume scheme to discretize the governing equations. Precisely, we employ a TPFA method in space on a structured mesh where the orthogonality condition is satisfied. We apply an implicit Euler scheme for the time discretization, which allows us to take large time steps and thus reduce the CPU time. An upwind scheme is used to approximate the convective flux. In this context, we have developed and integrated a two phase diffuse module in the framework of DuMu X . The validity of the implemented module is examined through the study of the numerical convergence of the scheme. Moreover, a 2D simulations on different grid resolution prove the ability of the developed module to predict the position of the interface in the presence of the transition zones. This has been discussed and presented in section 4.3.

Thereafter, we propose to compare the numerical results of our 2D module with those obtained with the classical 3D model for miscible displacements. The 2D model yields information on the position of the salt front whereas the 3D model describes the evolution of the salt concentration. Numerical simulations of both models are afterwards compared to those presented in [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF]. The numerical results showed that this approach gives physically consistent and performance results. Finally, section 4.4 presents some concluding comments.

Mathematical model of the sharp-diffuse interfaces problem

In this section, we present the mathematical modeling of seawater intrusion using a mixed sharp-diffuse interfaces approach [START_REF] Choquet | Mathematical analysis of a sharp-diffuse interfaces model for seawater intrusion[END_REF][START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. Let us mention that the derivation of this model is done in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. Therefore, we only present here the outlines to obtain the governing equations modeling saltwater intrusion using the mixed approach. Freshwater and saltwater fluxes are modeled by mass conservation and Darcy's laws for each fluid. By neglecting the vertical variation of the fluids, the flow is assumed to be horizontal. The problem is therefore reduced to a 2D model that takes into account the dynamics of the transition zones. Neglecting the saltwater storage coefficient, the governing equations modeling the freshwater and saltwater flows are as follows

S f b f ∂u ∂t -div (b f D∇u) + F (u) -F (Z) = Q f , (4.2.1) -div ((1 + γ)b s D∇v) + F (Z) = Q s . (4.2.2)
where u is the freshwater head, Z is the salt front elevation and v is the saltwater head. The other notations are summarized in For modeling the transition zones, a phase field model developed in [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening[END_REF] is employed. For more details on the approximation of the fluxes across the transition zones, the authors can refer to [START_REF] Choquet | Mathematical analysis of a sharp-diffuse interfaces model for seawater intrusion[END_REF]. The expression for the approximate fluxes describing the dynamics of diffusive transition zones is as follows:

F (u) = φβ ∂u ∂t -δdiv (D∇u) . F (Z) = φ ∂Z ∂t -δ∆Z .
where φ and δ are respectively the porosity of the medium and the width of the transition zones. The elevation of the salt front is written γZ = (1 + γ)v -u, while the parameter β is used to combine the two cases: free and confined aquifers, β = 1, if the aquifer is free, 0, if the aquifer is confined.

Let Ω be an open polygonal domain, ∂Ω = Γ, x = (x 1 , x 2 ) ∈ R 2 and ]0, T [ a time interval.

Let Ω T =]0, T [×Ω and ∂Ω T =]0, T [×∂Ω. Summing equations (4.2.1)-(4.2.2) and invoking the expression of the approximated fluxes, the mathematical system is given by

                 (S f b f + βφ) ∂u ∂t -div (D(b s + b f )∇u) -div (γDb s ∇Z) -βδdiv (φD∇u) = Q f + Q s in Ω T , φ ∂Z ∂t -div (Db s ∇u) -div (γDb s ∇Z) -δdiv (φD∇Z) = Q s in Ω T , D∇u • n = 0, D∇Z • n = 0, Db s ∇u • n = 0 on ∂Ω T , D(b s + bf )∇u • n = 0, Db s ∇Z • n = 0 on ∂Ω T , u(., 0) = u 0 , Z(., 0) = Z 0 in Ω,
where the main unknowns are, from now on, the freshwater head u and the salt front elevation Z. The initial conditions u 0 (x) and Z 0 (x) are given for u and Z at (t = 0) to describe the initiale state of the aquifer while n is the unit outward normal to ∂Ω. The other physical parameters used in the system are defined above. interfaces model

The problem consists of a coupled system of two nonlinear parabolic partial differential equations simulating two immiscible fluids and tacking into account the dynamics of the transition zones.

Numerical simulations of the sharp-diffuse interfaces model

Let us mention that the finite volume numerical scheme for the sharp-diffuse interfaces model using the TPFA finite volume method is similar to that presented in chapter 3. Therefore, this subsection discusses numerical simulations of this model using the TPFA method to predict seawater intrusion in coastal aquifers.

The numerical scheme afterwards is implemented and integrated in the framework DuMu X . We validate our scheme on different test cases. First, we procceed to the numerical convergence of the scheme in the case of a free aquifer. Besides, we have tested our 2pdiff-SWI module against an analytical solution derived from Ferris model [START_REF] Ferris | Cyclic fluctuations of water level as a basis for determining aquifer transmissibility[END_REF]. The obtained results match well with the analytical solution and are in good agreement with those presented in [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF][START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF]. Our results are completely equivalent to the results given in [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF][START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF] thus validate our methodology and demonstrate the usability and capacity of our module to simulate marine intrusion into coastal aquifers.

Test 1: Pumping of freshwater

Here, we consider a test case described in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] dealing with the effect of pumping on the displacement of the salt front in an unconfined aquifer where diffusion in the transition zones is not neglected. First, we proved the validity of our sharp-diffuse interfaces model implemented in DuMu X by studying the numerical convergence of the TPFA scheme. Different meshes resolution are used to visualize the convergence of the solutions to the reference solution. Numerical results, proving the efficiency and accuracy of our model to predict saltwater intrusion, are presented.

Test description

We consider a free aquifer of thickness 10 m and horizontal surface Ω = ]-50 m, 50 m[ × ]-20 m, 20 m[. Homogeneous Neumann boundary conditions are imposed on the boundaries to let the interfaces evolve freely. We pump a constant amount of freshwater with different rates ranging from 0.8 m 3 /day to 0 according to the position (x 1 , x 2 ) for three days. The main pumping point is located at (x 1 = 15 m, x 2 = 0 m) with a maximum rate of 0.8 m 3 /day. The rate decreases as it approaches the boundaries. .

Numerically, the pumping (or injection ) is modeled by a term source imposed in a specified part of the domain. In this case, let the term source of the freshwater equation

Q f = -0.8 exp -0.01 (x -15) 2 + y 2 .
The initial position of the saltwater/freshwater interface is computed with an analytical solution proposed in [START_REF] Keulegan | An example report on model laws for density current[END_REF]. The position of the intersection of the interface with the base of the aquifer is initially set at 30 m.

The physical parameters and properties of the aquifer are summarized in

Numerical convergence

First, we proceed to the study of the numerical convergence. We calculate the solution of the problem for different cells. The initial mesh is refined several times. For example by a factor 2. We consider the following grids as presented in Table 4 We calculate the reference solution on a fine grid corresponding to 40960 = 320 × 128 cells. We start by making several sections in space and time. In Figure 4.2, we plot the section in space from the point (-50, 0) to the point (50, 0). The section passes through the main part of the pumping. Figure 4.2 shows the evolution of several quantities in a free aquifer subjected to pumping for 3 days. We observe the convergence of the solution calculated on different grids towards the reference solution.

In addition, we consider a plot selection over time taken at point (-10, 0) that belongs to the main part of the pumping. Figure 4.3 shows the evolution of the freshwater head, the saltwater head and the elevation of the salt front over time. We see that during 3 days of pumping, the solution calculated on different grids converges towards the reference solution. interfaces model Finally, we make an analysis of the error starting from an initial value 40 cells, then we will calculate the numerical convergence order. The tables 4.3-4.5 below represent the standard l 2 and l ∞ relative norms as well as the order of convergence. Cells 40 160 640 2560 10240 Freshwater head u 1.05e -2 4.92e -3 2.4e -3

1.16e -3 5.22e -4 Saltwater head v 8.9e -3 4.23e -3 2.07e -3 1.01e -3 4.52e -4 Elevation Z 6.03e -2 2.96e -2 1.46e -2 7.13e -3 3.18e -3 By refining the initial grid several times, the program takes more time to complete the process. The CPU time therefore gradually increase.

Simulations have been performed with an initial time step of 10 s and a maximum time step of 0.5 day. The time step size increases if Newton's method converges with less iterations in the last time iteration and decreases otherwise. This process is manipulated during the simulation according to the function SuggestStepTimeSize() implemented in DuMu X . The tolerances for the Newton method and the BICGSTAB solver are respectively 10 -8 and 10 -6 . In this test, Newton's method converges in less than 4 iterations.

Numerical results

This subsection aims to illustrate the displacement of the interface of the salt front and the ground water level in a free aquifer subjected to pumping during 3 days. In this test, we consider an uniform rectangular resolution of 80×32 in the direction x 1 ×x 2 . In Figure 4.5, we can see that the "upconing" is appeared and grow up quickly and a local depression is located at the main zone of pumping. The presence of the freshwater/saltwater transition zone and the saturated/unsaturated zone promotes diffusion and advancement of the free surface and salt front.

The velocity of freshwater and saltwater flows during 3 days of pumping is illustrated in Figure 4.6. The flow velocity of both fluids is maximum in the main pumping zone. Moreover the direction of the freshwater flow tends downward while tends to upper for the saltwater fluid. Furthermore, the direction of freshwater flow tends downward, corresponding to the decrease in freshwater level, while the saltwater flow tends toward the upper surface. This proves the ability of our model to predict saltwater intrusion into an unconfined aquifer and its ability to provide physically satisfactory results.

Chapter 4. A finite volume method for numerical simulation of the sharp-diffuse interfaces model

Test 2: Tidal effects

The purpose of this test is to assess the validity of our implemented module in the case of a free aquifer. For this purpose, we consider a test case presented in [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF][START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF], which mainly focuses on the interface displacement coupling with tidal fluctuation effects. We propose to compare our numerical results with an analytical solution derived from the Ferris model [START_REF] Ferris | Cyclic fluctuations of water level as a basis for determining aquifer transmissibility[END_REF]. The obtained results are promising and show a good agreement for both calculations. Afterwards, we compare our obtained results to those presented in [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF][START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] in different situations. The numerical simulations prove that our module is efficient and able to predict sewater intrusion in a free aquifer taking into account the width of the transition zones and sea fluctuations effects.

Test description

We consider a free aquifer represented by ]0, 50[×]0, 50[. The depth of the aquifer is fixed at 10 m. We combine the problem, of rotary character, with tidal effects. The oscillations of the sea progressively generate pressure waves close to the aquifer. For this test case, we use the parameters proposed in [START_REF] Cooper | A hypothesis concerning the dynamic balance of fresh water and salt water in a coastal aquifer[END_REF] after a scaling of our aquifer. We take D = 39.024m/day for the hydraulic head of freshwater and φ = 0.3 for the porosity. We impose a Dirichlet boundary condition at the seashore, which corresponds to {x 1 = -50}, for the salt front Z (while the freshwater head u evolves freely on the impermeable side). The Dirichlet value is calculated using an analytical solution proposed in [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF][START_REF] Ferris | Cyclic fluctuations of water level as a basis for determining aquifer transmissibility[END_REF]. Homogeneous boundary conditions are imposed on the rest of the boundary. The amplitude of the transition zone δ is evaluated from the equation for the width of the dispersive zone (eq. ( 2) in [START_REF] Cooper | A hypothesis concerning the dynamic balance of fresh water and salt water in a coastal aquifer[END_REF]). Its value is approximately δ = 0.1. We performed the simulations with an initial time step of 0.01 s and a maximum time step of 1 day.

Numerical results

In order to ensure the validity of our developed module in the case of a free aquifer, we propose to compare our numerical results issued from our developed module with an analytical solution computed with Ferris model. In this case, the aquifer is not subjected to any external forces. For these simulations, we consider two uniform rectangular mesh with mesh 1 = 100 × 100 and mesh 2 = 200 × 200 in the x 1 × x 2 direction. Figure 4.7 illustrate the position of the interface obtained with this present model compared to the analytical solution after a period and a half-period sea fluctuations.

T =0.5 day, mesh 1 T = 1 day, mesh 1 T =0.5 day, mesh 2 T =1 day, mesh 2 By refining the mesh twice, we can observe the convergence of the numerical solution to the analytical solution derived from the Ferris model at both times T = 0.5, 1 days. After one day, the interface evolves freely on the landside and returns to its initial position on the seaward side. The wave pressure causes oscillations on the seaside and accelerates the movement of the interface in the whole domain. Next, we propose to compare our obtained results with those presented in [START_REF] Abudawia | Analyse numérique d'une approximation élément fini pour un modèle d'intrusion saline dans les aquifères côtiers[END_REF] in two different cases. The objective here is to evaluate the ability of our module to handle different physical phenomena. For the spatial discretization, we use the resolution of mesh 1. We illustrate in Figure 4.8 the evolution of the interface at different times without any external force while Figure 4.9 presents its position during a pumping scenario. interfaces model We place a hydrograph artifacts at x = -35 to illustrate the behavior of the salt interface over time. Figure 4.10 shows the results predicted by the hydrograph artifacts for the sharp interface model (δ = 0.0) and the sharp-diffuse interfaces model (δ = 0). The obtained results for both models are then compared to the analytical solution derived from Ferris. Afterwards, our numerical results for both models (sharp and sharp-diffuse interfaces model) are compared to those presented in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. The Dirichlet value at the seaside evolves over time and take into account the tidal effect. The impact of tidal fluctuations is illustrated in both situations. The interface evolves freely as in natural conditions despite the oscillations at the seashore. In addition, pumping effects are observed at the well position.

Implemented in DuMu X Choquet and al [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] Figure 4.10: Oscillations on the salt front elevation, predicted by the hydrograph artifacts (δ = 0 (left) and δ = 0(right)) and comparison between both models.

The first comment we can make is the similarity of our results implemented in DuMu X (left) with those presented in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. The oscillations are more noticeable in the model with (δ = 0). Comparing the two models to the analytical solution derived from the Ferris model, we observe that the model with a sharp-diffuse interfaces matches the analytical solution well. The evolution of the freshwater head u in this case for both models is illustrated in Figure 4.11.

Now, we concentrate on the behavior of the freshwater head when setting up seaside oscillations and with Neumann homogeneous BCs for both variables u and Z. In Figure 4.12, we present the result presented by the hydro-graph over time. interfaces model Implemented in DuMu X Choquet and al [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. Implemented in DuMu X Choquet and al [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. Figure 4.12: Evolution of the freshwater head u with and without oscillations in the boundary.

We distinguish here two cases: with and without oscillations on the boundary. By imposing the Dirichlet value calculated with the Ferris model, we can observe sinusoidal oscillations caused by tidal effects at the hydrograph located at x 1 = -35. Using homogeneous Neumann boundary conditions, we observe that the water table level decreases and stabilizes after 0.25 days.

Test 3: 3D variable density test

We consider a test case described in [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF] that aims to numerically compare the 3D variable density model with the 2D sharp-diffuse interfaces model. The 3D model describes the evolution of the salt concentration while the 2D model gives the position of the salt front and takes into account the width of the transition zone. For the validation of our 2D and 3D models implemented in DuMu X , we propose to compare our numerical results obtained for both models with those presented in [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF]. The CPU time required for the 2D simulations is less than one minute, while it takes about 2 hours and 30 minutes for the 3D model.

Numerical results

In this sub-section, we will present numerical results for each model (3D and 2D models) implemented in DuMu X compared to those of [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF]. We consider the data and the physical parameters described above. Different cases dealing with saltwater pollution of an aquifer are analyzed.

Case 1: rotating interface problem

In this test, we focus on the behavior of the interface at normal conditions where the aquifer is not subjected to any external forces (i.e., Q 1 = Q 2 = 0.0 ). To do so, we consider that the boundaries of the domain are impermeable. The storage coefficient S f is neglected and we fix the density contrast γ = 0.025. The figure 4.15 shows the evolution of the salt concentration and the salt front obtained by the 3D and 2D models.

Cherfils and al [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF]. 3D density variable.

Sharp-diffuse interfaces model 2D. The first remark we can make concerns the consistency between our results and those presented in [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF]. The simulations start with a quasi-horizontal interface and, under the density contrast the interface evolves freely to reach its equilibrium state after 10 days.

On the other hand, we can see that the diffuse interfaces issued from the 2D model is consistent with the 50 % contour of the transition zone in the 3D model. This can be confirmed with the result presented in figure 4.16. The dotted lines indicate 5 and 95 % salinity contours in the 3D model while the solid line represents the salt front elevation issued from the 2D model. We can see a good agreement between the results issued from both models. We illustrate in Figure 4.17 the gradient of the freshwater head. In Figure 4.19, the upconing shown by the sharp-diffuse 2D model is not a pike as it is known for the classical sharp interface model. This is due to the existence of the width of the transition zone in the system of equations.

Case 3: injection process

In the following experiment, we add some form of injection to the problem. This is numerically modeled by a positive value imposed in a specific part of the domain. In our case, suppose that Q 2 = 0.7χ The pumping effect is visualized in both models. Two upconings are formed and tend to the surface under the high pumping rate. The 2D model is able to visualize the pumping effect as does the 3D model and the transition zone is large due to the force applied in this region.

Case 5: pumping and injection experiment

For the last test, we combine the effect of injection and pumping at the same time. We choose a production well at the left and an injection well at the right. The source term in this case is given by Q 1 =0.0 and The pumping effect does not appear too much and a stabilization of the interface in the whole domain is observed due to the injection process. Also in this case, we can see that the 2D model yields physically efficient and accurate results.

Q 2 = -1 • χ [-5.
Remark 4.3.1. We have compared our results with those presented in [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF] for both 2D and 3D modes for a test case handling the contamination of a confined aquifer by saltwater in different situations. The obtained results proved a good agreement for both models, which ensures the validity of our methodology. This comparison is used to show the ability of the 2D model to predict seawater intrusion in different cases. The numerical results proved that the sharp-diffuse interfaces 2D model provides physically efficient results and take into account the width of the transition zone.

Conclusion

This chapter deals with a finite volume numerical method for modeling seawater intrusion problem using the mixed sharp-diffuse interfaces approach in coastal aquifer. First, we outlined the main steps of the derivation of the mathematical model. The construction of the numerical scheme is established by integrating the governing equations on each control volume and evaluating the flows on each edge. We adopt a fully coupled, fully implicit TPFA finite volume method for the discretization of these equations in space while we use an implicit Euler scheme for the discretization in time. An upwind scheme is used to approximate the convective fluxes.

Afterwards, the numerical scheme is integrated and implemented in the framework DuMu X . The developed module has been validated by performing a numerical convergence of the scheme. The numerical simulations showed that the existence of a diffuse interfaces model freshwater/saltwater and saturated/unsaturated zones favors the advance and diffusion of the salt front and water surface level. Moreover, we tested the ability of our developed module to take into account physical and hydrogological phenomena such as tidal fluctuations. To do so, we consider a test case described in [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF][START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] that aims to illustrate the impact of tidal fluctuations caused by waves motion on the displacement of interfaces. Our numerical results are compared to an analytical solution derived from Ferris model and showed good agreement with the results presented in [START_REF] Abudawia | A finite element method for a seawater intrusion problem in unconfined aquifers[END_REF][START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] in different cases.

Besides, the validity of our module is examined by comparing our results obtained for the 2D model with the classical 3D model for miscible fluid flows. Our results for both models (2D and 3D) are then compared with those presented in [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF] and presents good correspondence.

From these series of simulations, it can be seen that the 2D sharp-diffuse interfaces model can provide an efficient and accurate results and take into account various hydraulic and physical phenomena with a reasonable CPU time. The numerical analysis of an extension of this method to a more efficient and accurate scheme to deal with heterogeneity and anisotropy on an unstructured mesh will be discussed in the next chapter.

Introduction

This chapter aims to present the convergence of the finite volume Multi-Point Flux Approximation (MPFA) scheme for a sharp-diffuse seawater intrusion problem. The derivation of the mathematical model is presented in Chapter 4. Here, we consider a sharpdiffuse interfaces model in an anisotropic and heterogeneous free aquifer. Such a process 90 Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model is governed by a coupled system of two parabolic partial differential equations simulating respectively the freshwater head and the salt front elevation. To handle such a system numerically, we propose a fully implicit MPFA finite volume approach that solves the governing equations at each time step using Newton's method. In fact, this method tackles the non-linearity of the problem, the anisotropy and the heterogeneity of the aquifer on unstructured mesh, which improves the efficiency and accuracy of the numerical computations of such system.

The present chapter has several objectives. First, we will develop a fully coupled fully implicit approach which combines advantages of the MPFA method to accurately solve the diffusion terms with an upwind method for space discretization on unstructured grids.

The second objective of this paper is to present a proof of the convergence of this scheme which is established due to some compactness properties. Moreover, we have demonstrate the weak convergence of the discrete gradient of the freshwater head and the salt front elevation. Passing to the limit, we have proved that the approximate solution is a weak solution to the continuous problem. Lastly, we have developed and implemented this scheme in a new module in the context of the parallel open source platform DuMu X [START_REF]DuMu X DUNE for multi-{ Phase, Component, Scale, Physics[END_REF], based on the Distributed and Unified Numerics Environment (DUNE) [2], allowing simulations for large-scale field applications involving seawater intrusion in coastal aquifers.

The overall objective of this part is the development of a new-generation framework and reservoir simulator suitable for massively parallel processors. Some numerical examples are presented, one of which is related to flows in a fractured porous medium.

The outline of the paper is as follows. Section 5.2 is devoted to the statement of the problem. Namely, in this section, we introduce the system of equations for the sharp-diffuse interfaces model in a free coastal aquifer [START_REF] Choquet | Mathematical analysis of a sharp-diffuse interfaces model for seawater intrusion[END_REF][START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]. Then we formulate the main assumptions on the data. In Section 5.3 we recall the general FV framework and formulate the fully coupled fully implicit MPFA FV scheme that will be studied. A discrete maximum principle, that will be used in the energy estimates for the scheme, is proven in Section 5.4. In Section 5.5.2, we first establish some energy estimates followed by the existence of discrete solution to the FV scheme. In Section 5.5.3, we establish some compactness results corresponding to space and time translate estimates. In Section 5.6, the convergence of the scheme is proved by compactness arguments. In Section 5.7, a description of the implementation of our strategy in DuMu X is given. Then, to validate our approach, we consider two test cases. The first test case is a simulation in a homogeneous domain, while the second one is in a heterogenous free surface aquifer where two intersecting fractures are considered. Finally, concluding remarks are summarized in section 5.8.

Mathematical model

We consider an open bounded domain Ω of R 2 , describing the projection of the aquifer on the horizontal surface, over a time period ]0, T [. Let Ω T =]0, T [×Ω, ∂Ω T =]0, T [×∂Ω and n the unit outward normal to ∂Ω. We consider the sharp-diffuse interfaces model for seawater intrusion problem in coastal free aquifers describing two immiscible phase seawater/freshwater flow and tacking into account the width of transition zones. Then the mathematical model is given by a coupled nonlinear system of parabolic PDEs (see [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF]) where the main unknowns are the freshwater head u [m] and the salt front elevation

Z [m]: (P )                  φ ∂u ∂t -div (D∇ϕ(u)) -div (D∇ϕ(Z)) -div (φδD∇u) = 0 in Ω T , φ ∂Z ∂t -div (Db s (Z)∇u) -div (D∇ϕ(Z)) -div (φδD∇Z) = 0 in Ω T , D∇u • n = 0, D∇Z • n = 0, Db s (Z)∇u • n = 0 on ∂Ω T , D∇ϕ (u) • n = 0, D∇ϕ (Z) • n = 0
on ∂Ω T , u(0, .) = u 0 (.), Z(0, .) = Z 0 (.)

in Ω.

(Z) = Z -Z B [m] the thickness of the freshwater, ϕ(u) = 1 2 (u -Z B ) 2
, u 0 and Z 0 are given functions. We have to mention, that for simplicity, we have neglected the source terms and taken γ = 1.

The main assumption on the data are as follows:

(A.1) The tensor D is assumed to be bounded symmetric and uniformly elliptic in Ω.

There exists constants

D + ≥ D -> 0 such that 0 < D -|ξ| 2 ≤ D(x)ξ • ξ ≤ D + |ξ| 2 < ∞ for any ξ ∈ R 2 and ξ = 0, a.e. in Ω. (A.
2) The porosity belongs to L ∞ (Ω), and there exists two positives values φ -and φ + such that 0 < φ -≤ φ(x) ≤ φ + < ∞ a.e. in Ω.

(A.3) u 0 and Z 0 are in L 2 (Ω) and the initial thicknesses satisfy: 0

< u 0 -Z 0 ≤ M, 0 < Z 0 -Z B ≤ M, a.e.
in Ω where M is a positive constant.

Remark 5.2.1. The assumptions (A.1)-(A.3) are classical and physically meaningful for seawater intrusion into coastal free aquifers. They are similar to the assumptions made in [START_REF] Choquet | Mathematical analysis of a sharp-diffuse interfaces model for seawater intrusion[END_REF][START_REF] Choquet | Well posedness of general cross-diffusion systems[END_REF] that dealt with the existence of a weak solution. Let us mention that with the assumption (A.3) and by introducing "sufficiently pumping" source terms, a maximum principle is shown in [START_REF] Choquet | Well posedness of general cross-diffusion systems[END_REF] that is the thicknesses u -Z and Z -Z B are bounded and positive.

Finite volume discretization

We now introduce our notation and the general form of the FV approximations that will be considered in the rest of the paper. We formulate a fully coupled fully implicit cell-centered FV scheme for the spatial discretization and a first-order (backward Euler) scheme for the integration in time. We present an MPFA [3,[START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] method which has high potential to improve efficiency and accuracy of numerical simulations for seawater intrusion in anisotropic and heterogeneous porous media. Furthermore, the orthogonality condition of the mesh is no more required which is necessary for the TPFA, see for instance [START_REF] Eymard | Convergence of finite volume schemes for semilinear convection diffusion equations[END_REF].

|K r |∇ Kr ϕ(Z δt h ) • ∇ Kr ϕ u δt h .
To estimate B 2N and B 3N , we use the same idea as done in A 2N and A 3N , to obtain

B 2N = N -1 n=0 δt K∈T h D K r=i,j,k |K r |∇ Kr ϕ Z n+1 K • ∇ Kr Z n+1 K , 1 λ N -1 n=0 δt K∈T h D K r=i,j,k |K r |∇ Kr ϕ(Z δt h ) • ∇ Kr ϕ Z δt h , and
B 3N = N -1 n=0 δt K∈T h D K δφ K r=i,j,k |K r |∇ Kr Z n+1 K • ∇ Kr Z n+1 K δ φ - λ 2 N -1 n=0 δt K∈T h D K r=i,j,k |K r |∇ Kr ϕ(Z δt h ) • ∇ Kr ϕ Z δt h .
Adding the previous inequalities and using equation (5.5.4), we get the following estimate

1 λ + δ φ - λ 2 N -1 n=0 δt K∈T h D K r=i,j,k |K r |∇ Kr ϕ(Z δt h ) • ∇ Kr ϕ Z δt h 104
Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

+ 1 λ N -1 n=0 δt K∈T h D K r=i,j,k |K r |∇ Kr ϕ(Z δt h ) • ∇ Kr ϕ u δt h 1 2 φ + Z 0 2 L 2 (Ω) .
(5.5.5)

Using the estimates (5.5.3) and (5.5.5), we obtain

N -1 n=0 δt K∈T h r=i,j,k |K r |D K ∇ Kr ϕ(u δt h ) • ∇ Kr ϕ u δt h + N -1 n=0 δt K∈T h r=i,j,k |K r |D K ∇ Kr ϕ(Z δt h ) • ∇ Kr ϕ Z δt h + N -1 n=0 δt K∈T h r=i,j,k |K r |D K ∇ Kr ϕ(Z δt h ) • ∇ Kr ϕ u δt h ≤ 1 min(δφ -/λ, 1) λφ + 2 Z 0 2 L 2 (Ω) + u 0 2 L 2 (Ω) .
Finally, using the following inequality, for any

a, b ∈ R 2 , 1 2 a 2 + b 2 ( a 2 + b 2 + a.b) in the previous inequality (with a = ∇ Kr ϕ(Z δt h ), b = ∇ Kr ϕ(u δt h )), we obtain N -1 n=0 δt K∈T h r=i,j,k |K r |D K ∇ Kr ϕ(u δt h ) • ∇ Kr ϕ u δt h + N -1 n=0 δt K∈T h r=i,j,k |K r |D K ∇ Kr ϕ(Z δt h ) • ∇ Kr ϕ Z δt h ≤ C where C = λφ + min(δφ -/λ, 1) Z 0 2 L 2 (Ω) + u 0 2 L 2 (Ω)
. This concludes the proof of Proposition 5.5.1

Existence of the discrete solution

In this subsection, we will prove an existence result for the nonlinear system of equations (5.3.19)- (5.3.25). To this end, we follow the strategy proposed in [START_REF] Cao | Convergence of an mpfa finite volume scheme for a two-phase porous media flow model with dynamic capillarity[END_REF][START_REF] Michel | A finite volume scheme for two-phase immiscible flow in porous media[END_REF] and we use Leray-Schauder's fixed point theorem.

U = u n+1 K , u n+1 σ 1 Kr , u n+1 σ 2 Kr K∈T h ,r∈{i,j,k},n∈[0,N -1] and Z = Z n+1 K , Z n+1 σ 1 Kr , Z n+1 σ 2 Kr K∈T h ,r∈{i,j,k},n∈[0,N -1]
.

Proof. Let (U, Z) be a solution of the system (5.3.19)- (5.3.25). Let E = R [0,N -1]×T h 2 and we define the mapping

G : E → E (U, Z) → Ũ, Z ,
such that Ũ, Z is the solution of the following set of equations: for all K ∈ T h , Lemma 5.5.2. Under our standing assumptions, for every s ∈ R + , the following uniform in s holds true:

T -s 0 Ω ϕ(u δt h (t + s, x)) -ϕ(u δt h (t, x)) 2 dx dt ≤ Cs (5.5.6)
where C is a constant independent of the discretization parameters and s.

Proof. Let s ∈]0, T [. Since ϕ is a non-decreasing, Lipschitz function with λ its Lipschitz constant, the following inequality holds

T -s 0 Ω ϕ(u δt h (t + s, x)) -ϕ(u δt h (t, x)) 2 dx dt ≤ λ T -s 0 A(t) dt, with
A(t) = Ω ϕ(u δt h (t + s, x)) -ϕ(u δt h (t, x)) u δt h (t + s, x) -u δt h (t, x) dx, ∀ t ∈]0, T -s[. Let n 0 (t) ∈ {0, ..., N -1} such that t n 0 (t) ≤ t ≤ t n 0 (t) + 1 and n 1 (t) ∈ {0, ..., N -1} such that t n 1 (t) ≤ t + s ≤ t n 1 (t) + 1.
The following equality yields

A(t) = K∈T h ϕ u n 1 (t) K -ϕ u n 0 (t) K × N -1 n=0 χ n (t, t + s) u n+1 K -u n K |K|, where χ n (t, t + s) = 1, if t n+1 belongs to [t, t + s], 0, otherwise.
Using the definition of the FV scheme (5.3.19)-(5.3.25) and changing the order of the summations between n and K, one gets

A(t) = δt N -1 n=0 χ n (t, t + s) K∈T h ϕ u n 1 (t) K -ϕ u n 0 (t) K × r=i,j,k Id=1,2 |σ Id Kr | D K φ K ∇ Kr ϕ u δt h + ∇ Kr ϕ u δt h + δφ K ∇ Kr u δt h • n σ Id Kr .
This can be written as follows

A(t) = A 1 (n 1 (t)) -A 1 (n 0 (t)) + A 2 (n 1 (t)) -A 2 (n 0 (t)) + A 3 (n 1 (t)) -A 3 (n 0 (t)). Let ρ = n 1 (t) or ρ = n 0 (t). One has A 1 (ρ) = δt N -1 n=0 χ n (t, t + s) K∈T h D K φ K r=i,j,k Id=1,2 |σ Id Kr | ∇ Kr ϕ u δt h • n σ Id Kr ϕ u ρ K , A 2 (ρ) = δt N -1 n=0 χ n (t, t + s) K∈T h D K φ K r=i,j,k Id=1,2 |σ Id Kr | ∇ Kr ϕ Z δt h • n σ Id Kr ϕ u ρ K , A 3 (ρ) = δt N -1 n=0 χ n (t, t + s) K∈T h δD K r=i,j,k Id=1,2 |σ Id Kr | ∇ Kr u n+1 K • n σ Id Kr ϕ u ρ K .
To facilitate the computations, we rewrite the above inequality as follows

A(t) ≤ D + φ -(2 + δφ -)(E 0 + E 1 ) + (1 + C)E 2 + E 3 , with E 0 = δt N -1 n=0 χ n (t, t + s) K∈T h r=i,j,k Id=1,2 |K r | ϕ u n 0 (t) σ Id Kr -ϕ u n 0 (t) K 2 | µ σ Id Kr | 2 , E 1 = δt N -1 n=0 χ n (t, t + s) K∈T h r=i,j,k Id=1,2 |K r | ϕ u n 1 (t) σ Id Kr -ϕ u n 1 (t) K 2 | µ σ Id Kr | 2 , E 2 = δt N -1 n=0 χ n (t, t + s) K∈T h r=i,j,k Id=1,2 |K r | ϕ u n+1 σ Id Kr -ϕ u n+1 K 2 | µ σ Id Kr | 2 , E 3 = δt N -1 n=0 χ n (t, t + s) K∈T h r=i,j,k Id=1,2 |K r | ϕ Z n+1 σ Id Kr -ϕ Z n+1 K 2 | µ σ Id Kr | 2 . We have t ∈]0, T -s[, n 0 (t) ∈ {0, ..., N -1}, m = n 0 (t) if and only if t ∈ [t m , t m+1 ]. As a consequence T -s 0 E 0 dt ≤ N -1 m=0 t m+1 tm δt   N -1 n=0 χ n (t, t + s) K∈T h r=i,j,k Id=1,2 |K r | ϕ u m σ Id Kr -ϕ (u m K ) 2 | µ σ Id Kr | 2   dt, ≤ N -1 m=0 δt t m+1 tm N -1 n=0 χ n (t, t + s) dt K∈T h r=i,j,k Id=1,2 |K r | ϕ u m σ Id Kr -ϕ (u m K ) 2 | µ σ Id Kr | 2 . Since N -1 m=0 t m+1 tm N -1 n=0 χ n (t, t + s) dt ≤ s, the following inequalities hold T -s 0 E 0 dt ≤ s ϕ u δt h 2 1,T h ,δt and
T -s 0 E 1 dt ≤ s ϕ u δt h 2 1,T h ,δt .
In the same way from

T -s 0 χ n (t, t + s) dt ≤ s, we obtain T -s 0 E 2 dt ≤ s ϕ u δt h 2 1,T h ,δt , and
T -s 0 E 3 dt ≤ s ϕ u δt h 2 1,T h ,δt .
Adding the previous inequalities and using Proposition 5.5.1 lead to the desired result

T -s 0 Ω ϕ(u δt h (t + s, x)) -ϕ(u δt h (t, x)) 2 dx dt ≤ Cs,
where C is a constant independent of the discretization parameters and s. This completes the proof of Lemma 5.5.2.

Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

Convergence results

In this section, we present some convergence results based on the estimates established in the last sections and Kolmogorov's theorem and we show that the limit of a sub-sequence of the approximate solution (5. (5.3.25). Then, we can extract a sub-sequence which converges in an adequate sense to a weak solution (u, Z), when h m and δt m tends to zero as m → ∞, and satisfying the following results:

1. Convergence results ϕ u δtm hm → ϕ (u) strongly in L 2 (0, T ; Ω) , u δtm hm → u a.e. in Ω T , ϕ Z δt m hm → ϕ (Z) strongly in L 2 (0, T ; Ω) , Z δtm hm → Z a.e. in Ω T , ∇ D ϕ u δtm hm ∇ϕ (u) weakly in L 2 (0, T ; Ω) 2 , ∇ D ϕ Z δtm hm ∇ϕ (Z) weakly in L 2 (0, T ; Ω) 2 .
2. The limit (u, Z) is a weak solution to the continuous problem (P) in the following sens

T 0 Ω φ ∂ψ ∂t u dx dt - T 0 Ω D ∇ϕ(u) + ∇ϕ(Z) + δφ∇u • ∇ψ dx dt + Ω u 0 (x)ψ(0, x) dx = 0, T 0 Ω φ ∂ζ ∂t Z dx dt - T 0 Ω D b s (Z)∇u + ∇ϕ(Z) + δφ∇Z • ∇ζ dx dt + Ω Z 0 (x)ζ(0, x) dx = 0, ∀(ψ, ζ) ∈ C ∞ 0 (0, T ; Ω) 2 .
Proof of Theorem 5.6.1

Step 1. Convergence results. From the space and time translations Lemmas 5.5.1 and 5.5.2 and due to Kolmogorov's theorem, the sequence ϕ u δtm hm m∈N * is relatively compact in L 2 (0, T ; Ω). Therefore, there exists a sub-sequence, still denoted by ϕ u δtm hm m∈N * , which converges strongly to ϕ in L 2 (0, T ; Ω), as m → ∞. The function ϕ being a strictly non-decreasing, so there is a continuous reciprocal function ϕ -1 such that u δtm hm → ϕ -1 ϕ a.e. in Ω T .

Thanks to the L ∞ bound given in Proposition 5.4.1, we deduce from Lebesgue's dominated convergence theorem that u δtm hm → u = ϕ -1 ϕ strongly in L 2 (0, T ; Ω) as h m , δt m → 0.

We follow the same method as done for ϕ u δtm hm . We can prove the space and time translation estimates for ϕ Z δtm hm . Then, we apply the same arguments to obtain the compactness of ϕ Z δtm hm using Kolmogorov's theorm. This yields to the existence of a sub-sequence, still denoted by ϕ Z δtm hm , wich converges, as m → ∞, strongly in L 2 (0, T ; Ω). Using again Lebesgue's dominated convergence theorem, we deduce that Z δtm hm → Z strongly in L 2 (0, T ; Ω) as h m , δt m → 0 when m → ∞. Now, we demonstrate a weak convergence of the discrete gradient. Proposition 5.5.1 gives an estimate of the discrete gradient in the space L 2 (0, T ; Ω)

2 , one has

Nm-1 n=0 δt m ∇ D ϕ u δtm hm 2 2 + Nm-1 n=0 δt m ∇ D ϕ Z δtm hm 2 2 ≤ C.
From this estimate, one can deduce that

∇ D ϕ u δtm hm resp. ∇ D ϕ Z δtm hm is bounded uniformly in L 2 (0, T ; Ω) 2 .
Then, there exists a sub-sequence such that

∇ D ϕ u δtm hm ζ 1 and ∇ D ϕ Z δtm hm ζ 2 weakly in L 2 (0, T ; Ω) 2 as m → ∞.
Using the strong convergence of ϕ u δtm hm and the weak convergence of the discrete gradient ∇ D ϕ u δtm hm , we obtain by passing to the limit

T 0 Ω ∇ D ϕ u δtm hm • V + ϕ u δtm hm • divV d xdt → T 0 Ω ζ 1 • V + ϕ (u) • divV d xdt, as m → ∞ for any V ∈ C ∞ 0 (Ω × (0, T ))
2 . This prove by the way that ∇ϕ(u) = ζ 1 (resp.

∇ϕ(Z) = ζ 2 ) in the sense of distribution and ϕ (u) , ϕ (Z) ∈ L 2 0, T ; H 1 (Ω) .

Step 2. Passage to the limit. Let ψ ∈ A T = {ψ ∈ C ∞ 0 (0, T ; Ω) | ∇ψ • n = 0 on ∂Ω T and ψ(., T ) = 0}. We denote

ψ n+1 K = ψ(x K , t n+1 ) and ψ n+1 σ Id Kr = ψ(x σ Id Kr , t n+1 )) for Id = 1, 2.
We multiply equation (5.3.19)

T 1m + T 2m + T 3m + T 4m = 0, with T 1m = Nm-1 n=0 K∈T hm |K|φ K u n+1 K -u n K ψ n+1 K , T 2m = Nm-1 n=0 δt m K∈T hm r=i,j,k |K r |D K ∇ Kr ϕ u δtm hm • ∇ Kr ψ δtm hm , T 3m = Nm-1 n=0 δt m K∈T hm r=i,j,k |K r |D K ∇ Kr ϕ Z δtm hm • ∇ Kr ψ δtm hm , T 4m = Nm-1 n=0 δt m K∈T hm r=i,j,k |K r |δφ K D K ∇ Kr u δtm hm • ∇ Kr ψ δtm hm .

scheme for a sharp-diffuse interfaces model

We rearrange the summation of T 1m and by using ψ Nm K = ψ(x K , T ) = 0, we get

T 1m = - Nm-1 n=0 K∈T hm |K|φ K u n+1 K ψ n+1 K -ψ n K + K∈T hm |K|φ K u Nm K ψ Nm K -u 0 K ψ 0 K , = - Nm-1 n=0 K∈T hm t n+1 tn K u n+1 K ∂ψ(x K , t) ∂t dx dt - K∈T hm |K|φ K u 0 K ψ 0 K . Now, by the strong convergence of ψ(x K , 0) K∈T hm , one has lim m→∞ K∈T hm |K|φ K u 0 K ψ 0 K = Ω φ(x)u 0 (x)ψ(x, 0) dx.
By the regularity of the function ψ, the ∂ t ψ(x K , .) K∈T hm converges uniformly towards ∂ t ψ(x, .). We deduce from Lebesgue's dominated convergence theorem that

lim m→∞ T 1m = - Ω φ(x)u 0 (x)ψ(x, 0) dx - T 0 Ω φ(x)u(t, x)∂ t ψ(t, x) dxdt.
We next prove the following limit

E 1m + E 2m + E 3m + E 4m = 0,

Convergence results

113

with

E 1m = Nm-1 n=0 K∈T |K|φ K Z n+1 K -Z n K ζ n+1 K , E 2m = Nm-1 n=0 δt m K∈T hm r=i,j,k |K r |D K ∇ Kr b s (Z δtm hm ) up σ Id Kr ∇ Kr u δtm hm • ∇ Kr ζ δtm hm , E 3m = Nm-1 n=0 δt m K∈T hm r=i,j,k |K r |D K ∇ Kr ϕ Z δtm hm • ∇ Kr ζ δtm hm , E 4m = Nm-1 n=0 δt m K∈T hm r=i,j,k |K r |δφ K D K ∇ Kr Z δtm hm • ∇ Kr ζ δtm hm .
As done previously, we have the following limits

lim m→∞ E 1m = - Ω φ(x)Z 0 (x)ζ(x, 0) dx - T 0 Ω φ(x)Z(t, x) ∂ζ(t, x) ∂t dxdt, lim m→∞ E 3m = T 0 Ω D∇ϕ (Z) • ∇ψ dx dt and lim m→∞ E 4m T 0 Ω δφD∇u • ∇ψ dx dt.
We next determinate the limit of E 2m , To do so, we consider the residual function R 2m defined by

R 2m = T 0 Ω Db s (Z δtm hm )∇ D u δtm hm • ∇ D ζ δtm T hm dx dt - T 0 Ω Db s (Z)∇u • ∇ζ dx dt, = R 21m + R 22m , with R 21m = T 0 Ω D b s (Z δtm hm )∇ D u δtm hm • ∇ D ζ δtm T hm -b s (Z)∇ D u δtm hm • ∇ D ζ T hm ,δtm dx dt, R 22m = T 0 Ω b s (Z)D ∇ D u δtm hm • ∇ D ζ δtm T hm -∇u • ∇ζ dx dt.
The goal here is to show that R 2m tends to zero when passing to the limit. Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

Numerical results

In this section, we present the numerical results for two test cases modeling different scenarios of the sharp-diffuse interfaces model for seawater intrusion. The first test case is a simulation in a homogeneous domain, while the second one is in a heterogenous free surface aquifer where two intersecting fractures are considered.

All our developments have been implemented in DuMu X . It provides many tools to solve numerically PDEs 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 schemes.

In this context, we have developed and implemented in DuMu X a new module, named 2p -SW I, which allows to solve the coupled system (5.3.19)- (5.3.25). The primary variables are the depth of the saltwater/freshwater interface Z and the hydraulic head of freshwater u. The spatial discretization employs a TPFA or MPFA finite volume method. The time discretization is done by an implicit Euler scheme. The nonlinear system is solved by the Newton method and a BiConjugate Gradient STABilized (BiCGSTAB) method preconditioned by an Algebraic Multigrid (AMG) solver, is used to solve the linear systems. This solver is integrated in the ISTL-Library of DUNE. Dynamic time stepping is used when integrating the equations involved. Time step sizes during transient simulations are dynamically recalculated depending on the convergence behavior of the nonlinear procedure (Newton method). The simulations start with a certain time step which can be increased or reduced, depending on the number of iterations allowed in each nonlinear iteration. Let us mention that throughout all numerical experiments, we observed that in no instance more than a maximum of 10 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. Therefore, there is no need to use other types of local time-stepping strategies proposed in the literature.

All computations were performed on a Del PC with Intel Core(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.

Our approach has been validated by solving several tests, including numerical convergence of the schemes [4]. The numerical results are satisfactory and replicated to those in the literature. The results of these simulations are omitted since nothing startling was found. Instead, we concentrate on the results obtained in two test cases. Finally a remarkable property of the scheme is that the discrete maximum principle (nonnegativity of the thickness of fresh and salt water in the aquifer) is satisfied wich is crucial to obtain physically meaningful approximate solutions. This has been verified in all our simulations.

Studied domain and physical data

We consider two test cases by adapting the geometry and data of an example proposed in [START_REF] Choquet | Derivation of a sharp-diffuse interfaces model for seawater intrusion in a free aquifer. numerical simulations[END_REF] dealing with the evolution of transition zones under significant solicitations The example presented above serves to visualize the numerical behavior of the transition zone considering two situations. In the first situation, the aquifer is supposed to be homogeneous, whereas it is taken to be heterogeneous in the second situation. In both situations, two pumping wells, with different flow rates, are also involved. The freshwater extraction rate in the pumping wells is assumed to be constant per day: When dealing with the homogeneous case, the horizontal surface is discretized by a structured mesh with the orthogonality condition being satisfied. The numerical scheme is established using a TPFA method. In the heterogeneous case, the current problem is coupled with a fracture possessing a high conductivity. The domain in this case is discretized using an unstructured triangulation mesh. To deal numerically with such a system, we use the MPFA method. The impact of heterogeneity on the displacement of transition zones is illustrated by numerical investigations. In the following, we present the obtained numerical simulations of each case.

Q f = Q 1 χ ]-
We performed the simulations with an initial time step of 10 seconds and a maximum time step of 1 hour. 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 4 iterations.

Test case 1: homogeneous aquifer

In this case, we focus on the evolution of the transition zones in an homogeneous free aquifer subjected to over-exploitation. The hydraulic conductivity in this case is equal to D = 4.5 10 -4 m/s. Figure 5.2 shows the structured mesh of 21164 cells used in the numerical simulations. To visualize the extent and advancement of the saltwater intrusion as a function of the location, vertical cross-section through the two pumping wells are performed. Figure 5.3 shows the position of the freshwater/saltwater transition zone at the vertical cross-section 1 through the two pumping wells, in the free aquifer subjected to pumping during 3 and 10 days. freshwater/saltwater interface is shown in Figure 5.5. It can be observed that "upconings" create beneath the pumping wells and grow up quickly. A local depression, located at the main zone of pumping, is also developed. The presence of the freshwater/saltwater transition zone and the saturated/unsaturated zone promotes diffusion and advancement of the free surface and salt front as it can be seen in Figure 5.6.

Test case 2: fractured aquifer

In this case, the impact of heterogeniety on the transition zone is examined. The aquifer described previously is no more homogeneous but crossed by a network of two fractures. The main one is situated in the middle of the aquifer in the length direction, whereas the To monitor the location and the extent of the transition zone, three vertical cross-sections are chosen. The first one passes through the two wells (cross-section 1), the second, instead, is lying along the main fracture (y=0, cross-section 2) and the third one is lying in the direction of the inclined fracture (cross-section 3). of the transition zone due to the greatest amount constitutes a threat for the quality of freshwater extrude in the pumping wells. The existence of the fracture accelerate the freshwater and saltwater flows and contributes to the advancement of the transition zone (Figures 5.9-5.11). After 10 days of over-exploitation, the salt bevel converges to the equilibrium state.

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

Conclusion

In this chapter, we have studied the convergence of an MPFA finite volume method for a sharp-diffuse interfaces model in an heterogeneous free aquifer. A fully implicit, fully coupled finite volume method has been developed to discretize the governing equations. This approach deals well with the anisotropy and heterogeneity of the soil and provides accurate and efficient approximations on an unstructured mesh. The construction of the numerical scheme is established by integrating the governing equations on each control volume and evaluating the fluxes on each half edge. To linearise the coupled system at each time step, we employ a Newton's method wherein each iteration a linear system is solved by an iterative Krylov method.

After that, we have shown that the numerical scheme preserves the positivity of the freshwater and the saltwater thicknesses. Moreover, we have proved an existence result 120 Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model to the nonlinear fully implicit numerical scheme using a fixed point theorem. The strong convergence of the approximate solution to the weak solution of the continuous problem is established via some recent compactness tools corresponding to the space and time translate estimates. In addition, we have proved the weak convergence of the discrete gradients of the freshwater head and salt front elevation.

We have developed and implemented this scheme in the DuMu X framework. To ensure the validity of our method, we have considered two test cases. The first test case is a simulation in a homogeneous domain, while the second one is in a heterogenous free surface aquifer where two intersecting fractures are considered. The obtained results proved that our new module based on recent numerical tools is accurate, efficient and able to solve a 2D seawater intrusion model tacking into account the dynamics of the transition zones.

Chapter 6

Conclusion and perspectives

In this chapter, we summarize the main contributions and achievements made in this work in the context of simulating seawater intrusion in coastal aquifers. We have presented accurate and efficient algorithms to solve a seawater intrusion problem using different approaches. A brief presentation of the state of the art of existing approaches is proposed. First of all, we have presented the derivation of the mathematical model using the sharp interface model in a free and confined aquifer. To numerically deal with such a system, we have considered a fully coupled, fully implicit finite volume TPFA method on a triangulation mesh where the orthogonality condition is satisfied. The nonlinear system is linearized using Newton's method, in which a linear system is solved at each iteration using iterative Krylov solvers, provided by the external DUNE library. To predict the extent of salt front elevation and water table level in large-scale aquifers, we have integrated the numerical scheme in the framework DuMu X . The implemented module is based on recent and efficient numerical tools allowing accurate predictions. Different 1D, 2D test cases are used for the validations of our module, including real tests. The obtained results correspond well to those presented in [START_REF] Abudawia | Numerical analysis for a seawater intrusion problem in a confined aquifer[END_REF][START_REF] Marion | Numerical simulations for a seawater intrusion problem in a free aquifer[END_REF][START_REF] Mehdizadeh | Vertical leakage in sharp-interface seawater intrusion models of layered coastal aquifers[END_REF]. The numerical module is then applied to the Souss Chtouka plain and the Tripoli aquifer. Numerical solutions based on available hydrogeological data proved the ability of our methodology to predict the progress and the extent of the salt-wadge for several years of over-exploitation of the nappe.

The sharp interface approach assumes that the two fluids are separated by an abrupt interface. This encouraged us to adopt a mixed sharp-diffuse interfaces approach, introduced recently in [START_REF] Choquet | Mathematical analysis of a sharp-diffuse interfaces model for seawater intrusion[END_REF] that takes into account the dynamics of transition zones. In this context, we have developed a fully implicit, fully coupled finite volume method for the numerical resolution of such system. One again, the numerical scheme is implemented in the DuMu X framework. The developed module is validated on test cases dealing with the progression of the salt front coupled with sea fluctuations. After that, we have compared the numerical result of the 2D model with those obtained with the classical 3D model for miscible displacements. The numerical results are satisfactory and showed good agreement with those presented in [START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF].

At the real scale, the aquifer is heterogeneous and anisotropic, which is why it is necessary to use a more suitable numerical method on unstructured meshes. In this work, we have therefore developed a fully coupled, fully implicit MPFA finite volume method for the numerical approximation of the sharp-diffuse model. The numerical analysis of this scheme is presented. We have proved the nonnegativity of the approximate freshwater and

3. 2

 2 Discretisation by a cell-centred finite volume method. . . . . . . . . . . . 3.3 Position of the interface obtained with the present model and the analytical solution at different simulation times. . . . . . . . . . . . . . . . . . . . . 3.4 The progress of the numerical and analytical toe. . . . . . . . . . . . . . . 3.5 Evolution of the salt interface at different times. . . . . . . . . . . . . . . 3.6 Rotary character of the freshwater/saltwater interface at times t = 10, 20, 100, 400 and 800 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Evolution of the groundwater surface Z T and the salt front Z in an injection (left) and pumping (right) process at time: t=1.5 and 10 days. . . . . . . 3.8 Left: Injection process. Right: pumping scenario. . . . . . . . . . . . . . 3.9 Schematic of the well screen in an unconfined aquifer [85]. . . . . . . . . . 3.10 Convergence of solutions on different grids to the reference solution during 30 years of activity on the main section (0, 10) -(500, 10). . . . . . . . . 3.11 Evolution of the solutions on each grid and convergence to the reference solution for 30 years at the point of pumping (350, 10). . . . . . . . . . . 3.12 Variation of the salt front while changing the pumping rate and fixing the position of the well. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Evolution of the salt front by varying the depth of the well and fixing the pumping rate at Q t = 0.1 m 3 /day. . . . . . . . . . . . . . . . . . . . . . . 3.14 Evolution of the salt front by varying the depth of the well and fixing the pumping rate at Q t = 0.07 m 3 /day. . . . . . . . . . . . . . . . . . . . . . . 3.15 Evolution of the salt front for Sc-1, Sc-7, Sc-8 and Sc-9 scenarios, while changing the longitudinal positions of the wells. . . . . . . . . . . . . . . .

3. 16

 16 The situation of the Souss-Massa basin and the Souss-Chtouka groundwater. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.17 3D representation of the studied domain. . . . . . . . . . . . . . . . . . . 47 3.18 Left: geometry and boundaries of the aquifer and well locations. Right: mesh for the aquifer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.19 Left: hydraulic conductivity (D f ) of the aquifer. Right: porosity (φ) of the soil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.20 Left: the freshwater storage coefficient S f . Right: comparison between our developed module 2p-SWI model and the measured piezometric map of 1968. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.21 Topography of the aquifer bottom. . . . . . . . . . . . . . . . . . . . . . . 50 3.22 Plan view of the piezometry of the plain Souss-Chtouka before (left) and after (right) solicitation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.23 Contours representing the freshwater potential of the Souss-Chtouka plain illustrated on the bottom of the aquifer in 1968 (left) and 2048 (right). . . 53 3.24 Plan view of the piezometry of the saltwater in its initial state in 1968 (left) and after 80 years of exploitation (right). . . . . . . . . . . . . . . . . . . 53 3.25 Position of the interface freshwater/saltwater in 1968 (black line) and in 2048 (red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.26 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 time: cut on the level 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.27 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 time: cut on the levels 2 (left) and 3 (right). 55 3.28 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). 56 3.29 Evolution in time of the salt front (top right) freshwater head (left) and saltwater head (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.30 Position of the Tripoli aquifer, study zone and wells position (cf. [57]). . 58 3.31 Geological map and stratigraphy of Tripoli (cf. [57]). . . . . . . . . . . . 59 3.32 Geometry and Boundaries of the aquifer (left). The mesh of the studied domain (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.33 Location of pumping zones (cf. [57]). The rectangle in blue corresponds to the private wells and the public wells are represented by a yellow zone. The yellow zone in the south of the aquifer contain 4 wells in Bhsas while the two other wells are in El-jisr and Malouleh respectively. . . . . . . . 61 3.34 Evolution of the saltwater/freshwater interface in 2008. "NOP" scenario . 62

3. 38

 38 Evolution of the piezometric head of freshwater in 2008. . . . . . . . . . . 3.39 Position of the freshwater head scenarios. a) After 10 years of pumping with 150 L/day: "REC" scenario. b) After 10 years of pumping with 250 L/day: "10Y" scenario. c) After 20 years of pumping with 250 L/day: "20Y" scenario. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Schematic representation of a coastal aquifer. . . . . . . . . . . . . . . . . 4.2 Convergence of solutions on each grid to the reference solution during 3 days of pumping on cross-section (-50, 0) -(50, 0). . . . . . . . . . . . . 4.3 Evolution of the solutions on each grid and convergence to the reference solution for 3 days at the point (-10, 0). . . . . . . . . . . . . . . . . . . 4.4 Conservation of the total amount of freshwater pumped . . . . . . . . . . 4.5 Evolution of the elevation of sharp-diffuse interfaces (left) and the freshwater head (right) in 1, 2 and 3 days. . . . . . . . . . . . . . . . . . . . . . . 4.6 Velocity of freshwater (left) and saltwater (right) flows during 3 days of pumping. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Convergence of the interface obtained with the present model towards the analytical solution derived from Ferris model. Times T = 0.5 day (left) and T = 1 day (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Evolution of the freshwater head and the salt front elevation in differnt times without any external forces. T= 0.1, 0.5 and 1 day. . . . . . . . . . 4.9 Displacement of the freshwater head and the salt front elevation in different times during a pumping scenario. T=0.1, 0.5, 1 day. . . . . . . . . . . . . 4.10 Oscillations on the salt front elevation, predicted by the hydrograph artifacts (δ = 0 (left) and δ = 0(right)) and comparison between both models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Sea oscillations effects on the freshwater head u. Simulations with δ = 0 and δ = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Evolution of the freshwater head u with and without oscillations in the boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Initial conditions : initial salt concentration for the 3D model (left) and initial position of the transition zone between freshwater/saltwater for the 2D model (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 3D presentation of the initial position of the interface. . . . . . . . . . . . 4.15 Evolution of the freshwater and saltwater flows in the normal conditions during T = 10 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4. 16 92 5. 2 1 1. 2 2 1. 3 3 1. 3 . 1 4 1. 3 . 2 4 1. 4

 16922122333143244 Representation of sharp-diffuse interfaces and 3D density variable results (50 % salinity contour and 5 % and 95% salinity contour) for T=5 days. .[START_REF] Sebben | Seawater intrusion in fractured coastal aquifers: A preliminary numerical investigation using a fractured henry problem[END_REF] 4.17 The velocity field in the rotating interface problem at t = 10 days. . . . . 83 4.18 Numerical results for pumping scenario: x 2 = 0.5 (top) x 2 = 1.5 (center) and x 2 = 2.5 (below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.19 Evolution of the concentration (3D model) and the salt front (2D model) under a large amount of pumping: slice on x 2 = 1.5. . . . . . . . . . . . . 84 4.20 Numerical results during an injection experiment. x 2 = 0.5 (top), x 2 = 1.5 m (center) and x 2 = 0.5 (bottom) . . . . . . . . . . . . . . . . . . . . . . . 85 4.21 Initial conditions: initial salt concentration for 3D model (left) and initial position of the salt front for the 2D model (right). . . . . . . . . . . . . . 86 4.22 Numerical results for double pumping wells. Slice x 2 = 0.5 m (top), slice x 2 = 1.5 m (center) and slice x 2 = 2.5 m. . . . . . . . . . . . . . . . . . . 86 4.23 Numerical results: injection and pumping, T = 10.5 days γ = 0.00025. . . 87 5.1 The control volume K and associated notations. . . . . . . . . . . . . . . Structured mesh of a homogeneous aquifer (21164 cells). . . . . . . . . . . 116 5.3 Displacement of the freshwater/saltwater transition zone at vertical crosssection through pumping wells in the homogeneous case after 3 days (left) and 10 days (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.4 Freshwater flow velocity in the homogeneous case after 10 days. . . . . . . 116 5.5 Displacement of the freshwater/saltwater interface 2D plot in the homogeneous case after 3 days (left) and 10 days (right). . . . . . . . . . . . . . . 117 5.6 Evolution of the freshwater head 2D plot in the homogeneous case after 3 days (left) and 10 days (right). . . . . . . . . . . . . . . . . . . . . . . . 117 5.7 Unstructured mesh in the heterogenous case with two intersecting fractures (20651 cells). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.8 Displacement of the freshwater/saltwater transition zone for different sections in the heterogeneous case after 3 days (left) and 10 days (right). . . 118 5.9 Freshwater flow velocity in the heterogeneous case after 10 days. . . . . . 119 5.10 Displacement of the freshwater/saltwater interface 2D plot in the heterogeneous case after 3 days (left) and 10 days (right). . . . . . . . . . . . . 119 5.11 Evolution of the freshwater head 2D plot in the heterogeneous case after 3 days (left) and 10 days (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Motivation and problem setting . . . . . . . . . . . . . . . . . . . . . . . . . The objective of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . The sharp interface module . . . . . . . . . . . . . . . . . . . . . . . The sharp-diffuse interfaces module . . . . . . . . . . . . . . . . . . Thesis structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

Figure 1 .

 1 1 represents seawater intrusion phenomena into coastal aquifers.

Figure 1 . 1 :

 11 Figure 1.1: Saltwater intrusion into coastal water wells. Source: http://imag.slid.com

  Time variable [s], x: Space variable [m], i = f, s: Phase index f= freshwater and s= saltwater. The main unknowns of the system u: The freshwater head [m], v: The saltwater head [m], -Z: The salt front elevation [m]. Physical parameters φ: porosity of medium [-], κ: Intrinsic permeability tensor m 2 , µ i : Dynamic viscosity kg.m -1 .s -1 , ρ i : Density of the fluid i kg • m -3 , q i : Darcy's velocity of the fluid i m.s -1 , -Q i : Quantity of water injected or pumped per unit of VER s -1 , g: Gravitational acceleration m.s -2 , -D : Hydraulic conductivity of freshwater [m/day] -D i : Hydraulic conductivity of the fluid i [m/day]

Figure 2 . 1 :

 21 Figure 2.1: Sharp interface model. Source: https://www.solinst.com

Figure 2 . 2 :

 22 Figure 2.2: Groundwater flow patterns and the dispersion zone in a homogeneous coastal aquifer. Source: http://water.usgs.gov/

where κ: intrinsic permeability tensor m 2

 2 , µ: dynamic viscosity kg.m -1 .s -1 , g: gravitational acceleration m.s -2 , x 3 : the vertical dimension [m] .

Figure 4 . 1 )

 41 . The function Z B [m] describes its lower topography while Z T [m] represents the elevation of the top of the aquifer.

Figure 2 . 3 :

 23 Figure 2.3: Schematic representation of a coastal aquifer [71].

Contents 3 . 1

 31 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 A fully implicit finite volume numerical scheme . . . . . . . . . . . . . . . . 24 3.2.1 Discretization and basic notations . . . . . . . . . . . . . . . . . . . 25 3.2.2 Resolution of the nonlinear system . . . . . . . . . . . . . . . . . . . 28 3.3 Implementation of the numerical scheme . . . . . . . . . . . . . . . . . . . . 31 3.4 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4.1 Test 1: The rotating interface problem . . . . . . . . . . . . . . . . . 32 3.4.2 Test 2: a field-scale free aquifer problem . . . . . . . . . . . . . . . 36 3.4.3 Test 3: Souss-Chtouka case study . . . . . . . . . . . . . . . . . . . . 46 3.4.4 Test 4: Tripoli aquifer case study . . . . . . . . . . . . . . . . . . . . 58 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

Figure 3 . 1 :

 31 Figure 3.1: Geometry of the domain Ω.

Figure 3 . 2 :

 32 Figure 3.2: Discretisation by a cell-centred finite volume method.

  2.11) The derivatives of the local freshwater residue R f u n+1 k , v n+1 k on the element V k and at time t n+1 with respect to u n+1 k and v n+1 k are the following:

1 :

 1 Parameter values for the rotating interface problem.Simulations were conducted over a 20-day period from the initial position. The simulations start run with a resolution of 100 × 100. Figure3.3 shows the position of the interface obtained with this model (curves "Z-Ap1t') and the analytical solution (curves "Zexact') at time t = 1, 10 and 20 days.

Figure 3 . 3 :

 33 Figure 3.3: Position of the interface obtained with the present model and the analytical solution at different simulation times.

Figure 3 . 4 :

 34 Figure 3.4: The progress of the numerical and analytical toe.

  Initial position t = 0.0. t = 10.0 days. t = 20.0 days. t = 100 days. t = 200 days. t = 400 days.

Figure 3 . 5 :

 35 Figure 3.5: Evolution of the salt interface at different times.

Figure 3 . 6 :

 36 Figure 3.6: Rotary character of the freshwater/saltwater interface at times t = 10, 20, 100, 400 and 800 days.

Figure 3 . 7 :

 37 Figure 3.7: Evolution of the groundwater surface Z T and the salt front Z in an injection (left) and pumping (right) process at time: t=1.5 and 10 days.

Figure 3 .

 3 Figure 3.7 shows the movement of the interfaces Z T and Z during a pumping and injection processes. We observe a local depression of the water table and the appearance of a salt cone " upconing" in the pumping area. The salt interface rises globally towards the horizon, which can be dangerous for the users. As to the injection scenario, the effects of freshwater filling with a decrease in the depth of salt water are observed, especially in the injection zone.

Figure 3 . 8 :

 38 Figure 3.8: Left: Injection process. Right: pumping scenario.

Figure 3 .

 3 Figure3.8 shows a comparison of the interfaces Z T and Z in the case of a compressible soil (S f = 0.25, the Z-InjSft curves) and in the case of an incompressible environment (S f = 0, the Z-Injt curves) at time t = 1, 2, 3 days. We can see an overestimation of the groundwater interface level Z T and an underestimation of the salt front interface elevation Z when S f is zero. Usually, the injection process is used to decrease the salt cone and limit the marine intrusion problem. However, we remark the opposite phenomenon during a pumping process. Therefore, we see here that neglecting the compressibility of the soil can produce estimation errors on the evolution of the Z T and Z interfaces, particularly in the injection area.

4 :

 4 Grids considered in the study of the numerical convergence.

Figure 3 . 10 :

 310 Figure 3.10: Convergence of solutions on different grids to the reference solution during 30 years of activity on the main section (0, 10) -(500, 10).

Figure 3 . 11 :

 311 Figure 3.11: Evolution of the solutions on each grid and convergence to the reference solution for 30 years at the point of pumping (350, 10).

Figure 3 . 12 :

 312 Figure 3.12: Variation of the salt front while changing the pumping rate and fixing the position of the well.

Figure 3 . 13 :

 313 Figure 3.13: Evolution of the salt front by varying the depth of the well and fixing the pumping rate at Q t = 0.1 m 3 /day.

Figure 3 . 14 :

 314 Figure 3.14: Evolution of the salt front by varying the depth of the well and fixing the pumping rate at Q t = 0.07 m 3 /day.

  c) Case 3: varying the longitudinal position of the wellFour scenarios are considered to study the effect of the distance of the well from the coast on the evolution of the salt front. To do so, we fixed the amount of freshwater at Q t = 0.1 m 3 /day and we varied the distance from the coast x w . The obtained results in this case are presented in Figure3.15.

Figure 3 . 15 :

 315 Figure 3.15: Evolution of the salt front for Sc-1, Sc-7, Sc-8 and Sc-9 scenarios, while changing the longitudinal positions of the wells.

Figure 3 .

 3 Figure 3.15 shows the evolution of the interface saltwater/freshwater for different values of x w . (Sc-8) corresponds to the greatest distance from the coast compared to the other simulations, it can be seen that the saltwater is too far from the well position, whereas it is completely polluted in the nearest one ( Sc-9). The saltwater cone is too close in

Figure 3 . 16 :

 316 Figure 3.16: The situation of the Souss-Massa basin and the Souss-Chtouka groundwater.

Figure 3 .Figure 3 . 18 :

 3318 Figure 3.17: 3D representation of the studied domain.

Figure 3 .

 3 Figure 3.19: Left: hydraulic conductivity (D f ) of the aquifer. Right: porosity (φ) of the soil.

Figure 3 . 20 :

 320 Figure 3.20: Left: the freshwater storage coefficient S f . Right: comparison between our developed module 2p-SWI model and the measured piezometric map of 1968.

Figure 3 .

 3 Figure 3.21 presents different values of the topography of the aquifer. The depth of the aquifer varies between -690 m and 35 m where the sea level is used as a reference. Negative values expressed the depth below sea level. The deeper values are located to the northwest of the domain while the higher values are located at the southeast.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[5], 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.

Figure 3 . 21 :

 321 Figure 3.21: Topography of the aquifer bottom.

Figure 3 . 22 :

 322 Figure 3.22: Plan view of the piezometry of the plain Souss-Chtouka before (left) and after (right) solicitation.

  .23).

Figure 3 . 23 :

 323 Figure 3.23: Contours representing the freshwater potential of the Souss-Chtouka plain illustrated on the bottom of the aquifer in 1968 (left) and 2048 (right).

Figure 3 .

 3 Figure 3.24 shows the piezometric of saltwater in its initial state, corresponding to 1968, and after 80 years of service.

Figure 3 . 24 :

 324 Figure 3.24: Plan view of the piezometry of the saltwater in its initial state in 1968 (left) and after 80 years of exploitation (right).

Chapter 3 .

 3 A finite volume method for numerical simulation of the sharp interface model

Figure 3 . 25 :

 325 Figure 3.25: Position of the interface freshwater/saltwater in 1968 (black line) and in 2048 (red line).

Figure 3 . 28 :

 328 Figure 3.28: 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 3 . 29 :

 329 Figure 3.29: Evolution in time of the salt front (top right) freshwater head (left) and saltwater head (right).

Figure 3 . 32 :

 332 Figure 3.32: Geometry and Boundaries of the aquifer (left). The mesh of the studied domain (right).

Figure 3 . 33 :

 333 Figure 3.33: Location of pumping zones (cf. [57]). The rectangle in blue corresponds to the private wells and the public wells are represented by a yellow zone. The yellow zone in the south of the aquifer contain 4 wells in Bhsas while the two other wells are in El-jisr and Malouleh respectively.

Figure 3 . 34 :

 334 Figure 3.34: Evolution of the saltwater/freshwater interface in 2008. "NOP" scenario

Figure 3 . 35 :

 335 Figure 3.35: Evolution of the saltwater/freshwater interface in 2018 with a rate of 150 L/day as recommended by the water authorities.

Figure 3 . 36 :

 336 Figure 3.36: Evolution of the saltwater/freshwater interface in 2018 with a rate of 250 L/day.

Figure 3 . 37 :

 337 Figure 3.37: Evolution of the saltwater/freshwater interface in 2028 with a rate of 250 L/day.

Figure 3 . 38 :

 338 Figure 3.38: Evolution of the piezometric head of freshwater in 2008.

Figure 3 . 39 :

 339 Figure 3.39: Position of the freshwater head scenarios. a) After 10 years of pumping with 150 L/day: "REC" scenario. b) After 10 years of pumping with 250 L/day: "10Y" scenario. c) After 20 years of pumping with 250 L/day: "20Y" scenario.

Figure 4 .

 4 1 and given as follows S f : the storage coefficient of freshwater, b f = u -Z: thickness of the freshwater. b s = Z -Z B : thickness of the saltwater.Q i : quantity of the fluid i = f, s pumped or injected per unit of surface. D: the hydraulic conductivity of freshwater. γ = ρs-ρ f ρ f : the density contrast and ρ i is the density of the fluid i = f, s. F (.): the approximated flux on the saturated/unsaturated and freshwater/saltwater zones. The aquifer is presented by a domain Ω × (Z B , Z T ) with Ω ⊂ R 2 represents a part of its horizontal surface. The function Z B [resp. Z T ] describes its lower [resp. upper] topography.

Figure 4 . 1 :

 41 Figure 4.1: Schematic representation of a coastal aquifer.

Figure 4 . 2 :

 42 Figure 4.2: Convergence of solutions on each grid to the reference solution during 3 days of pumping on cross-section (-50, 0) -(50, 0).

Figure 4 . 3 :

 43 Figure 4.3: Evolution of the solutions on each grid and convergence to the reference solution for 3 days at the point (-10, 0).

Figure 4 .

 4 Figure 4.4 shows the exact total amount of freshwater pumped against the numerical amount. The aquifer is closed. The total amount of freshwater pumped over 3 days is founded exactly using the numerical model. As a consequence, the mass conservation law is satisfied. Numerically, we represent the pumping process in the term source by a negative value and the injection by a positive value.

Figure 4 . 4 :

 44 Figure 4.4: Conservation of the total amount of freshwater pumped

Figure 4 . 7 :

 47 Figure 4.7: Convergence of the interface obtained with the present model towards the analytical solution derived from Ferris model. Times T = 0.5 day (left) and T = 1 day (right).

Figure 4 . 8 :

 48 Figure 4.8: Evolution of the freshwater head and the salt front elevation in differnt times without any external forces. T= 0.1, 0.5 and 1 day.

Figure 4 . 9 :

 49 Figure 4.9: Displacement of the freshwater head and the salt front elevation in different times during a pumping scenario. T=0.1, 0.5, 1 day.

Figure 4 . 11 :

 411 Figure 4.11: Sea oscillations effects on the freshwater head u. Simulations with δ = 0 and δ = 0.

4. 3 . 3 . 1 Figure 4 . 13 :Figure 4 . 14 :

 331413414 Figure 4.13: Initial conditions : initial salt concentration for the 3D model (left) and initial position of the transition zone between freshwater/saltwater for the 2D model (right).

Figure 4 . 15 :

 415 Figure 4.15: Evolution of the freshwater and saltwater flows in the normal conditions during T = 10 days.

Figure 4 . 16 :

 416 Figure 4.16: Representation of sharp-diffuse interfaces and 3D density variable results (50 % salinity contour and 5 % and 95% salinity contour) for T=5 days.

Figure 4 . 17 :

 417 Figure 4.17: The velocity field in the rotating interface problem at t = 10 days.

Figure 4 . 18 :

 418 Figure 4.18: Numerical results for pumping scenario: x 2 = 0.5 (top) x 2 = 1.5 (center) and x 2 = 2.5 (below).

Figure 4 .

 4 Figure 4.18 shows the evolution of the concentration and the salt front in a confined aquifer subjected to pumping. To visualize such a behavior in this case, we made three slices. Then, we consider the slice defined on x 2 = 1.5 and on x 2 = 2.5. The pumping effect appears more in slice x 2 = 1.5. Under the impact of a significant pumping flux, an "upconing" is is formed and tends towards the horizon. The same behavior is shown in the 2D model in both slices. This leads to the conclusion that the 2D model fits well with the realistic 3D model. We perform the same experiment for a higher pumping rate and a higher density contrast. Let Q 2 = -4.0χ [-2,-3]×[1,2]×[-0.5,0] and γ = 0.025. The objective of this experiment is to assess the validity of the 2D sharp-diffuse interfaces model with a large amount of pumping.

Figure 4 . 19 :

 419 Figure 4.19: Evolution of the concentration (3D model) and the salt front (2D model) under a large amount of pumping: slice on x 2 = 1.5.

  [START_REF]DuMu X DUNE for multi-{ Phase, Component, Scale, Physics[END_REF]2]×[START_REF]DuMu X DUNE for multi-{ Phase, Component, Scale, Physics[END_REF]2]×[-0.5,0] and γ = 0.00025 Cherfils and al[START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF] Sharp-diffuse 2D model 3D density variable.

Figure 4 . 20 :Figure 4 .Case 4 :

 42044 Figure 4.20: Numerical results during an injection experiment. x 2 = 0.5 (top), x 2 = 1.5 m (center) and x 2 = 0.5 (bottom)

Figure 4 . 21 :

 421 Figure 4.21: Initial conditions: initial salt concentration for 3D model (left) and initial position of the salt front for the 2D model (right).

Figure 4 .

 4 [START_REF] Cheng | Coastal Aquifer Management-Monitoring, Modeling, and Case Studies[END_REF] presents numerical results for a double pumping wells during 7.5 days.Cherfils and al [23]Sharp-diffuse 2D model 3D density variable.

Figure 4 . 22 :

 422 Figure 4.22: Numerical results for double pumping wells. Slice x 2 = 0.5 m (top), slice x 2 = 1.5 m (center) and slice x 2 = 2.5 m.

  5,-4.5]×[2.3]×[-0.5.0] +0.4χ [4.5.5]×[2.3]×[-0.5.0] . Figure 4.23 represents the numerical results of the corresponding test.Cherfils and al[START_REF] Cherfils | Numerical validation of an upscaled sharp-diffuse interface model for stratified miscible flows[END_REF] Sharp-diffuse 2D model 3D density variable.

Figure 4 . 23 :

 423 Figure 4.23: Numerical results: injection and pumping, T = 10.5 days γ = 0.00025.

φ

  [%] is the porosity of the medium, D [m • day -1 ] the hydraulic conductivity of the freshwater, ρ s , ρ f [kg • m -3 ] the densities of the fresh and the saltwater, γ = ρs-ρ f ρ f [%] the density contrast, δ [m] the width of the transition zone, Z B [m] the bottom of the aquifer, b s (u) = u -Z B [m] the thickness of the aquifer, b s

Proposition 5 . 5 . 2 .

 552 Let assumptions (A.1)-(A.3) be fulfilled. The coupled system (5.3.19)-(5.3.25) admits a discrete solution

2 .

 2 D∇ϕ (u) • ∇ψ dx dtUsing the definition of the discrete gradient on Ω, we getT 2m = T 0 Ω D K ∇ D ϕ u δtm hm • ∇ D ψ δtm hm dx dt.Following the above , the discrete gradient ∇ D ϕ u δtm hm converges weakly towards ∇ϕ(u) in L 2 (0, T ; Ω) Due to the regularity property of the function ψ, ∇ D ψ δtm T hm converges uniformly to ∇ψ. Under the assumption (A.1), we deduce by passing to the limit thatlim m→∞ T 2m = T 0 Ω D∇ϕ (u) • ∇ψ dx dt.In the same way, we obtain the following limits lim m→∞ • ∇ψ dx dt.For the saltwater equation, we follow the same steps as previously. For any ζ ∈ A T , we multiply equation (5.3.20) by ζ n+1 K = ζ(x K , t n+1 ) and equations (5.3.25), (5.3.22) and (5.3.24) by |σ Id Kr |ζ n+1 σ Id Kr . Adding the equalities and summing over K ∈ T hm and over n = 0, • • • , N m -1, we obtain

0 Ω

 0 Due to the strong convergence of Z δtm hm towards Z and the continuity of the function b s , one has b s (Z δtm hm ) → b s (Z) as m → ∞. By using the regularity of the test function ζ and the compactness criterion of ∇ D u δtm hm , we deduce that R 21m 0, as m → ∞. By the compactness estimate of the discrete gradient and the weak convergence of ∇ D u δtm hm ∇u when m → ∞, we obtain R 22m 0, as m → ∞. Db s (Z δtm hm )∇ D u δtm hm • ∇ D ζ δtm T hm dx dt = T Db s (Z)∇u • ∇ζ dx dt, which completes the proof of Theorem 5.6.1.

  . The example considers a 10 m thick free surface aquifer which the horizontal surface is represented by Ω = ]-50 m, 50 m[ × ]-20 m, 20 m[. Homogeneous Neumann boundary conditions are imposed on the boundaries to allow the interfaces evolve freely. For the initial conditions, a fixed hydraulic head of freshwater u = -1 m is imposed and a vertical interface separates the freshwater and saltwater zones.

Figure 5 .Chapter 5 .Figure 5 . 2 :

 5552 Figure 5.2: Structured mesh of a homogeneous aquifer (21164 cells).

Figure 5 . 3 :Figure 5 . 4 :

 5354 Figure 5.3: Displacement of the freshwater/saltwater transition zone at vertical crosssection through pumping wells in the homogeneous case after 3 days (left) and 10 days (right).

Figure 5 . 5 :

 55 Figure 5.5: Displacement of the freshwater/saltwater interface 2D plot in the homogeneous case after 3 days (left) and 10 days (right).

Figure 5 . 6 :

 56 Figure 5.6: Evolution of the freshwater head 2D plot in the homogeneous case after 3 days (left) and 10 days (right).

Figure 5 . 7 :

 57 Figure 5.7: Unstructured mesh in the heterogenous case with two intersecting fractures (20651 cells).

Figure 5 .Chapter 5 .Figure 5 . 8 :

 5558 Figure 5.8: Displacement of the freshwater/saltwater transition zone for different sections in the heterogeneous case after 3 days (left) and 10 days (right).

Figure 5 . 9 :

 59 Figure 5.9: Freshwater flow velocity in the heterogeneous case after 10 days.

Figure 5 . 10 :

 510 Figure 5.10: Displacement of the freshwater/saltwater interface 2D plot in the heterogeneous case after 3 days (left) and 10 days (right).

Figure 5 . 11 :

 511 Figure 5.11: Evolution of the freshwater head 2D plot in the heterogeneous case after 3 days (left) and 10 days (right).

Table of Figures

 of 1.1 Saltwater intrusion into coastal water wells. Source: http://imag.slid.com 2.1 Sharp interface model. Source: https://www.solinst.com . . . . . . . . . . 2.2 Groundwater flow patterns and the dispersion zone in a homogeneous coastal aquifer. Source: http://water.usgs.gov/ . . . . . . . . . . . . . . .

2.3 Schematic representation of a coastal aquifer [71]. . . . . . . . . . . . . . 3.1 Geometry of the domain Ω. . . . . . . . . . . . . . . . . . . . . . . . . . .

  Table of Figures xv 3.35 Evolution of the saltwater/freshwater interface in 2018 with a rate of 150 L/day as recommended by the water authorities. . . . . . . . . . . . . . . 3.36 Evolution of the saltwater/freshwater interface in 2018 with a rate of 250 L/day. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.37 Evolution of the saltwater/freshwater interface in 2028 with a rate of 250 L/day. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 3 . 2 .

 32 Chapter 3. A finite volume method for numerical simulation of the sharp interface model

	Parameters α D f D s	φ	ρ f	ρ s	S f S s
	Values	40 1.0 41.0 0.35 1000 1025 0.0 0.0

Table 3 .

 3 2: Parameters and properties of the aquifer.Different scenarios of the pumping rate and the wells position are presented in Table3.3. Q t represents the amount of water pumped, x w is the distance between the well and the sea while z w is the depth of the well.

	Scenarios	Sc-1 Sc-2 Sc-3 Sc-4 Sc-5 Sc-6 Sc-7 Sc-8 Sc-9 Sc-10 Sc-11
	Q t m 3 /day	0.1	0.05 0.07 0.15 0.1	0.1	0.1	0.1	0.1	0.07	0.07
	x w [m]	150 150 150 150 150 150 200 300 100 150	150
	z w [m]	15	15	15	15	0	25	15	15	15	0	25

Table 3 . 3

 33

: Different pumping scenarios.

Table 3 .

 3 5: Numerical convergence of the scheme in l 2 -relative norm.

	-5

Table 3 .

 3 6: Numerical convergence of the scheme in l ∞ -relative norm.

Table 3 .

 3 7: Order of convergence for the l 2 -relative norm. model

		400 1600 6400 25600 102400
	Freshwater head u -	1.06 1.2	1.3	1.46
	Saltwater head v	-	0.88 1.09 1.23	1.34
	Elevation Z	-	1.02 1.12 1.2	1.3

Table 3 .

 3 8: Order of convergence for the l ∞ -relative norm.Finally, we present in Table3.9 the value of the CPU time corresponding to each simulation.

	Number of cells	400 1600 6400 25600 102400 409600
	Time CPU (min) 0.06 0.3	1.3	4.06	17.46	80.49

Table 3 . 9

 39

: CPU time taken for each simulation.

Table 3 .

 3 10: Well position and associated rates.

		P1	P2	P3	P4	P5	P6
	X [m]	99811	100277 101250 100750 102000 99507
	Y [m]	384285 383966 374280 374280 375400 375246
	Rate m 3 /day	2918.84 2686.99 1011.10 673.31 710.22 1113.03
		P7	P8	P9
		101710 100062 100000
		374431 372311 374996
		1533.89 1402.35 1488.54

Table 3 .

 3 11.

	Parameters	symbol value Unit
	Bottom of the aquifer	Z B	0	[m]
	Porosity of the medeuim φ	0.15	[-]
	Sea level	v 0	200	[m]
	Hydraulic conductivity	D	1.156 [m/day]
	Density of freshwater	ρ f	1000 [kg • m -3 ]
	Density of saltwater	ρ s	1025 [kg • m -3 ]
	Density contrast	γ	0.025 [-]
	Well raduis	--	10;40 [m]
	Storage coefficient	S f	0.0	[-]
	Freshwater inflow	Q inf low 26.0	m 3 /day

Table 3 .

 3

11: Parameters and properties of the aquifer (cf.

Table 3 .

 3 [START_REF] Badaruddin | Water table salinization due to seawater intrusion[END_REF] summarized the associated pumping rates for public wells.

	Wells	Number of wells Pumping rate [m 3 /day] Well raduis [m]
	Private wells -	21 000	-
	Bahsas	4	4680	40
	Malouleh	1	3264	10
	El-Jisr	1	4896	10

Table

Table 4

 4

	.1.

Table 4 .

 4 1: Parameter values for the pumping scenario.

Table 4 .

 4

	.2

2: Different grids used in the numerical convergence.

Table 4 .

 4 3: Numerical convergence of the scheme in l 2 -relative norm.

	Cells	40	160	640	2560	10240
	Freshwater head u 8.35e -2 4.18e -2 1.95e -2 8.29e -3 2.79e -3
	Saltwater head v	5.74e -2 2.73e -2 1.23e -2 5.3e -3	1.8e -3
	Elevation Z	2.03	1.08	0.48	0.2	6.98e -2

Table 4 .

 4 4: Numerical convergence of the scheme in l ∞ norm

	Number of cells	40 160 640 2560 10240
	Freshwater head u -1.09 1.03 1.04 1.15
	Saltwater head v	-1.07 1.03 1.03 1.16
	Elevation Z	-1.02 1.01 1.03 1.16

Table 4 . 5

 45

: Order of convergence

Table 4 .

 4 [START_REF] Abudawia | Analyse numérique d'une approximation élément fini pour un modèle d'intrusion saline dans les aquifères côtiers[END_REF] summarizes the CPU time and the number of iterations calculated on different grids.

	Cells	40	160	640	2560 10240 40960
	Time CPU (s) 0.06 0.249 0.971 4.06 17.46 80.29

Table 4 .

 4 6: CPU time and the number of iterations required by Newton's method to converge

  3.19)-(5.3.25) tends to a weak solution of the continuous problem, as discretization parameters go to zero. Theorem 5.6.1. Let assumptions (A.1)-(A.3) be fulfilled. Let (T hm ) m∈N * , be a sequence of an admissible triangulation of Ω, h m being the size of the mesh T hm and δt m is the time step size. Let (u δtm hm , Z δtm hm ) be a sequence of discrete solutions to the FV scheme (5.3.19)-

Table 5 .

 5 [START_REF] Cooper | A hypothesis concerning the dynamic balance of fresh water and salt water in a coastal aquifer[END_REF][START_REF] Dagan | Solving the problem of local interface upcoming in a coastal aquifer by the method of small perturbations[END_REF][×]-10, -9.5[ + Q 2 χ ]30, 31.5[×]12, 12.5[ , where χ is the characteristic function of a subdomain. The properties and parameters of the aquifer are summarized in Table5.1. 1: Parameters and properties of the aquifer.

	Values	0.025	0.1	0.3	3	20

Parameters

γ [ ] δ [m] φ [ ] Q 1 [m 3 /day] Q 2 [m 3 /day]

Remerciements

Le travail présenté dans ce mémoire rentre dans le cadre d'une co-direction de thèse entre l'Université Moulay Ismaïl, École Nationale Supérieure d'Arts et Métiers, Meknès, Maroc et l'Université de Pau et des Pays de l'Adour, Laboratoire de Mathématiques et leurs Applications, Pau, France. Ce travail a bénéficié du soutien financier du CNRS via l'IRN (International Research Networks) GE2MI (Groupement Euro-Maghrébin de Mathématiques et de leurs Interactions.

We note finally that we encountered difficulties to find reliable and well documented heterogeneous and anisotropic benchmarks. In many articles, some data are missing. We think that a well documented benchmark for a seawater intrusion problem in heterogeneous and anisotropic aquifers would be very useful for the community. Further work on these important issues is in progress.

Level positions

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. In fact, in the concerned area, 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 3.26 (left)). Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

Discretization and notations

For sake of simplicity, we consider an uniform partition {t 0 , • • • , t N } of [0, T ] with δt = t n+1 -t n the time step size and t n = nδt for all n ∈ {0, ..., N }. Furthermore, let (T h ) h>0 be an admissible triangulation in the sense of [START_REF] Eymard | Finite volume methods[END_REF], such that Ω = ∪ K∈T h K. ∂K = K \ K the boundary of K. Let us note by E = E int ∪E ext the interior and exterior edges of the T h . We denote σ = K \ L the interface between the adjacent triangles K and L. Let |K| be the 2-dimensional Lebesgue measure of K and |σ| the 1-dimensional Lebesgue measure of σ. P = (x K ) K∈T h , (x σ ) σ∈E is a set of points in Ω with x K represents the barycentric point of the control volume K while x σ is the point of continuity on the edge σ. V is the set of vertices of the mesh T h and we note by V K the set of vertices of the triangle K. The size of the mesh T h can be given by h = max K∈T h diam (K). We construct for each K ∈ T h three sub-volume controls K r as presented in Figure 5.1. We note P i , P j and P k the vertices of the triangle K. The points P ij , P ki and P jk are the midpoints of the edges (P i P j ), (P i P k ) and (P j P k ) respectively. We denote σ 1 K i (resp. σ 2 K i ) the segment made by the points P i and P ij [resp. P ki ] with a measure noted |σ 1

] the outward unit normal vector to

For any vertex r ∈ {i, j, k}, the sub-volume control K r is the quadrilateral made by x K , P .,r , P r and P r,. as shown in Figure 5.1 . Let us state the following notations:

The discrete space-time solution corresponds to an approximate solution of the continuous problem (P) given by

Following [START_REF] Cao | Convergence of an mpfa finite volume scheme for a two-phase porous media flow model with dynamic capillarity[END_REF], the point P i/2,j is set on the segment (P i P ij ) such [START_REF] Cao | Convergence of an mpfa finite volume scheme for a two-phase porous media flow model with dynamic capillarity[END_REF] by

(5.3.3)

Fully coupled fully implicit finite volume scheme

In this subsection, we provide a description of the fully coupled fully implicit FV scheme for the system (P ) modeling sharp-diffuse interfaces for seawater intrusion in heterogeneous porous media. It is based on the MPFA method for fluxes and diffusive terms and upstream for advective terms with implicit Euler's time discretization. The main idea of the MPFA scheme [3,[START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] is to obtain a consistent approximation of the flux on each edge of the control volume. Moreover, This scheme gives a linear formulation of the flux across the edges.

The main unknowns are the freshwater head (u n+1 K ) K∈T h and the salt elevation (Z n+1 K ) K∈T h for any n ∈ {0, ..., N -1}. By integrating the coupled system (P ) in [t n , t n+1 ] × K, we obtain the following scheme:

with n K,σ is the external normal to σ. Using a fully implicit approach and

Kr

, yields to the following equations

Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

where n σ Id Kr is the external normal to σ Id Kr and the exact diffusive sub-fluxes on each half edge σ Id Kr (r = i, j, k, Id = 1, 2) are given by the following expressions: In addition, the exact convective sub-flux is as follows

Now, the construction of the MPFA scheme requires a consistent approximation of the sub-fluxes on each half edge of the control volume K. However, each vertex r = {i, j, k} is associated with two half edges σ 1 Kr and σ 2 Kr . The exact sub-fluxes (5.3.8)-(5.3.12) are then approximated by the numerical sub-fluxes on each half edge of the control volume. By using the discrete gradient operator (5.3.2)-(5.3.3) on K r in the diffusive term, the numerical sub-fluxes on each inner half edge are given, for all K ∈ T h and around each vertex r = {i, j, k}, by

and

Kr

• n σ Id Kr .

(5.3.16)

For the convective flux, we evaluate the function b s (Z) in the equation (5.3.12) by an upwind scheme. Using the discrete gradient definition in the equation (5.3.2), the numerical convective sub-flux on each half-edge is expressed as follows

L is the neighboring element of K which has in common the edge σ Id Kr . Finally, the scheme of the problem (P ) is expressed by the following nonlinear coupled (5.3.19)-(5.3.25) system:

(5.3.20)

The flux continuity on each half edge reads

Lr

• n σ = 0, (5.3.21)

Lr

• n σ = 0, (5.3.22)

Lr

• n σ = 0, (5.3.23)

Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

Lr

• n σ = 0, (5.3.24)

The discrete initial conditions are given by

(5.3.26)

Discrete maximum principle

In this section, we prove a discrete maximum principle for the approximate solutions freshwater head and the salt front elevation respectively. To do so, we suppose the existence of at least a solution to the numerical scheme. The existence result is proved in the next section. The nonnegativity of the thickness of the freshwater and saltwater issued from the numerical scheme is presented in the following proposition.

Proposition 5.4.1. Let assumptions (A.1)-(A.3) be fulfilled. Then the following inequalities hold:

The approximate freshwater and saltwater thicknesses are non-negative.

Proof. We demonstrate this result by induction. This statement is hold up for n=0. We assume that this property is valid for any k ≤ n, and we will demonstrate the validity of this property for k = n + 1. We select a control volume K such that

we want to show the non-negativity of the saltwater thickness.

To do so, we note by

-Z B and we assume that ω n+1 K < 0. We multiply the saltwater equation (5.3.20) by

We start by proving that the convective term A 1 is positive. Using the upstream values of b s (Z δt h ) up , we get

where we used the fact that the function b s (Z n+1 K ) is extended by 0 whenever

Combining the two cases, we get

Following the same steps we get that A 2 is also positive. As a result

Since ω n+1 K -is negative and using the induction assumption of Z n K , we get

To prove that these terms are positives, we follow the same idea as done in the previous step. To compute B 1 , we multiply equations (5.3.23) and (5. Substrating the first one from the second and adding the result to B 1 , we get

Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

Using the non-decreasing property of the function G * σ Id Kr with respect to u n+1 K -Z n+1 K , we have

This yields to the following inequality

It remains to prove that B 2 is positive. We have

Using the following inequality

Combining both cases, we get that B 2 is positive. As a result, we have

Therefore, we have

This completes the proof of Proposition 5.4.1.

Energy estimates and compactness results

In this section we first present some energy estimates followed by the existence of discrete solutions to the FV scheme and a compactness result which will be used for the convergence of the numerical scheme.

A priori estimates

In this subsection, we present a set of preliminary results which are needed for the analysis of the scheme.

Proposition 5.5.1. Let assumptions (A.1)-(A.3) be fulfilled. Then there exists a constant C > 0 independent of h and δt such that

Proof. First, we multiply equation ( 5 . Adding the equalities and summing over K ∈ T h , we get

with

We use the following formulas from [START_REF] Cao | Convergence of an mpfa finite volume scheme for a two-phase porous media flow model with dynamic capillarity[END_REF],

Kr

, in the last terms, we obtain after factoring

Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

Summing the result over n = 0 to N -1, equation (5.5.1) can be written as

.5.2)

where

Next, using the identity: for any

A 0N can be written as

From the fact that ϕ is a non-decreasing Lipschitz continuous function with λ its Lipschitz constant, we obtain the following estimate for the term A 1N :

In an analogue way, we get

Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

Adding the last inequalities and using equation (5.5.2), we obtain

We now pass to the proof of the discrete gradient estimates for the saltwater equation. We proceed in the same way. We multiply equation (5.

where

Using the identity: for any

we have

Lr

• n σ = 0,

Lr

• n σ = 0.

106 Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

The initial conditions are given by ũ0

This system of equations is linear and has a unique solution since the associated matrix is symmetric and diagonally dominant. Therefore, the function G is well defined in E. Moreover, the functions ϕ and b s are continuous and the diffusion tensor D is symmetric bounded and uniformly elliptic, we obtain that G is a continuous mapping in E. Now, by construction, for any α ∈ [0, 1], the problem (U, Z) = αG (U, Z) has exactly the same solutions as the numerical scheme (5.3.19)-(5.3.25) with αϕ, αu 0 and αZ 0 instead of ϕ, u 0 and Z 0 . Moreover, for any α ∈ [0, 1], we have, αu 0 u 0 , αZ 0 Z 0 , and λ is also a Lipschitz constant for αϕ. Then the same a priori estimate in Proposition 5.5.1 is also satisfied for any α ∈ [0, 1] and for any solution of (U, Z) = αG(U, Z).

Finally, we have verified all the assumptions of the Leray Schauder theorem [START_REF] Cao | Convergence of an mpfa finite volume scheme for a two-phase porous media flow model with dynamic capillarity[END_REF][START_REF] Deimling | Nonlinear functional analysis[END_REF][START_REF] Michel | A finite volume scheme for two-phase immiscible flow in porous media[END_REF]. The function G has then a fixed point and the numerical scheme (5.3.19)-(5.3.25) admits at least a solution. Proposition 5.5.2 is proved.

Space and time translates estimates

In this subsection, we derive compactness tools corresponding to space and time translate estimates for the approximate solutions, which allow the application of Kolmogorov's theorem [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF][START_REF] Eymard | Finite volume methods[END_REF]. We start with the definition of the discrete gradient in Ω as introduced in [START_REF] Cao | Convergence of an mpfa finite volume scheme for a two-phase porous media flow model with dynamic capillarity[END_REF]:

The discrete semi-norm in the space L 2 0, T ; H 1 (Ω) is defined by

. Now we present the following lemma about space translations.

Lemma 5.5.1. Under our standing assumptions, we extend u δt h by 0 outside of [0, T ] × Ω. Then, there exists a constant C independent of the discretization parameters such that the following inequality holds

Proof. The proof can be done by arguments similar to those from Lemma 4.3 in [START_REF] Cao | Convergence of an mpfa finite volume scheme for a two-phase porous media flow model with dynamic capillarity[END_REF], therefore it is omitted here.

The estimate for the time translation is given by the following lemma.

Chapter 5. Numerical analysis of a multi-point flux approximation finite volume scheme for a sharp-diffuse interfaces model

In order to obtain an estimate for A 1 (ρ), we use the local conservation property of the discrete sub-fluxes on the half-edges. We multiply equation (5.

Chapter 6. Conclusion and perspectives

saltwater thicknesses. After that, an existence result of the discrete solutions is established by proving a fixed point to the nonlinear system. In order to satisfy the assumptions of Kolmogorov's theorem for the convergence of the approximate solution, we proved some compactness tools corresponding to space and time translate estimates. Afterwards, we have presented the theorem of the strong convergence of the discrete solutions. Passing to the limit, we have demonstrate that the approximate solution is a weak solution to the continuous problem. The MPFA numerical scheme is implemented and integrated in the DuMu X framework. Let us mention here that all our development and improvement are done in versions 2.12 and 3.0. Our new module is tested and validated in a strongly heterogeneous free aquifer.