Get Complete Project Material File(s) Now! »

## Seawater Intrusion in coastal aquifers

The population density of coastal regions around the world is nearly three times the average global population density (Small and Nicholls, 2003). The inhabitants of coastal regions rely greatly on groundwater resources as freshwater supplies around the world. This reliance is increased in dry countries and regions where the supply of surface water does not suffice the demand of a growing population (esp. close to the urban areas). Anthropogenic activity (e.g. over-pumping) (Pokhrel et al., 2015) and climate-change (Oude Essink et al., 2010; Ranjan et al., 2006) can cause the depletion of groundwater resources. The depletion of groundwater in coastal regions can result in, or further aggravate, hazards such as seawater intrusion. Seawater intrusion is reported to be the most significant threat to the groundwater resort quality in coastal regions (Werner et al., 2013).

Seawater intrusion is the subterranean movement of saltwater into the freshwater resorts due to the effect of density and dispersion (Jiao and Post, 2019; Werner et al., 2013). It has been declared a major threat by environmental protection agencies for decades (e.g. EEA, EPA, USGS, BRGM). This phenomenon can lead to the long-term salinization of freshwater resorts and soil in coastal regions which can endanger agricultural activity and impose great economic losses in these regions. Over-exploitation of groundwater resorts due to anthropogenic activities (e.g. over-pumping) or the local heterogeneity (e.g. stratification and existence of fractures) are of the most influential contributors to the extension of the saltwater wedge (i.e. the intruding saltwater front).

Seawater intrusion is usually tackled by numerical simulation. There are two types of models to simulate seawater intrusion: sharp interface approximation and variable-density flow. In sharp interface approximation, freshwater and saltwater are assumed as two immiscible fluids, having a common interface along which the pressures of the two fluids are continuous. In variable-density flow, there is a transition zone with finite thickness, in which the density of the fluid and the concentration of the salt varies continuously from saltwater to the freshwater (Werner et al., 2013). Variable density flow is usually considered more realistic than sharp interface approximation because it has the privilege of considering the transition zone between the freshwater and saltwater, known as the mixing zone. Only the variable-density flow model provides salinity distribution that can be compared and applied on field measurements (Werner et al., 2013).

To properly model seawater intrusion using variable-density flow model, a coupled system of variable-density flow and solute transport equations need to be discretized and solved. The system of equations is fundamentally nonlinear. Loosely speaking, the salt concentration results in a density difference between saltwater and freshwater. The density difference causes a convective flow which is highly sensitive to the density difference (e.g. on average there is only around 2.5% difference in density between the seawater and freshwater). The flow causes dispersion (i.e. velocity-dependent dispersion) and dispersion changes the concentration gradient and results in a new state of density distribution. These aforementioned interactions subjected to the presence of heterogeneity can result in a complex and in most cases a complicated, mathematical problem.

Numerical modeling of seawater intrusion with realistic assumptions is a challenging task. Common numerical techniques such as Finite Difference and Finite Element methods perform well while solving the flow equation, however, they result in large numerical dispersion when applied to solve the solute transport equation. One solution for such a problem is to refine the grid meshing, but the trade-off between the accuracy of the solution and the computational time consumption limits the refinement process (Werner et al., 2013). Therefore, there is a need for accurate numerical methods with proper performance to model seawater intrusion in a realistic configuration. In this part of the work, we contribute to the development of new semi-analytical solutions that can be used in the assessment of newly developed numerical techniques. We were also interested in seawater intrusion in fractured coastal aquifers. For such a case, it is well known that uncertainties related to fractures characteristics and topology have a significant impact on the model outputs. Thus, we contribute an efficient technique of uncertainty analysis based on surrogate modeling. These contributions are discussed in the next sub-sections.

**An overview of Chapter II: A Generalized Semi-Analytical Solution for the Dispersive Henry Problem**

Analytical solutions are of great interest due to their ability to benchmark new numerical solutions and the insights they provide for understanding the physics of problems and the inter-relation of coexisting parameters. Due to their high performance and accuracy, they are ideally suited for the interpretation of physical processes and performing sensitivity analysis and parameter estimation. Analytical solutions of seawater intrusion are predominantly based on sharp interface approximation in the literature (Bear et al., 1999; Bruggeman, 1999; Kacimov, 2001; Kacimov et al., 2006; Reilly and Goodman, 1985; Werner et al., 2013). For variable-density flow model there is no analytical solution due to the complexities associated with the mathematical model and complex boundary conditions. Semi-analytical solutions are alternatives to analytical solutions. They provide high accuracy and performance as analytical solutions and they are more suitable to deal with nonlinearity and complex boundary conditions, as numerical methods. This Henry problem (Henry, 1964), is the only existing semi-analytical solution for seawater intrusion, with the variable density flow model. After more than five decades of its first introduction, the Henry problem and its variants are still widely used as simplified abstractions of seawater intrusion. The popularity of the Henry problem is due to its simplicity and existence of a closed-form unique semi-analytical solution. The first semi-analytical solution was proposed by H. Henry (1964) and was developed for a high uniform diffusion coefficient with oversimplifying density effects assumptions (Ségol, 1994; Simpson and Clement, 2004, 2003; Voss et al., 2010). Several further studies have contributed more realistic solutions with lower diffusion coefficients or velocity-dependent dispersion. However, the previous extensions of the Henry problem did not deliver a robust analytical solution with the inclusion of hydrodynamic dispersion, anisotropy and heterogeneity. This restricts the applicability of the solution to more realistic configurations.

Our contribution in this context is the development of a new semi-analytical solution, based on the Fourier series solutions, for an anisotropic and stratified (i.e. layered heterogeneity) configuration of the dispersive Henry problem. In addition to velocity-dependent dispersion, which was previously included in the solution by Fahs et al. (2016), the effects of anisotropy and stratification are also added to the solution. The solution is obtained using the Fourier series method. A special model is used to describe hydraulic conductivity–depth in stratified heterogeneity. An efficient technique is developed to solve the flow and transport equations in the spectral space. The results provide a better understanding of the combined effect of anisotropy and stratification on seawater intrusion and quantifying metrics such as saltwater toe and flux. The results also helped to explain about some contradictory results stated in previous works in the literature. The proposed solution can be used for benchmarking purposes or explaining the inter-relation of the parameters associated with the physics of the problem.

**An overview of Chapter III: Semi-analytical solution of contaminant transport in coastal aquifers**

In addition to seawater intrusion, coastal aquifers are susceptible to contamination coming from industrial and anthropogenic activities that are situated near the coastal regions. We were interested in contaminant transport in coastal aquifers. Accurate closed-form solutions of this problem are scarce (Bolster et al., 2007). Part of the reason because contaminant transport in coastal aquifers is inherently complex. This complexity mainly comes from the highly variable velocity field near the sea boundary due to the effect of seawater intrusion. Some analytical solutions have been developed for non-uniform velocity with simplifying assumptions or theoretical velocity fields (Bolster et al., 2007; Craig and Heidlauf, 2009; Tartakovsky and Di Federico, 1997). In Chapter III, a robust implementation of the Fourier series method is developed to deliver a semi-analytical solution for contaminant transport in coastal aquifers under variable-velocity field.

This work is inspired by the implementation of Fourier series method on density-driven flow (Fahs et al., 2014; van Reeuwijk et al., 2009). Two scenarios of contaminant sources are considered in this problem. The first scenario is defined by considering a contaminant source at the aquifer top surface. In the second scenario, the contaminant source is located upstream of the freshwater influx boundary. The considered configuration, the domain and the seawater intrusion boundary conditions are inspired by the Henry problem.

**An overview of Chapter IV: Uncertainty analysis for seawater intrusion in fractured coastal aquifers: application to Clashnessie Bay, UK**

Fractured Coastal Aquifers are found globally. Several examples can be found France, USA, UK and all around the Mediterranean region (Arfib and Charlier, 2016; Bakalowicz et al., 2008; Chen et al., 2017; Comte et al., 2018; Doukou and Karatzas, 2012; Perriquet et al., 2014; Xu et al., 2018). Fractures play an important role in controlling seawater intrusion as they represent preferential pathways for water flow and can enhance seawater intrusion (Bear et al., 1999). However, fractures properties and topology are usually uncertain. These uncertainties lead to uncertainty in the model output and question the reliability of its predictive results. Global sensitivity analysis can be applied as an effective step towards understating uncertainty propagation through models. In seawater intrusion, global sensitivity analysis has been applied to study the effects of hydrodynamic parameters in homogeneous coastal aquifers (Herckenrath et al., 2011; Rajabi et al., 2015a; Rajabi et al., 2015b; Rajabi and Ataie-Ashtiani, 2014; Riva et al., 2015; Dell’Oca et al.,2017). However, to the best of our knowledge, global sensitivity analysis has never been applied to seawater intrusion in fractured coastal aquifers. In this study, we present an efficient approach to evaluate uncertainty related to fractures on model outputs. This approach is based on the use of the Polynomial Chaos Expansion as a surrogate model of seawater intrusion. We analyze the uncertainties imposed by some key parameters associated with fractured coastal aquifers (e.g. fracture location, density, hydraulic conductivity, aperture and dispersivity) on the seawater intrusion and relevant metrics (e.g. saltwater toe, salinized area and intruding saltwater flux).

A coupled approach based on a Discrete Fracture Network model and variable-density flow is implemented with COMSOL Multiphysics® to act as model input for generating the metamodel used for the uncertainty quantification. The metamodel is generated by Polynomial Chaos Expansion. In addition to the conventional full order polynomial chaos expansion, we used a sparse polynomial chaos expansion algorithm that allows for reducing the number of simulations required to generate the metamodel. The uncertainty measures used to quantify and discuss the sensitivity of the model output to the model inputs are Sobol’ indices. Sobol’ Indices are used to identify the key variable imposing uncertainty on the various seawater intrusion metrics. Spatial analysis of the Sobol’ indices also helps in locating the most sensitive parts of domain salinity distribution in regards to the input parameters. Additionally, marginal effects are calculated to further understand the contribution of each parameter to the model output. The findings in this chapter can contribute to various environmental applications such as preliminary assessment for field data collection and monitoring purposes associated with seawater intrusion in fractured coastal aquifers. They also show the potential for quality control and management of fresh groundwater resorts in coastal aquifers. The results related to this work have been published as a journal article in the Journal of Hydrology (Koohbor et al., 2019).

The developed strategy for uncertainty analysis has been applied for a hypothetical case. Further investigation of this strategy is its applicability to field studies. The field chosen for the application is in Clashnessie Bay (Scotland, UK) which is a fractured basement-rock coastal aquifer (Lewisian gneiss complex aged Precambrian) intersected by a densely fractured regional fault zone. In the context of this Ph.D., we started this application by developing a 3D seawater intrusion model for the aquifer. Geophysical data (i.e. electrical resistivity tomography) have been used to estimate variations in fracture density and certain equivalent hydrodynamic properties. The model is developed with a discrete fracture approach. This work is a preliminary step to global sensitivity application. It has been developed during my research visit to the School of Geosciences at the University of Aberdeen.

**Variably saturated flow in fractured domains: Application to El Assal aquifer in Lebanon**

Accurate modeling of variably saturated flow in fractured porous media with relatively low computational time remains a challenge. This challenge is partly due to the nonlinear nature of Richards’ Equation (RE), used to describe the flow in unsaturated porous media, and partly due to the complexity in heterogeneity imposed by fractures as conduits. The numerical solution of RE, even in unfractured domains, is among the most challenging problems in hydrogeology (Miller et al., 2013). Mathematical characteristics and strong nonlinearity induced by the constitutive relationships (Suk and Park, 2019) impose numerical difficulties in the process of discretization/resolving RE using conventional methods. The existence of fracture networks can add-up to this challenge. Therefore, numerical solutions of RE in explicitly fractured domains are not well-developed in the literature.

**Table of contents :**

**Chapter I: Introduction **

1.1. Overview

1.2. Seawater Intrusion in coastal aquifers

1.2.1 An overview of Chapter II: A Generalized Semi-Analytical Solution for the Dispersive Henry Problem

1.2.2. An overview of Chapter III: Semi-analytical solution of contaminant transport in coastal aquifers

1.2.3. An overview of Chapter IV: Uncertainty analysis for seawater intrusion in fractured coastal aquifers: application to Clashnessie Bay, UK

1.3. Variably saturated flow in fractured domains: Application to El Assal aquifer in Lebanon

**Chapitre I : Introduction **

1.1. Prolégomènes

1.2. Intrusion d’eau de mer dans les aquifères côtiers

1.2.1 Aperçu du chapitre II : solution semi-analytique généralisée au problème de dispersion de Henry

1.2.2. Aperçu du chapitre III : solution semi-analytique du transport de contaminants dans les aquifères côtiers

1.2.3. Aperçu du chapitre IV : analyse d’incertitudes pour l’intrusion d’eau de mer dans les aquifères côtiers fracturés : application à la baie de Clashnessie, Royaume-Uni

1.3. Ecoulements variablement saturés dans les domaines fracturés : application à l’aquifère d’El Assal au Liban

**Chapter II: Dispersive Henry problem: a generalization of the semianalytical solution to anisotropic and layered coastal aquifers **

2.1. Introduction

2.2. The mathematical model and boundary conditions

2.3. Semianalytical solution

2.4. New technique for solving the equations in the spectral space

2.5. Results and discussions

2.5.1. Verification: Stability of the Fourier series solution and comparison against numerical solution

2.5.2. Effect of anisotropy on seawater intrusion in homogenous aquifer

2.5.3. Coupled effect of anisotropy and stratified heterogeneity on seawater intrusion

2.6. Conclusion

**Chapter III: Semi-analytical solutions for contaminant transport under variable velocity field in a coastal aquifer **

3.1. Introduction

3.2. Problem description and methodology

3.3. Governing Equations

3.4. The semi-analytical solution

3.4.1. Adaptation of the FG method

3.4.2. Implementation

3.5. Evaluation of the contaminant transport characteristics

3.6. Results: test examples, verification and comparison against numerical solution

3.6.1. Stability of the semi-analytical solution and effect of Péclet number

3.6.2. Comparison against numerical solution: verification, efficiency of the FG implementation and benchmarking issues

3.7. Effect of seawater intrusion on contaminant transport

3.8. Conclusion

**Chapter IV: Uncertainty analysis for seawater intrusion in fractured coastal aquifers: Effects of fracture location, aperture, density and hydrodynamic parameters **

4.1. Introduction

4.2. Material and methods

4.2.1. Conceptual model: Fractured Henry Problem

4.2.2. DFMM-VDF mathematical model:

4.2.3. DFMM-VDF finite element model: COMSOL Multiphysics®:

4.2.4. Metrics Design:

4.3. Global sensitivity analysis

4.3.1. Sobol’ indices

4.3.2. Polynomials Chaos Expansion (PCE)

4.3.3. Sparse polynomial chaos expansion

4.4. Validations: COMSOL model and Boussinesq approximation

4.5. Global sensitivity Analysis: results and discussion

4.5.1. The single horizontal fracture configuration (SHF)

4.5.2. The network of orthogonal fractures configuration (NOF)

4.6. Conclusion

4.7. Field case study: Application to Clashnessie Bay, UK

**Chapter V: An advanced discrete fracture model for variably saturated flow in fractured porous media**

5.2. Governing Equations of VSF in fractured domains

5.3. Numerical solution: MHFE method, ML technique, and MOL

5.3.1. Trial functions of the MHFE method for the matrix and fractures

5.3.2. Flux discretization in the matrix elements (MHFE method and ML technique)

5.3.3. Flux discretization for a fracture element

5.3.4. Hybridization: mass conservation on edges

5.3.5. The new ML technique for fractures

5.3.6. High-order adaptive time integration

5.4. Results: Verification and advantages of the new developed numerical scheme

5.4.1. Verifications: Fractured Vauclin test case

5.4.2. Advantages of the ML technique for fractures

5.5. Results: Field Scale Applications

5.5.1. Objectives and overall presentation of the site

5.5.2. Methodology for simulation and analysis

5.5.3. Simulations with real data from 2013 to 2019

5.5.4. Predictions under climate change from 2019 to 2099

5.6. Conclusion

**Chapter VI: Conclusion and perspectives **

6.1. General conclusion

6.2. Perspective

**Chapitre VI : Conclusion et perspectives **

6.1. Conclusion générale

6.2. Perspectives

**References**