# Semianalytical solutions for contaminant transport under variable velocity field in a coastal aquifer

Behshad Koohbor, Marwan Fahs, Behzad Ataie-Ashtiani, Craig T. Simmons, Anis Younes

► To cite this version:

Behshad Koohbor, Marwan Fahs, Behzad Ataie-Ashtiani, Craig T. Simmons, Anis Younes. Semianalytical solutions for contaminant transport under variable velocity field in a coastal aquifer. Journal of Hydrology, 2018, 560, pp.434-450. 10.1016/j.jhydrol.2018.03.048. hal-02625795

HAL Id: hal-02625795

https://hal.inrae.fr/hal-02625795v1

Submitted on 20 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Flinders
UNIVERSITY

Archived at the Flinders Academic Commons:

http://dspace.flinders.edu.au/dspace/

'This is the peer reviewed version of the following article:
Koohbor, B., Fahs, M., Ataie-Ashtiani, B., Simmons, C. T., &
Younes, A. (2018). Semianalytical solutions for
contaminant transport under variable velocity field in a
coastal aquifer. Journal of Hydrology, 560, 434-450.
https://doi.org/10.1016/j.jhydrol.2018.03.048

which has been published in final form at
http://dx.doi.org/10.1016/j.jhydrol.2018.03.048

© 2017 Elsevier. This manuscript version is made
available under the CC-BY-NC-ND 4.0 license:
http://creativecommons.org/licenses/by-nc-nd/4.0/
Research papers

Semianalytical solutions for contaminant transport under variable velocity field
in a coastal aquifer

Behshad Koohbor, Marwan Fahs, Behzad Ataie-Ashtiani, Craig T. Simmons,
Anis Younes

PII: S0022-1694(18)30214-2

DOI: https://doi.org/10.1016/j.jhydrol.2018.03.048

Reference: HYDROL 22681

To appear in: *Journal of Hydrology*

Received Date: 4 October 2017

Revised Date: 15 March 2018

Accepted Date: 17 March 2018

Please cite this article as: Koohbor, B., Fahs, M., Ataie-Ashtiani, B., Simmons, C.T., Younes, A., Semianalytical solutions for contaminant transport under variable velocity field in a coastal aquifer, *Journal of Hydrology* (2018), doi: https://doi.org/10.1016/j.jhydrol.2018.03.048

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

<figure><img src="image_3.png" /></figure>
# Semianalytical solutions for contaminant transport under variable velocity field in a coastal aquifer

Behshad Koohbor¹, Marwan Fahs¹*, Behzad Ataie-Ashtiani²,³, Craig T. Simmons³, Anis Younes¹,⁴,⁵

¹Laboratoire d'Hydrologie et Geochemie de Strasbourg, University of Strasbourg/EOST/ENGAGES, CNRS, 1 rue Blessig 67084 Strasbourg, France.

²Department of Civil Engineering, Sharif University of Technology, PO Box 11155-9313, Tehran, Iran.

³National Centre for Groundwater Research & Training and School of the Environment, Flinders University, GPO Box 2100, Adelaide, South Australia 5001, Australia.
4 IRD UMR LISAH, F-92761 Montpellier, France.

⁵Laboratoire de Modélisation en Hydraulique et Environnement, Ecole Nationale d'Ingénieurs de Tunis, Tunisia

**Submitted to Journal of Hydrology**

*Contact person: Marwan Fahs

E-mail: fahs@unistra.fr
**Abstract**

Existing closed-form solutions of contaminant transport problems are limited by the mathematically convenient assumption of uniform flow. These solutions cannot be used to investigate contaminant transport in coastal aquifers where seawater intrusion induces a variable velocity field. An adaptation of the Fourier-Galerkin method is introduced to obtain semi-analytical solutions for contaminant transport in a confined coastal aquifer in which the saltwater wedge is in equilibrium with a freshwater discharge flow. Two scenarios dealing with contaminant leakage from the aquifer top surface and contaminant migration from a source at the landward boundary are considered. Robust implementation of the Fourier-Galerkin method is developed to efficiently solve the coupled flow, salt and contaminant transport equations. Various illustrative examples are generated and the semi-analytical solutions are compared against an in-house numerical code. The Fourier series are used to evaluate relevant metrics characterizing contaminant transport such as the discharge flux to the sea, amount of contaminant persisting in the groundwater and solute flux from the source. These metrics represent quantitative data for numerical code validation and are relevant to understand the effect of seawater intrusion on contaminant transport. It is observed that, for the surface contamination scenario, seawater intrusion limits the spread of the contaminant but intensifies the contaminant discharge to the sea. For the landward contamination scenario, moderate seawater intrusion affects only the spatial distribution of the contaminant plume while extreme seawater intrusion can increase the contaminant discharge to the sea. The developed semi-analytical solution presents an efficient tool for the verification of numerical models. It provides a clear interpretation of the contaminant transport processes in coastal aquifers subject to seawater
intrusion. For practical usage in further studies, the full open source semi-analytical codes are made available at the website https://lhyges.unistra.fr/FAHS-Marwan.

**Keywords:** Contaminant transport, Coastal aquifer, Seawater intrusion, Semi-analytical solution, Fourier series
# 1. Introduction

Coastal zones are the most populated regions in the world. In Europe, for example, according to Eurostat, 40.8% of the population lives in coastal regions which cover 40% of the European Union territory (Collet and Engelbert, 2013). In many coastal ecosystems, groundwater is considered as the major source of freshwater. However, due to the higher population density and anthropogenic perturbations, coastal aquifers are vulnerable to serious contamination risks and groundwater is increasingly under threat from pollution. The protection of the groundwater reserves demands accurate predictions and assessments of the pollution impact on the aquifers. Therefore, the study of contaminant transport in coastal aquifers (CTCA) is essential for the management of freshwater resources and planning in coastal areas.

Contaminant transport processes in coastal aquifers are inherently complex. Groundwater flow in these aquifers is significantly influenced by various processes at the seaward boundary such as variable density flow, sea-level rise and tidal fluctuations (e.g. Ataie-Ashtiani et al., 1999, 2013; Abd-Elhamid and Javadi, 2011; Bakhtyar et al. 2013; Anwar et al., 2014; Ketabchi et al., 2016, Lu et al., 2015). Coastal aquifers are primarily characterized by variations of groundwater salinity in space and time as they are often subjected to seawater intrusion (Werner et al. 2013). Due to variable density effects, this salinity may affect the flow patterns and consequently the mass transport processes (Volker et al., 2002). CTCA has been investigated using laboratory experiments (Zhang et al., 2002; Chang and Clement, 2013; Liu et al., 2014; Oz et al., 2015). Comprehensive numerical models have been developed to account for the aforementioned complexities (e.g. Ataie-Ashtiani, 2007; Miller et al., 2013). Despite the advances in numerical modeling, the development of analytical solutions is still necessary to verify the computer models and to perform scoping calculations (e.g. Bolster et al., 2007; Chen et al., 2012; Hayek et
al., 2012; Shahkarami et al., 2015; Parker and Kim, 2015). Analytical solutions are ideally suited for the interpretation of the physical processes. They are suitable for conducting sensitivity and uncertainty analysis as well as for parameter estimation because they can be computed more efficiently than numerical solutions. However, due to the complexity of the mathematical model, analytical solutions of contaminant transport can be only obtained under simplified conditions.

For CTCA, and in general for 2D contaminant transport in variable velocity fields, there is a scarcity of analytical solutions (Bolster et al., 2007). Some analytical solutions have been developed for non-uniform velocity by transforming the variable coefficients of the advection-dispersion equation to constant coefficients (Tartakovskiy and Federico, 1997; Tartakovskiy, 2000; Craig and Heidlauf, 2009). These solutions have been obtained for cases dealing with theoretical variable velocity fields. Bolster et al. (2007) developed analytical solutions for contaminant transport in coastal aquifers on the basis of a perturbation expansion in the parameter that describes the coupling between the flow and salt transport processes. The perturbation method may result in poor accuracy when the perturbation parameter is not sufficiently small (Yeh et al., 2010). Semi-analytical methods take advantage of both analytical methods and numerical techniques and circumvent many limitations of the analytical methods allowing them to be applied in more complex problems (Suk, 2016; Ameli and Craig, 2014).

Recently, a semi-analytical solution for solute transport in a coastal aquifer subject to tidal fluctuations was developed by Suk (2017) based on the classical integral transformation technique. In this solution, the saltwater wedge and density differences were not considered.

The Fourier-Galerkin (FG) method has been widely used to derive high accuracy solutions for problems involving fluid flow, mass and/or heat transfer. This method combines the exactness of the analytical methods with an important extent of generality in describing the geometry and
boundary conditions as with the numerical methods (BniLam and Al-Khoury, 2017). In the frame
of flow and transfer in porous media, FG method has been used for density-driven flow (Fahs et
al., 2014; van Reeuwijk et al., 2009), thermal natural convection (Fahs et al., 2015; Nield and
Kuznetsov 2012), double-diffusive convection (Shao et al., 2016) and heat flow in geothermal
systems (BniLam and Al-Khoury, 2017). In the context of coastal aquifers, the FG method has
been used to obtain the semi-analytical solutions of the popular Henry seawater intrusion
problem (Henry, 1964; Ségol, 1994; Simpson and Clement, 2004; Zidane et al., 2012). Recently
Fahs et al. (2016) and Younes and Fahs (2014) have developed efficient implementations of the
FG method to solve the Henry problem with increasing levels of complexity and realism.

In view of the scarcity of analytical solutions for CTCA, the objective of this work was i) to
derive semi-analytical solutions for problems dealing with contaminant transport in coastal
aquifers using the FG method and ii) to take advantage of these solutions to investigate the effect
of seawater intrusion on contaminant transport. The semi-analytical solution is developed for two
contamination scenarios dealing with contaminant leakage from the aquifer top surface and
contaminant migration from a source at the landward boundary. An adaptation of the FG method
is introduced to solve coupled flow, salt transport and contaminant transport equations. Efficient
implementation is presented to deal with discontinuous contaminant boundary conditions. To the
best of our knowledge, such boundary conditions type has never been solved with the FG
method. The Fourier series are used to evaluate several metrics assessing the contaminant
transport such as the discharge flux to the sea, the amount of contaminant persisting in the
groundwater, the area of the contaminated zone and the solute flux as measured by the Sherwood
number. These metrics are used to investigate the effect of seawater intrusion in contaminant
plume migration. Necessary assumptions (i.e. vertical coastal boundary, constant salt
concentration at the sea boundary and constant dispersion tensor) were made to obtain the semi-analytical solutions. Such assumptions are common practice in coastal aquifers modeling. Despite these simplifications, the results of this study highlighted the worthiness of the developed semi-analytical solutions for benchmarking numerical models and provided a clear and new understanding of the processes of plume migration in coastal aquifers.

## 2. Problem description and governing equations

The configuration under consideration is inspired from the Henry problem which has a long history in the investigation of seawater intrusion in coastal aquifers (Henry, 1964). The problem considers a rectangular aquifer of length ($l$) and depth ($d$). The domain is subjected to a constant flux of freshwater on its left vertical edge. The right vertical edge is assumed to be in contact with the sea. The top and bottom horizontal surfaces are impermeable (Figure 1). This is an idealization of coastal aquifers that made the assumption of vertical coastal boundary and neglects the effect of freshwater discharge on the seawater salinity. These assumptions are necessary to obtain the semi-analytical solution. Yang et al. (2013) have shown how the idealized vertical beach can be used to model real coastal aquifers. Despite these idealized assumptions, Henry problem continues to be used as a surrogate for the understanding of seawater intrusion mechanisms (Abarca et al., 2007; Lu et al., 2009; Kerrou and Renard 2010; Yang et al., 2013; Nick et al., 2013; Rajabi and Ataie-Ashtiani, 2014; Riva et al., 2015; Sebben et al., 2015; Fahs et al., 2016).

Two contamination scenarios were studied in this work. The first scenario is defined by considering a contaminant source at the aquifer top surface (Figure 1a). This scenario corresponds to practical applications where the contaminant reaches the aquifer surface by diffusive leakage through the upper confining unit. Typical examples of this scenario could be
the contaminants leaking through the confining unit from waste disposal and landfill sites, septic tanks, fuel storage tanks and pesticide spills (Woumeni and Vauclin, 2006; Kopp et al., 2010; Bakhtyar et al., 2013; Geng et al., 2016). An analytical solution pertaining to a similar scenario in an inland confined aquifer under uniform velocity field has been developed by Emami-Meybodi et al. (2015). In the second scenario, the contaminant source is located upstream of the freshwater flow (Figure 1b). This source may represent leaching from relatively deep waste disposal facility (J.-S. Chen et al., 2016a). This scenario is widely used for testing numerical algorithms of transport in inland aquifers (Younes and Ackerer, 2008; Burdakov et al., 2012; Finkel et al., 2016) and for understanding transport processes (J.-S. Chen et al., 2016b; K. Chen et al., 2016). Analytical solutions for this scenario in an inland aquifer have been presented by Leij and Dane (1990) and J.-S. Chen et al. (2016b). The combined cases of these contamination scenarios, together with seawater intrusion, do not exist in current literature. They are clearly interesting, important and significant physical problems that warrant attention.

For both scenarios, the contaminant is considered as a tracer. Such a configuration was investigated experimentally by Chang and Clement (2013) and analytically by Bolster et al. (2007). This assumption is made in this work to decouple saltwater intrusion and contaminant transport systems which alleviates computationally the Fourier series solution procedure.

Water flow in the domain can be modeled using the continuity equation and Darcy's law. Under steady-state conditions, the governing equations can be written as follows (Guevara Morel et al., 2015):

$$ \nabla \cdot \mathbf{q} = 0 \qquad (1) $$

$$ \mathbf{q} = -\frac{\rho_0 g}{\mu} k \left( \nabla h + \frac{\rho - \rho_0}{\rho_0} \nabla z \right) \qquad (2) $$
where **q** is the Darcy's velocity [$LT^{-1}$], **h** is the equivalent freshwater head [$L$], **g** is the gravitational acceleration [$LT^{-2}$], $\rho$ is the density of the fluid [$ML^{-3}$], $\mu$ is the dynamic viscosity of the fluid [$ML^{-1}T^{-1}$], $k$ is the permeability [$L^2$], $\rho_0$ is the freshwater density [$ML^{-3}$] and $z$ is depth [$L$].

Steady-state conditions are assumed for both saltwater and passive contaminant mass transport. Such a condition is common for the investigation of seawater intrusion as salt transport in coastal aquifers is a long-term natural process. The leakage of contaminants is usually accidental and transient. Steady-state conditions can be reached if the contaminant source is permanent. Practical examples corresponding to this configuration could be the permanent landfills located near shorelines (Cox et al., 1996; Njue et al., 2012). In this work, contaminant transport steady-state conditions are assumed because, at this stage of research, FG method cannot be applied to transient systems. Furthermore, steady-state solutions are relevant to understand the asymptotic behavior of transient problems. The salt and contaminant transport equations can be written as:

$$ \mathbf{q} \nabla s - \nabla \cdot \mathbf{D}^s \nabla s = 0 \qquad (3) $$

$$ \mathbf{q} \nabla c - \nabla \cdot \mathbf{D}^c \nabla c = 0 \qquad (4) $$

where $c$ (resp. $s$) is the relative solute concentration of contaminant (resp. saltwater) [–] which corresponds to the solute concentration normalized by the concentration of the source contamination (resp. seawater). $D^c$ (resp. $D^s$) is the diffusion-dispersion coefficient for contaminant (resp. salt) [$L^2 T^{-1}$].

The saltwater transport equation is coupled with the flow system via the linear mixture density equation as follows:

$$ \rho = \rho_0 + (\rho_1 - \rho_0)s \qquad (5) $$
where $\rho_1$ is the seawater density [$ML^{-3}$].

A constant isotropic dispersion tensor is considered. For the sake of ultimate accuracy, this assumption is necessary to obtain the semi-analytical solution without recourse to numerical integration required for the treatment of the velocity-dependent dispersion terms. It is also necessary to render the solution affordable without excessive computational requirement (Fahs et al., 2016). The assumption of constant dispersion tensor was firstly suggested by Henry (1964) and widely considered in subsequent studies related to the Henry problem. Kalejaiye and Cardoso (2005) have shown that, in the case where gravity is the main driving force (natural convection) this assumption is legitimate for flows where the Rayleigh number is less than 1000.

Bolster et al. (2007) have generalized this assumption to contaminant transport in coastal aquifers. In such a case the Rayleigh number quantifies the relative importance of the buoyancy and dispersive forces within the saltwater wedge. As defined by Kalejaiye and Cardoso (2005) and Bolster et al. (2007), the Rayleigh number is given by:

$$Ra = \frac{(\rho_1 - \rho_0) \cdot k \cdot d \cdot g}{\mu D^s} \qquad (6)$$

The boundary conditions for flow and salt transport are similar to the standard Henry problem.

The contaminant transport is subject to the boundary conditions as in Eqs (7) and (8) for the top surface and landward contamination scenarios, respectively.

$$ \begin{aligned} \frac{\partial c}{\partial x} &= 0, \quad x = 0, \quad x = \ell \\ c &= 0, \quad z = 0 \\ c &= f_{\text{top}}(x), \quad z = d \end{aligned} \qquad (7) $$

$$ \begin{aligned} \frac{\partial c}{\partial z} &= 0, \quad z = 0, \quad z = d \\ c &= f_{\text{land}}(z), \quad x = 0 \\ c &= 0, \quad x = \ell \end{aligned} \qquad (8) $$
where $f_{Top}(x)$ and $f_{Land}(z)$ are functions simulating the distribution of the contaminant source.
They will be introduced and discussed further in the text.

# 3. The semi-analytical solution

## 3.1 Adaptation of the FG method

The semi-analytical solutions of both contamination scenarios are derived using the FG method applied to the stream function form of the governing equations. Specific changes of variables are used for each contamination scenario in order to obtain periodic boundary conditions that are essential for the FG method. Using the stream function, the non-dimensional form of the flow and saltwater transport equations are given as follows (Fahs et al., 2016):

$$ \frac{1}{N_g} \left( \frac{\partial^2 \Psi}{\partial X^2} + \frac{\partial^2 \Psi}{\partial Z^2} \right) - \frac{\partial S}{\partial X} - \frac{1}{\xi} = 0 \qquad (9) $$

$$ \frac{1}{Pe^s} \left( \frac{\partial^2 S}{\partial X^2} + \frac{\partial^2 S}{\partial Z^2} \right) - \frac{\partial \Psi}{\partial Z} \frac{\partial S}{\partial X} + \frac{\partial \Psi}{\partial X} \frac{\partial S}{\partial Z} - \frac{1}{\xi} \frac{\partial \Psi}{\partial Z} - \frac{\partial S}{\partial X} - \frac{1}{\xi} = 0 \qquad (10) $$

where $N_g = \frac{kg(\rho_1 - \rho_0)d}{\mu q_d}$ is the gravity number representing the ratio between free convection (or buoyancy) and advection effects (Riva et al., 2015), $Pe^s = \frac{q_d}{D^s}$ is the saltwater Péclet number quantifying the relative importance of the advection and dispersion mechanisms (Bolster et al., 2007), $q_d$ [$L^2 T^{-1}$] is the recharge of freshwater from the landward boundary into the domain, $\xi = l/d$ is the aspect ratio, $X = \frac{x}{d}$ and $Z = \frac{z}{d}$. $\Psi = \frac{\psi}{q_d} - Z$ and $S = s - \frac{X}{\xi}$ are the new non-dimensional forms of the stream function and salt concentration after the change of variables required to obtain periodic boundary conditions.

Fourier series used to solve Eqs. (9) and (10) are given by:
$$ \Psi = \sum_{m=1}^{Nm} \sum_{n=0}^{Nn} A_{m,n} \sin(m\pi Z) \cos\left(n\pi \frac{X}{\xi}\right) \quad (11) $$

$$ S = \sum_{r=0}^{Nr} \sum_{s=1}^{Ns} B_{r,s} \cos(r\pi Z) \sin\left(s\pi \frac{X}{\xi}\right) \quad (12) $$

where $Nn$, $Nm$, $Ns$ and $Nr$ are the number of Fourier series modes (or the truncation orders) for the stream function and saltwater concentration in the $x$ and $z$ direction, pairwise and respectively. $A_{m,n}$ and $B_{r,s}$ are the Fourier series coefficients.

The Fourier series expansions are substituted into Eqs. (9) and (10) and the resulting equations are multiplied by the trial functions (Fourier modes) and integrated overall the domain. This yields a final system of nonlinear equations with the Fourier coefficients as unknowns.

$$ R_{g=1..Nm,h=0..Nn}^{F} = \frac{1}{N_g} \omega_h \pi^2 \xi \left( \frac{h^2}{\xi^2} + g^2 \right) A_{g,h} - \sum_{r=0}^{Nr} h \Lambda_{g,r} \tilde{B}_{r,h} - \frac{2}{\pi} \Lambda_{g,0} \delta_{h,0} = 0 \quad (13) $$

$$ R_{g=0..Nr,h=1..Ns}^{Ts} = \frac{1}{Pe^s} \omega_g \pi^2 \xi \left( \frac{h^2}{\xi^2} + g^2 \right) B_{g,h} - g \sum_{n=0}^{Nn} \Lambda_{h,n} \tilde{A}_{g,n} - \omega_g \sum_{s=1}^{Ns} s \Lambda_{h,s} B_{g,s} \\ - \frac{\pi}{4} \sum_{m=1}^{Nm} \sum_{n=0}^{Nn} \sum_{r=0}^{Nr} \sum_{s=1}^{Ns} \left( s.m. \eta_{g,m,r} \theta_{h,n,s} - r.n. \kappa_{g,m,r} \lambda_{h,n,s} \right) B_{r,s} A_{m,n} - \frac{2}{\pi} \Lambda_{h,0} \delta_{g,0} = 0 \quad (14) $$

where $R^F$ and $R^Ts$ are the residual vectors for the flow and salt transport equations, respectively.

The coefficients of this system are given in appendix I.

- The surface contamination scenario:

The contaminant transport equation is reformulated by applying a specific change of variables to obtain periodic boundary conditions and by deriving its non-dimensional form:

$$ \frac{1}{Pe^c} \left( \frac{\partial^2 C}{\partial X^2} + \frac{\partial^2 C}{\partial Z^2} + Z \cdot f_{Top}^* \right) - \left( \frac{\partial \Psi}{\partial Z} + 1 \right) \left( \frac{\partial C}{\partial X} + Z \cdot f_{Top} \right) + \frac{\partial \Psi}{\partial X} \left( \frac{\partial C}{\partial Z} + f_{Top} \right) = 0 \quad (15) $$

where $Pe^c = \frac{q_d}{D^c}$ is the Péclet number for contaminant transport, $C$ is the new contaminant concentration after the change of variable, defined by:
$$C = c - Z \cdot f_{Top}(X) \qquad (16)$$

The Fourier series used to solve Eq. (15) is:

$$C = \sum_{k=1}^{Nk} \sum_{l=0}^{Nl} E_{k,l} \sin(k\pi Z) \cos\left(l\pi \frac{X}{\xi}\right) \qquad (17)$$

where $Nl$ and $Nk$ are the number of Fourier series modes (truncation orders) in the $x$ and $z$ direction, respectively. $E_{k,l}$ represents the Fourier series coefficients.

The FG method can give high accuracy solutions with relatively few Fourier modes as it converges exponentially to the exact solution. However, in the case of sharp solutions, the number of Fourier modes should be considerably increased in order to prevent the occurrence of nonphysical oscillations related to the Gibbs phenomenon (Ameli et al., 2013). For the surface contamination scenario, there is a discontinuity in the boundary condition as the source of contamination is imposed at a given interval of the top surface. In order to control the sharpness of the transition, the source is approximated by a smoother function based on a hyperbolic tangent as previously employed in van Reeuwijk et al. (2009). The source function used in the semi-analytical solution is:

$$f_{Top}(X) = 1 - \tanh\left(\left(\frac{X - X_c}{L_c}\right)^2\right) \qquad (18)$$

Where $X_c$ is the dimensionless abscissa of the contaminant source and $L_c$ is the dimensionless source diameter, both normalized by the aquifer depth ($d$).

The justification behind defining this source function is that, depending on $X_c$ and $L_c$, the values of the first derivative can be zero on left and right boundaries ($f'_{Top}(0) = f'_{Top}(\xi) = 0$). Thus, the boundary conditions on the top and bottom surfaces remain periodic. Similar to the
flow and saltwater system, the Fourier expansions of the stream function (Ψ) and contaminant concentration (S) are substituted into Eq. (15). The arising equations are then multiplied by the trial function $\Theta = 4\sin(g\pi Z)\cos(\frac{h\pi X}{\xi})$ (for g=1,..., Nk and h=0,..., Nl) and integrated over the 2D space domain. All the resulting integrals can be evaluated analytically except the integrals corresponding to the contamination source function ($f_{top}(X)$) which are integrated numerically.

Hence the final system writes:

$$R_{g=1..Nk, h=0..Nl}^{Tc} = \frac{1}{Pe^c} \omega_h \pi^2 \xi \left(\frac{h^2}{\xi^2} + g^2\right) E_{g,h} \\ - \frac{\pi}{4} \sum_{m=1}^{Nm} \sum_{n=0}^{Nn} \sum_{k=1}^{Nk} \sum_{l=0}^{Nl} (k.n.\hat{\kappa}_{g,m,k} \hat{\lambda}_{h,n,l} - m.l.\hat{\eta}_{g,m,k} \hat{\theta}_{h,n,l}) E_{k,l} A_{m,n} + \sum_{l=0}^{Nl} l.E_{g,l} \Lambda_{l,h} \\ + \frac{2\pi}{\xi} \sum_{n=0}^{Nn} n \tilde{A}_{g,n} I_1^{BC} - \sum_{m=1}^{Nm} \sum_{n=0}^{Nn} m A_{m,n} \Upsilon_{g,m} I_2^{BC} + \frac{4}{\pi} \frac{(-1)^g}{g} I_3^{BC} = 0 \quad (19)$$

where $I_1^{BC}$, $I_2^{BC}$, $I_3^{BC}$ are the integrals on the x-direction. These integrals and the coefficients of Eq. (19) are given in Appendix I.

- The landward contamination scenario:

For this scenario, the governing non-dimensional differential equation for contaminant transport is given by:

$$\frac{1}{Pe^c} \left[ \frac{\partial^2 C}{\partial X^2} + \frac{\partial^2 C}{\partial Z^2} + f_{Land}'' \left( 1 - \frac{X}{\xi} \right) \right] - \left( \frac{\partial \Psi}{\partial Z} + 1 \right) \left( \frac{\partial C}{\partial X} - \frac{f_{Land}}{\xi} \right) \\ + \frac{\partial \Psi}{\partial X} \left( \frac{\partial C}{\partial Z} + f'_{Land} \left[ 1 - \frac{X}{\xi} \right] \right) = 0 \quad (20)$$

The change of variable applied on C to insure periodic boundary conditions is:

$$C = c - f_{Land}(Z) \left[ 1 - \frac{X}{\xi} \right] \quad (21)$$

The Fourier series expansion used to satisfy the boundary conditions is given by:
$$C = \sum_{k=0}^{Nk} \sum_{l=1}^{Nl} E_{k,l} \cos(k\pi Z) \sin\left(\frac{l\pi X}{\xi}\right) \quad (22)$$

As for the previous scenario, the contaminant source $f_{Land}(Z)$ is approximated using a smooth function in order to control the discontinuity. The contaminant source function is defined as follows (van Reeuwijk et al., 2009):

$$f_{Land}(Z) = 1 - \tanh\left(\left(\frac{Z - Z_c}{L_c}\right)^2\right) \quad (23)$$

where $L_c$ is the diameter and $Z_c$ is the location of the source center.

This function can (depending on $Z_c$ and $L_c$) satisfy the periodic boundary conditions given by
$$f'_{Land}(0) = f'_{Land}(1) = 0.$$

Eq. (21) is treated using the trial function $\Theta = 4\cos(g\pi Z)\sin\left(\frac{h\pi X}{\xi}\right)$ (for $g=0, \dots, Nk$ and $h=1, \dots, Nl$). The final system of equations for landward contamination scenario is:

$$R_{g=0..Nk,h=1..Nl}^{Tc} = \frac{1}{Pe^c} \bar{\omega}_g \pi^2 \xi \left(\frac{h^2}{\xi^2} + g^2\right) E_{g,h} - \bar{\omega}_g \sum_{l=1}^{Nl} l \Lambda_{h,l} E_{g,l} \\ - \frac{\pi}{4} \sum_{m=1}^{Nm} \sum_{n=0}^{Nn} \sum_{k=0}^{Nk} \sum_{l=1}^{Nl} (k \cdot n \cdot \bar{k}_{g,m,k} \bar{\lambda}_{h,n,l} - m \cdot l \cdot \bar{\eta}_{g,m,k} \bar{\theta}_{h,n,l}) E_{g,h} A_{m,n} \quad (24) \\ + 2 \sum_{m=1}^{Nm} \sum_{n=0}^{Nn} (m J_1^{BC} \Lambda_{h,n} + n J_2^{BC} \Sigma_{n,h}) A_{m,n} + \frac{2}{\pi} J_3^{BC} \Lambda_{h,0} - \frac{4}{\pi} \frac{\xi}{h} b^c J_4^{BC} = 0$$

Integrals and coefficients of Eq. (24) are listed in the Appendix I.

## 3.2 Implementation

For the surface contamination scenario, the semi-analytical solution is obtained by solving the system of Eqs (13), (14) and (19). This system is solved sequentially. Thus, Eqs. (13) and (14)
are firstly solved as a nonlinear system to obtain the coefficients $A_{m,n}$ and $B_{r,s}$. This step is
performed using the FORTRAN computer code developed by Fahs et al. (2016). Next, the
coefficients $A_{m,n}$ are substituted in Eq (19). This yields a linear system with the coefficients $E_{g,h}$
as unknowns. The same procedure is applied to the landward contamination scenario by
replacing Eq (19) by Eq (24). Two computer codes are developed in FORTRAN language to
obtain the Fourier series solution for both considered scenarios. In these codes, we implemented
the unifrontal/multifrontal parallel linear solver UMFPACK (Davis, 2004) which has shown
several advantages in the numerical simulation of water resources problems (Miller et al., 2013).
Preliminary runs of these codes have shown that despite the advantages of the linear solver, the
solution is CPU consuming and impractical. Scrutiny of this matter reveals that the most time-
consuming task of the solving procedure is the numerical evaluation of the integrals $I_1^{BC}$, $I_2^{BC}$ and
$I_3^{BC}$ (for the first scenario) and $J_1^{BC}$, $J_2^{BC}$, $J_3^{BC}$ and $J_4^{BC}$ (for the second scenario). These integrals
are independent of the parameters describing the flow and transport processes. Hence, we
evaluated them separately and used them directly in the codes. For the evaluation of the
numerical integrals, we used a specific integration routine from the IMSL libraries (QDAWO)
[http://www.roguewave.com/products-services/imsl-numerical-libraries/fortran-libraries]. This
procedure allows for considerable reduction in the number of operations required for the
numerical integration. For instance, if we consider a hypothetical case where
$Nm = Nn = Nk = Nl = 100$ and we use 1000 points for the numerical integration, the evaluation
of the term involving $I_2^{BC}$ in the final system requires $O(10^{12})$ operations. By applying the
procedure described here, the number of operations can be reduced to be $O(10^9)$. In addition,
the computational efficiency of the semi-analytical method is enhanced by an OPENMP parallel implementation for the evaluation of the matrix and the right-hand side of the system.

## 4. Evaluation of the contaminant transport characteristics

In our analysis, we addressed several practical concerns related to the migration of contaminants in a coastal aquifer. These include, for example, the mass rate of contaminant released in the domain, the location and the area of the contaminated zone and the contaminant discharge to the sea. Thus, we used the Fourier series to evaluate analytically i) the concentration contours, ii) the contaminant plume area ($A_c$), iii) the amount of contaminant remaining in the aquifer ($M_c$), iv) the average Sherwood number ($\overline{Sh}$) and v) the mass flux of contaminant discharge to the sea ($F_c$). $A_c$ is defined by the area of the domain where the contaminant concentration is greater than 0.1. It is calculated by numerical integration. The mass of the contaminant that persists in the aquifer at the steady-state regime is defined by $M_c = \int_0^1 \int_0^\xi c.dXdZ$. It can be expressed in terms of the Fourier coefficients as follows:

$$ M_c = \begin{cases} \frac{\xi}{\pi} \sum_{k=1}^{Nk} \frac{1-(-1)^k}{k} E_{k,0} + \frac{1}{2} \int_0^\xi f_{Top}(X)dX & \text{first scenario (surface)} \\ \frac{\xi}{\pi} \sum_{l=1}^{Nl} \frac{1-(-1)^l}{l} E_{0,l} + \frac{\xi}{2} \int_0^1 f_{Land}(Z)dZ & \text{Second scenario (landward)} \end{cases} \quad (25) $$

The Sherwood number is analogous to the Nusselt number in heat transfer. It represents the ratio of the mass rate transferred due to forced convection to the one transferred by dispersion at the boundary (Xie et al., 2012):

$$ Sh = \frac{h_m}{D^c / d} \quad (26) $$

where $h_m$ is the forced convection mass transfer coefficient.
For the first scenario, the mass transfer between the impermeable layer and the groundwater in the aquifer occurs by forced convection within the solute boundary layer along the aquifer top surface. In this case, the total mass flow rate can be calculated using the mass analogous Newton's cooling law. Conservation between the forced convection and dispersive fluxes at the top surface can be written as follows:

$$ \left| \frac{D^c}{d} \frac{\partial c}{\partial Z} \right|_{Z=1} = \left| h_m (c_{surf} - c_0) \right| \quad (27) $$

where $c_{surf}$ and $c_0$ are the contaminant concentration at the surface and within the entering freshwater.

In our case, we have $c_{surf} = f_{Top}(X)$ and $c_0 = 0$. Hence, using Eq. (27), the local Sherwood number can be expressed in terms of the concentration gradient at the aquifer top surface by:

$$ Sh = \frac{h_m}{D^c / d} = \frac{1}{f_{Top}(X)} \left| \frac{\partial c}{\partial Z} \right|_{Z=1} \quad (28) $$

$Sh$ can be calculated using the Fourier series expansion of the concentration as follows:

$$ Sh = \frac{\pi}{f_{Top}(X)} \sum_{k=1}^{Nk} \sum_{l=0}^{Nl} k(-1)^k E_{k,l} \cos \left( l\pi \frac{X}{\xi} \right) + 1 \quad (29) $$

The average Sherwood number ($\overline{Sh}$) is given by:

$$ \overline{Sh} = \frac{1}{L_c} \int_{X_c - L_c/2}^{X_c + L_c/2} Sh.dX \quad (30) $$

For the second scenario, the contaminant enters the aquifer by advection-dispersion processes. The mass conservation at the landward boundary does not allow $Sh$ to be expressed in terms of concentration. Thus, the Sherwood number cannot be evaluated.
For the surface contamination scenario, the contaminant discharge flux to the sea can be calculated as follows:

$$F_c = \int_{Z_l}^{1} Q_x c|_{x=\xi} dZ = \frac{1}{\pi} \sum_{k=1}^{Nk} \sum_{l=0}^{Nl} (-1)^l \frac{\cos(k\pi Z_l) - (-1)^k}{k} E_{k,l} + \frac{1}{2} \sum_{m=1}^{Nm} \sum_{n=0}^{Nn} \sum_{k=1}^{Nk} \sum_{l=0}^{Nl} m(-1)^{n+l} A_{m,n} E_{k,l} \left[ \frac{\cos(\pi(k+m)Z_l) - (-1)^{k+m}}{k+m} + \frac{\cos(\pi(k-m)Z_l) - (-1)^{k-m}}{k-m} \right] \quad (31)$$

where $Q_x = q_x d/q_d$ is the dimensionless horizontal velocity and $Z_l$ is the coordinate of the point at the seaside boundary separating the freshwater discharge zone and the seawater inland flow zone. $Z_l$ is evaluated analytically using the Fourier series.

In the landward contamination scenario, the contaminant discharge to the sea is given by:

$$F_c = \frac{1}{Pe^c} \int_{0}^{1} \left. \frac{\partial c}{\partial X} \right|_{x=\xi} dZ = \frac{\pi}{Pe^c \xi} \sum_{l=1}^{Nl} (-1)^l E_{0,l} - \frac{1}{Pe^c \xi} \int_{0}^{1} f_{Land}(Z) dZ \quad (32)$$

The integrals in Eqs. (30) and (32) are evaluated numerically using the IMSL library.

## 5. Results and discussions

Three test examples, dealing with different dispersion coefficients, are examined and discussed for each scenario. We based the discussion on test cases with different dispersion coefficients because the sharpness of the solution (and as a consequence its stability) is chiefly controlled by this parameter. In addition, these cases can be relevant to assess the effect of the Péclet number on the contaminant transport processes. The physical parameters used for the three test cases are $d=1m$, $l=2.5m$, $\rho_0=1000 \text{ kg.m}^{-3}$, $\rho_1=1025 \text{ kg.m}^{-3}$, $q_d=6.6 \times 10^{-5} \text{ m}^2 \text{ s}^{-1}$, $k=1.0204 \times 10^{-9} \text{ m}^2 \text{ s}^{-1}$, $\mu=10^{-3} \text{ kg.m}^{-1} \text{ s}^{-1}$, $\epsilon=0.35$ and $g=9.81 \text{ m.s}^{-2}$. These parameters are identical to those used in the original Henry problem, except the domain length. As in Zidane et al. (2012), $l$ is increased to be 2.5m (instead of 2m in the standard Henry problem) to reduce the inland boundary effects
on the semi-analytical solution. We assumed equal values for $D^c$ and $D^s$ as in Bolster et al. (2007). Three values extracted from the literature are considered: $6.601 \times 10^{-6} m^2 s^{-1}$, $1.3202 \times 10^{-6} m^2 s^{-1}$ and $3.3005 \times 10^{-7} m^2 s^{-1}$ (Henry 1964; Zidane et al., 2012; Younes and Fahs, 2014). The corresponding Péclet numbers are 10, 50 and 200, respectively. For the surface contamination scenario, we assumed $X_c = 1.5$ and $L_c = 0.4$. With these values the source can be as close as possible to the sea while satisfying the conditions imposed on the function $f_{Top}(X)$ ($f'_{Top}(0) = 4 \times 10^{-11}$ and $f'_{Top}(\xi) = -10^{-4}$). For the landward contamination scenario, we assumed $Z_c = 0.5$ and $L_c = 0.2$, as it is customary in the literature. The conditions on the function $f_{Land}(Z)$ are satisfied because $f'_{Land}(0) = -f'_{Land}(1) = 0.0003$. For all cases, the coupled flow and saltwater transport equations are solved as in Fahs et al. (2016) with 19,120 Fourier coefficients ($Nm = 40, Nn = 90, Nr = 42$ and $Ns = 360$).

## 5.1 Stability of the semi-analytical solution and effect of the Péclet number

As mention previously, the Fourier series solution may suffer from unphysical oscillations when the solution is relatively sharp. These oscillations are related to the Gibbs phenomenon (Durran, 1999; Ameli et al., 2013; Fahs et al., 2014). In such a case, a large number of Fourier modes should be used to obtain stable concentration contours. The evaluation of the Sherwood number ($\overline{Sh}$) for the first scenario and the discharge flux ($F_c$) for the second scenario involves the first order derivatives of the concentration Fourier series. Accurate values of ($\overline{Sh}$) and ($F_c$) requires more Fourier modes than that required to obtain stable concentration contours. In this work, we assumed that the solution is stable when ($\overline{Sh}$) and ($F_c$) becomes insensitive to the Fourier modes.
- Results for the surface contamination scenario:

For the first test case ($Pe = Pe^c = Pe^s = 10$), we obtained stable concentration contours with 210 Fourier coefficients ($Nk = 10$; $Nl = 20$). The stable value of the Sherwood number requires 1,550 Fourier coefficients ($Nk = 50$; $Nl = 30$). In the second and third test cases ($Pe = 50$ and 200), the concentration distributions are sharper than the first case. The stable solutions have been obtained, respectively, using 4,550 ($Nk = 50$; $Nl = 90$) and 7,770 ($Nk = 70$; $Nl = 110$) Fourier coefficients. The semi-analytical results (contaminant plume, velocity field and saltwater isochlors) are given in Figure 2. This figure displays different behaviours for the plume migration, saltwater wedge and contaminant discharge to the sea, depending on the Péclet number. For the first test case (small Péclet number), Figure 2a shows that the isochlors are relatively close to the contaminant source. Due to the high dispersivity, the contaminant transport occurs within the mixing zone. The contaminant mainly reaches the sea through the saltwater wedge. Figures 2b and 2c show that the increase of the Péclet number leads to a greater penetration of the seawater tongue at the aquifer bottom and a significant recession of the saltwater wedge at the top surface. The increase in $Pe$ is equivalent to the increase of the Rayleigh number ($Ra = N_g \cdot Pe$). This can be interpreted as an intensification of the buoyancy effects against dispersion mechanisms (within the saltwater wedge) and results in a more accelerated downward vertical flow near the sea boundary. Fluid mass conservation and boundary conditions convert the acceleration of the vertical velocity into horizontal flow near the bottom surface and cause greater penetration of the seawater tongue. At the top surface, advection dominates over dispersion. The freshwater discharge washes out the saltwater and prohibits the saltwater penetration in this zone. Contrary to the first test case, for intermediate
and high Péclet number regimes, the contaminant transport does not occur in the saltwater wedge. Indeed, when the plume approaches the saltwater interface it is uplifted by the flow and exits to the sea through the freshwater discharge zone.

Figure 2 also shows that the contaminated zone contracts and becomes nearly horizontal above the saltwater wedge when $Pe$ is increased. This is related to the freshwater flow intensification in front of the contaminant source which can be explained as follows. The increase in Péclet number pushes the saltwater wedge landward. This saltwater wedge forms a barrier narrowing the area through which the freshwater is discharged to the sea. This causes the acceleration of the freshwater flow above the saltwater wedge. As a consequence, the solute boundary shrinks and the contaminated zone contracts.

Table 1 summarizes the ultimate values of the metrics characterizing the contaminant transport ($A_c, M_c, F_c$ and $\overline{Sh}$) for the three test cases. As expected, this table shows a reduction of $A_c$ with the increase in $Pe$, which is consistent with the results shown in Figure 2. Table 1 also indicates that the amount of contaminant remaining in the aquifer ($M_c$) decreases as $Pe$ increases. $M_c$ depends on the contaminant release to the aquifer and the duration required to reach steady-state conditions. The former is attenuated with an increase in $Pe$ (as explained at the end of this paragraph). The duration required to reach stationary conditions depends mainly on the advection processes. It diminishes with an increase in $Pe$.

Table 1 shows that $\overline{Sh}$ increases with $Pe$. In fact, $\overline{Sh}$ is proportional to the concentration gradient. With the increase of $Pe$, the solute boundary layer shrinks and the concentration gradient increases. In general, $\overline{Sh}$ provides a measure of the mass transfer occurring at the surface. Its augmentation is often seen as an increase of the average forced convection mass
transfer coefficient ($\overline{h}_m = D^c \overline{S}h/d$) and indicates the amplification of the mass transfer from the contaminated surface to the aquifer. However, this assumption is not valid in our case because the dispersion coefficient is variable. To further elaborate the effect of $Pe$ on the amount of contaminant entering the domain, we calculated the rate of mass transfer between the contaminated surface and the aquifer. The values for the three test cases are calculated to be 0.09, 0.04 and 0.02, respectively. They indicate that despite the concentration gradient increase with $Pe$, the mass rate of contaminant entering the domain decreases. Finally, Table 1 shows that $F_c$ decreases as $Pe$ increases. This is related to the attenuation of the contaminant release because at steady-state the flux of contaminant discharge to the sea is equal to the mass flux entering the domain.

- Results for the landward contamination scenario

For the three test cases, stable solutions have been obtained with 2,790 ($Nk = 30; Nl = 90$), 6,150 ($Nk = 40; Nl = 150$) and 15,300 ($Nk = 50; Nl = 300$) Fourier coefficients, respectively. The semi-analytical results are plotted in Figure 3. The first observation is that the second scenario leads to more widespread contamination in the aquifer than the first one. For the three Péclet number regimes, the contaminant plume travels within the domain toward the freshwater discharge zone at the upper part of the seaside boundary. Figure 3a shows that the contaminant transport in the saltwater wedge is compelling but it is limited to the zone of low isochlors. For $Pe=10$, unlike the surface contamination scenario, the contaminant does not reach the sea through the saltwater wedge. It is transported by the upward velocity and discharged mainly at the top of the seaside boundary. For the intermediate value of $Pe$, the contaminant concentration contours closely follow the flow structure (Figure 3b). The contaminant transport still occurs in
the saltwater wedge (in contrast to the first scenario) but its occurrence is relatively attenuated.
For the highest value of Pe (test case 3), the contaminant transport does not occur in the
saltwater wedge (Figure 3c).

The ultimate values of $A_c$, $M_c$ and $F_c$ are listed in table 1. The conclusions drawn for the first
scenario concerning the variation of $A_c$, $M_c$ are still applicable. $F_c$ appears to be insensitive to
$Pe$. The behavior of $A_c$ is evident and only the variations of $M_c$ and $F_c$ will be discussed here.

The behavior of $F_c$ can be explained as follows. In the steady-state regime, the rate of contaminant discharge to the sea is equal to the rate of mass entering through the landward boundary which is mainly due to advection (as the concentration gradient is negligible). Hence, at constant gravity number, $F_c$ is slightly sensitive to $Pe$. We recall that, the value of $F_c$ is equal to the dimensionless mass rate of contaminant entering the domain ($= \int_0^1 f_{Land}(Z) dZ = 0.303$) which confirms a good mass balance obtained with the semi-analytical method. Despite the fact that the input and output fluxes are insensitive to $Pe$, $M_c$ slightly decreases. The reason for this variation is related to the time duration required to reach the steady-state condition which (as for the first scenario) decreases as $Pe$ increases.

## 5.2 Comparison against numerical solution: verification, efficiency of the FG implementation and benchmarking issues

Two computer codes have been developed to obtain the semi-analytical solutions for both considered scenarios. Experiments that deals with the conceptual model considered in this paper do not exist in the literature. Hence, to examine the correctness of these codes we compared their results with those obtained using an advanced in-house numerical model (Younes and Ackerer,
2008; Younes et al., 2009). This model is based on specific techniques for both space and time discretization. Its efficiency and accuracy have been highlighted in several applications related to density-driven flow and transport in porous media (Konz et al., 2009, Fahs et al., 2014, 2015, 2016; Shao et al., 2016). To compare transient numerical solutions against steady-state semi-analytical solutions, we performed transient simulations for a long duration. For ultimate accuracy, numerical simulations are performed using a fine mesh of 36,000 regular triangles.

The results of the numerical model have been compared against the semi-analytical solutions for seawater intrusion in Younes and Fahs (2014) and Fahs et al. (2016). For the sake of brevity, we present here the comparison only for the contaminant transport. Figure 4 illustrates the numerical and semi-analytical concentration contours (10%, 30%, 50%, 70% and 90%) for both scenarios. Close agreement between the semi-analytical and numerical solutions can be seen. Table 2 lists the metrics characterizing the contaminant transport obtained using the numerical model. These results also confirm the excellent agreement between the numerical and semi-analytical solutions. This comparison provides compelling evidence on the correctness of the new codes developed for the semi-analytical solutions and more confidence on the robustness of the numerical model. We should mention that the comparison is not only based on visual inspection, which can be somewhat subjective. It is also based on quantitative indicators that provide more rigorous evaluation.

Further numerical simulations are performed to highlight the relevance of the semi-analytical solution in avoiding time and space discretization artifacts. These artifacts may seriously affect the accuracy of the predictive results, the parameter estimation procedure and the sensitivity analysis outcomes (Nassar and Ginn, 2014; Esfandiar et al., 2015). The results (not presented for the sake of brevity) show that in the case of high Péclet number, the numerical solution is
sensitive to the level of mesh refinement. Coarser meshes introduce important numerical dispersion that can lead to significant anomalies when comparing numerical and semi-analytical solutions. Numerical simulations also show that careful attention must be paid to the time step size in order to avoid possible unphysical oscillations. In some cases, the oscillations are related to seawater intrusion because they appear at the interface between the freshwater and saltwater wedge.

Numerical simulations are also relevant to highlight the efficiency of the semi-analytical solution from a computational point of view. We compared the number of Fourier modes required to obtain the stable semi-analytical solution to the number of unknowns of the numerical model required to obtain equivalent accuracy. The results show that the semi-analytical method requires many less unknowns than the numerical model, irrespective of Pe and for both contamination scenarios. The superiority of the semi-analytical method becomes more pronounced as Pe is increased. For instance, in the high Péclet number case (surface contamination scenario), the semi-analytical solution requires about 7000 Fourier coefficients. Equivalent accuracy can be obtained with a numerical solution involving a mesh of about 20,000 elements corresponding to about 90,000 degrees of freedom.

To further highlight the computational efficiency of the semi-analytical method we evaluated the amount of improvement obtained by applying the specific technique implemented in the codes for the evaluation of the numerical integrals (see section 3). We considered, as an example, a case of the first scenario in which about 4,500 Fourier coefficients should be used to obtain the semi-analytical solution. We performed the calculations on a single computer with two Intel Xeon E5-2690 series (2×12 cores) sharing 64GB of memory. The results show that, with the standard procedure, the CPU time is measured to be more than 1 hour. For the same case, the
CPU time with the new technique for the evaluation of numerical integrals reduces to only 12min. This represents a gain of performance of about 450%. With the parallel implementation, the evaluation of this solution takes only 78 sec. This result highlights, on the one hand, the high efficiency of our implementation and on the other hand the flexibility of the semi-analytical method to enable parallel programming.

The developed semi-analytical solutions can be used for validating and testing numerical codes in the case of coupled flow and transport processes or only transport under known flow field conditions. In this context, the results presented in this paper (mainly Figure 4) can be directly used in further comparative studies. Tables 1 and 2 also provide highly accurate quantitative data that can be used for benchmarking purposes. Furthermore, to facilitate the validation and benchmarking tasks for other configurations we provide the full semi-analytical codes on the website https://lhyges.unistra.fr/FAHS-Marwan.

## 5.3 Effect of seawater intrusion on contaminant transport

In this section, we take advantage of the proposed semi-analytical solution to investigate the effect of $N_g$ on the contaminant transport processes. We aim to provide a better understanding of the influence of seawater intrusion on the contaminant plume migration and discharge to the sea. Sixteen runs were carried out for each test case of the contamination scenarios by varying the value of $N_g$. Similarly to Abarca et al. (2007), $N_g$ was varied between 0.625 and 20. In this discussion, for the sake of clarity, the variation of $N_g$ (at constant Pe) will be interpreted as a variation of the aquifer permeability. Thus, the increase of $N_g$ will be seen as an increase of the aquifer permeability. This leads to an intensification of the flow within the saltwater wedge and results in a greater intrusion of the saltwater tongue. The intensification of seawater intrusion, in turn, accelerates the freshwater main flow as it narrows the pathway through which the
freshwater is discharged to the sea. In the following, we will discuss the influence of this velocity variation on the contaminant plume characteristics.

- *Results for the surface contamination scenario*

Figure 5 reports the contaminant concentration contours (20%, 50% and 80%) for three selected values of $N_g$. It shows that, for the three Péclet regimes, the contaminant plume moves up toward the aquifer top surface as $N_g$ is increased. This is related to the shrinking of the solute boundary layer caused by the acceleration of the freshwater flow in front of the contaminant source. The effects of $N_g$ on $A_c, M_c, F_c$ and $\overline{Sh}$ are reported in Figure 6. A decreasing variation of $A_c$ is reported which is consistent with Figure 5. $M_c$ also decreases (Figure 6b) due to the reduction in the time required to reach the steady-state regime. The results reveal a strong correlation between $A_c$ and $M_c$ which fit very well with the bi-exponential regression model given by: $y_0 + A_1e^{-N_g/t_1} + A_2e^{-N_g/t_2}$. The coefficients $y_0, A_1, t_1, A_2$ and $t_2$ are given in Appendix II. This model explains that $A_c$ and $M_c$ are exponentially decaying when $N_g$ is increased but with different rates depending on $N_g$. Indeed, as we can see in Appendix II, $t_1$ is significantly less than $t_2$. This means that for small values of $N_g$, we have $N_g/t_2 \approx 0$. Hence, the exponential term involving $t_1$ dominates the one with $t_2$. The reverse is true for the largest values of $N_g$.

Figure 6c and 6d show that $\overline{Sh}$ and $F_c$ are perfectly correlated. In fact $\overline{Sh}$ measures the mass rate of contaminant release at the inlet which is equal, at steady-state conditions, to the rate of contaminant discharge to the sea. This good correlation highlights the accuracy of the
simulations and gives more confidence in the robustness of the analysis. Only the variation of $\overline{Sh}$ will be discussed here. Figure 6c indicates the existence of two regimes for the evolution of $\overline{Sh}$. We have almost constant value for small gravity number and increasing variation otherwise.

For small gravity number, the seawater wedge is very far from the contamination source. Hence, there is no influence of the seawater intrusion on the contaminant concentration distribution. For large values of $N_g$, the intense seawater intrusion intensifies the freshwater flow in front of the source and increases the concentration gradient. This result indicates that seawater intrusion can increase the contaminant leakage to the aquifer and as a consequence the risk of contaminant discharge to the sea. The regression analysis reveals a good fitting of $\overline{Sh}$ with a second-degree polynomial model, for the regime where $N_g$ is influential. The polynomial model is given by:

$$B_0 + B_1 N_g + B_2 (N_g)^2 . \text{ The coefficients are given in Appendix II.}$$

- Results for the landward contamination scenario

Figure 7 displays the contaminant concentration contours (20%, 50% and 80%) for three selected values of $N_g$. It shows that the plume migration is more sensitive to the seawater intrusion than the first scenario. The plume moves up toward the top surface of the aquifer with an increase in $N_g$. For the smallest values of $N_g$, seawater intrusion affects mainly the low concentration contours near the sea while for the largest values all concentration contours are impacted throughout the entire aquifer.

Figure 8 displays the variations of $A_c$, $M_c$ and $F_c$ as a function of $N_g$. Figure 8a shows a decreasing variation of $A_c$ at constant Péclet number. This indicates that seawater intrusion limits the spread of the contaminant. The same behaviour is recorded for $M_c$ (Figure 8b) which is also
related to the flow acceleration that reduces the duration required to reach the steady-state condition. Figure 8b also shows that, at small gravity number ($N_g < 1$) $M_c$ is almost insensitive to $Pe$. This is logical because, for $N_g < 1$, the seawater intrusion is very weak. Hence, the increase of $Pe$ has a slight impact on the flow velocity and by consequence has limited impact on the duration required to reach the steady-state condition.

The variation of $F_c$ exhibits two regimes depending on $N_g$ (Figure 8c). A better understanding of this behavior can be gained by investigating the effect of $N_g$ on the mass rate of contaminant entering the aquifer which is equal to $F_c$ (under time-stationary conditions). At small gravity numbers, the contaminant mainly enters the domain by advection due to the freshwater recharge flux. Hence $F_c$ is independent of $N_g$. For the extreme cases of seawater intrusion (largest values of $N_g$), the saltwater wedge can reach the landward boundary and generates an uplifting flow along the landward vertical boundary. The vertical component of the velocity can generate a concentration gradient and amplify the contaminant penetration in the aquifer because the diffusive flux becomes a significant component of the contaminant flux. Below the critical $N_g$ (corresponding to the transition between the two regimes of evolution) a further increase in $N_g$ leads to intense seawater intrusion and causes the increase of contaminant mass entering the domain and consequently the increase of $F_c$. Good regression for all metrics can be obtained using the second-degree polynomial model ($B_0 + B_1 N_g + B_2 (N_g)^2$). The coefficients are given in Appendix II. The good regression gives additional confidence in the physical plausibility of the results.
## 6. Conclusion

A semi-analytical solution is derived for the simulation of contaminant transport in a coastal aquifer subject to seawater intrusion. The problem geometry and configuration are inspired from the Henry problem. Two contamination scenarios, dealing with contaminant leakage from a source at the aquifer top surface and aquifer contamination from the landward boundary, are considered. The semi-analytical solutions are developed based on the FG method. An adaptation of this method is developed to deal with discontinuous contaminant boundary conditions. Specific changes of variables have been applied in order to include the contamination source according to the prescribed contamination scenarios. The Fourier series have been used to evaluate several metrics characterizing the contaminant transport. A sound implementation of the FG based on OpenMP parallel programming and an appropriate evaluation of the Fourier series has been developed to optimize the computational requirements of the semi-analytical method. Using this implementation, the semi-analytical solutions have been obtained for example cases dealing with high Péclet numbers. The main inconvenience of the FG solution is related to its rather complex implementation. For practical usage in further studies, the full open source semi-analytical code is made available at the website https://lhyges.unistra.fr/FAHS-Marwan.

The developed semi-analytical solution is compared against the numerical solution of an in-house code. An excellent agreement has been observed which gives confidence on the correctness of both semi-analytical and numerical solutions. This comparison shows that the semi-analytical solutions avoid several numerical artifacts related to space discretization and time integration. This highlights their benefits as useful mathematical tools for validating numerical models of CTCA and in general for contaminant transport under variable velocity field. In this context, the analytical values of the metrics characterizing the contaminant transport
$(A_c, M_c, F_c)$ and $\overline{Sh}$) provide high-quality, high accuracy quantitative indicators suitable for code benchmarking. The high robustness and computational efficiency of the developed semi-analytical method render it suitable for forward analysis (i.e. parameter estimation or global sensitivity analysis) where a large number of solutions must be evaluated.

Furthermore, results show that the semi-analytical solutions developed in this work help to improve our understanding of the contaminant transport mechanisms in coastal aquifers subjected to seawater intrusion. The findings and conclusions drawn are relevant for the management of groundwater pollution in coastal aquifers as the developed test cases can be found in several real applications related to coastal aquifer contamination. Results indicate that, in the case of surface contamination, seawater intrusion leads to a contracted contaminant plume but increases the contaminant discharge to the sea. In the case of landward contamination, moderate seawater intrusion only affects the spatial distribution of the contaminant plume and has no effect on the contaminant discharge to the sea. However, extreme seawater intrusion can increase the contaminant discharge to the sea. Analytical expressions of the contaminant transport metrics in terms of the gravity number (characterizing seawater intrusion) are derived.

The introduced semi-analytical method is not limited to the scenarios discussed in this work. It can be applied to other interesting scenarios of contaminant transport in coastal aquifers and a wide range of applications involving coupled flow and transport processes. We note that the presented semi-analytical solution can be extended to dense pollutant transport and heterogeneous coastal aquifers. Questions about how FG method can be effectively and efficiently used to solve velocity-dependent dispersion transport and/or transient problems is worthy of future inquiry.
Acknowledgement

Marwan Fahs acknowledges the financial support of the Partenariats Hubert Curien (PHC) via the project CEDRE (35369WE). Behzad Ataie-Ashtiani and Craig T. Simmons acknowledge support from the National Centre for Groundwater Research and Training, Australia. Behzad Ataie-Ashtiani also appreciates the support of the Research Office of the Sharif University of Technology, Iran.
References

Abarca, E., Carrera, J., Sánchez-Vila, X., Dentz, M., 2007. Anisotropic dispersive Henry problem. *Advances in Water Resources* 30, 913–926.
doi:10.1016/j.advwatres.2006.08.005

Abd-Elhamid, H.F., Javadi, A.A., 2011. Impact of sea level rise and over-pumping on seawater intrusion in coastal aquifers. *Journal of Water and Climate Change* 2, 19.
doi:10.2166/wcc.2011.053

Ameli, A.A., Craig, J.R., 2014. Semianalytical series solutions for three-dimensional groundwater-surface water interaction. *Water Resources Research* 50, 3893–3906.
doi:10.1002/2014WR015394

Ameli, A.A., Craig, J.R., Wong, S., 2013. Series solutions for saturated-unsaturated flow in multi-layer unconfined aquifers. *Advances in Water Resources* 60, 24–33.
doi:10.1016/j.advwatres.2013.07.004

Anwar, N., Robinson, C., Barry, D.A., 2014. Influence of tides and waves on the fate of nutrients in a nearshore aquifer: Numerical simulations. *Advances in Water Resources* 73, 203–213.
doi:10.1016/j.advwatres.2014.08.015

Ataie-Ashtiani, B., Volker, R.E., Lockington, D.A., 1999. Tidal effects on seawater intrusion in unconfined aquifers. *J. Hydrol.* 216 (1–2), 17–31.
doi:10.1016/S0022-1694(98)00275-3

Ataie-Ashtiani, B., 2007. MODSharp: Regional-scale numerical model for quantifying groundwater flux and contaminant discharge into the coastal zone. *Environmental Modelling & Software* 22, 1307–1315.
doi:10.1016/j.envsoft.2006.07.005

Ataie-Ashtiani, B., Werner, A.D., Simmons, C.T., Morgan, L.K., Lu, C., 2013. How important is the impact of land-surface inundation on seawater intrusion caused by sea-level rise? *Hydrogeology Journal* 21, 1673–1677.
doi:10.1007/s10040-013-1021-0

Bakhtyar, R., Brovelli, A., Barry, D.A., Robinson, C., Li, L., 2013. Transport of variable-density solute plumes in beach aquifers in response to oceanic forcing. *Advances in Water Resources* 53, 208–224.
doi:10.1016/j.advwatres.2012.11.009

BniLam, N., Al-Khoury, R., 2017. A spectral element model for nonhomogeneous heat flow in shallow geothermal systems. *International Journal of Heat and Mass Transfer* 104, 703–717.
doi:10.1016/j.ijheatmasstransfer.2016.08.055

Bolster, D.T., Tartakovsky, D.M., Dentz, M., 2007. Analytical models of contaminant transport in coastal aquifers. *Advances in Water Resources* 30, 1962–1972.
doi:10.1016/j.advwatres.2007.03.007

Burdakov, O., Kapyrin, I., Vassilevski, Y., 2012. Monotonicity recovering and accuracy preserving optimization methods for postprocessing finite element solutions. *Journal of Computational Physics* 231, 3126–3142.
doi:10.1016/j.jcp.2011.12.041

Chang, S.W., Clement, T.P., 2013. Laboratory and numerical investigation of transport processes occurring above and within a saltwater wedge. *Journal of Contaminant Hydrology* 147, 14–24.
doi:10.1016/j.jconhyd.2013.02.005

Chen, J.-S., Hsu, S.-Y., Li, M.-H., Liu, C.-W., 2016a. Assessing the performance of a permeable reactive barrier–aquifer system using a dual-domain solute transport model. *Journal of Hydrology* 543, 849–860.
doi:10.1016/j.jhydrol.2016.11.002

Chen, J.-S., Lai, K.-H., Liu, C.-W., Ni, C.-F., 2012. A novel method for analytically solving multi-species advective-dispersive transport equations sequentially coupled with first-order decay reactions. *Journal of Hydrology* 420–421, 191–204.
doi:10.1016/j.jhydrol.2011.12.001
Chen, J.-S., Liang, C.-P., Liu, C.-W., Li, L.Y., 2016b. An analytical model for simulating two-dimensional multispecies plume migration. Hydrology and Earth System Sciences 20, 733-753. doi:10.5194/hess-20-733-2016

Chen, K., Zhan, H., Zhou, R., 2016. Subsurface solute transport with one-, two-, and three-dimensional arbitrary shape sources. Journal of Contaminant Hydrology 190, 44-57. doi:10.1016/j.jconhyd.2016.04.004

Collet, I., Engelbert, A., 2013. Coastal regions: people living along the coastline, integration of NUTS 2010 and latest population grid [WWW Document]. URL http://ec.europa.eu/eurostat/web/products-statistics-in-focus/-/KS-SF-13-030 (accessed 3.9.17).

Cox, M.E., Hillier, J., Foster, L., Ellis, R., 1996. Effects Of A Rapidly Urbanising Environment On Groundwater, Brisbane, Queensland, Australia. Hydrogeology Journal 4, 30-47. doi:10.1007/s100400050084

Craig, J.R., Heidlauf, T., 2009. Coordinate mapping of analytical contaminant transport solutions to non-uniform flow fields. Advances in Water Resources 32, 353-360. doi:10.1016/j.advwatres.2008.11.013

Davis, T.A., 2004. Algorithm 832: UMFPACK V4.3: an unsymmetric-pattern multifrontal method. ACM Transactions on Mathematical Software 30, 196-199. doi:10.1145/992200.992206

Durran, D.R., 1999. Numerical Methods for Wave Equations in Geophysical Fluid Dynamics, Texts in Applied Mathematics. Springer New York, New York, NY. doi:10.1007/978-1-4757-3081-4

Emami-Meybodi, H., Hassanzadeh, H., Ennis-King, J., 2015. CO2 dissolution in the presence of background flow of deep saline aquifers: CO2 dissolution in the presence of background flow. Water Resources Research 51, 2595-2615. doi:10.1002/2014WR016659

Esfandiar, B., Porta, G., Perotto, S., Guadagnini, A., 2015. Impact of space-time mesh adaptation on solute transport modeling in porous media. Water Resources Research 51, 1315-1332. doi:10.1002/2014WR016569

Fahs, M., Ataie-Ashtiani, B., Younes, A., Simmons, C.T., Ackerer, P., 2016. The Henry problem: New semianalytical solution for velocity-dependent dispersion. Water Resources Research 52, 7382-7407. doi:10.1002/2016WR019288

Fahs, M., Younes, A., Makradi, A., 2015. A Reference Benchmark Solution for Free Convection in A Square Cavity Filled with A Heterogeneous Porous Medium. Numerical Heat Transfer, Part B: Fundamentals 67, 437-462. doi:10.1080/10407790.2014.977183

Fahs, M., Younes, A., Mara, T.A., 2014. A new benchmark semi-analytical solution for density-driven flow in porous media. Advances in Water Resources 70, 24-35. doi:10.1016/j.advwatres.2014.04.013

Finkel, M., Grathwohl, P., Cirpka, O.A., 2016. A travel time-based approach to model kinetic sorption in highly heterogeneous porous media via reactive hydrofacies. Water Resources Research 52, 9390-9411. doi:10.1002/2016WR019147

Geng, X., Boufadel, M.C., Cui, F., 2016. Numerical modeling of subsurface release and fate of benzene and toluene in coastal aquifers subjected to tides. Journal of Hydrology. doi:10.1016/j.jhydrol.2016.10.039

Guevara Morel, C.R., van Reeuwijk, M., Graf, T., 2015. Systematic investigation of non-Boussinesq effects in variable-density groundwater flow simulations. Journal of Contaminant Hydrology 183, 82-98. doi:10.1016/j.jconhyd.2015.10.004
Hayek, M., Kosakowski, G., Jakob, A., Churakov, S.V., 2012. A class of analytical solutions for multidimensional multispecies diffusive transport coupled with precipitation-dissolution reactions and porosity changes. Water Resources Research 48, WR3525. doi:10.1029/2011WR011663

Henry, H., 1964. Effects of dispersion on salt encroachment in coastal aquifer. U.S. Geol. Surv. Water Supply Pap. 1613, C70-C84.

Kalejaiye, B.O., Cardoso, S.S.S., 2005. Specification of the dispersion coefficient in the modeling of gravity-driven flow in porous media. Water Resources Research 41, n/a-n/a. doi:10.1029/2004WR003925

Kerrou, J., Renard, P., 2010. A numerical analysis of dimensionality and heterogeneity effects on advective dispersive seawater intrusion processes. Hydrogeology Journal 18, 55-72. doi:10.1007/s10040-009-0533-0

Ketabchi, H., Mahmoodzadeh, D., Ataie-Ashtiani, B., Simmons, C.T., 2016. Sea-level rise impacts on seawater intrusion in coastal aquifers: Review and integration. Journal of Hydrology 535, 235-255. doi:10.1016/j.jhydrol.2016.01.083

Konz, M., Younes, A., Ackerer, P., Fahs, M., Huggenberger, P., Zechner, E., 2009. Variable-density flow in heterogeneous porous media — Laboratory experiments and numerical simulations. Journal of Contaminant Hydrology 108, 168-175. doi:10.1016/j.jconhyd.2009.07.005

Kopp, A., Binning, P.J., Johannsen, K., Helmig, R., Class, H., 2010. A contribution to risk analysis for leakage through abandoned wells in geological CO2 storage. Advances in Water Resources 33, 867-879. doi:10.1016/j.advwatres.2010.05.001

Leij, F.J., Dane, J.H., 1990. Analytical solutions of the one-dimensional advection equation and two- or three-dimensional dispersion equation. Water Resources Research 26, 1475-1482. doi:10.1029/WR026i007p01475

Liu, Y., Mao, X., Chen, J., Barry, D.A., 2014. Influence of a coarse interlayer on seawater intrusion and contaminant migration in coastal aquifers. Hydrological Processes 28, 5162-5175. doi:10.1002/hyp.10002

Lu, C., Kitanidis, P.K., Luo, J., 2009. Effects of kinetic mass transfer and transient flow conditions on widening mixing zones in coastal aquifers. Water Resources Research 45, n/a-n/a. doi:10.1029/2008WR007643

Lu, C., Xin, P., Li, L., Luo, J., 2015. Seawater intrusion in response to sea-level rise in a coastal aquifer with a general-head inland boundary. Journal of Hydrology 522, 135-140. doi:10.1016/j.jhydrol.2014.12.053

Miller, C.T., Dawson, C.N., Farthing, M.W., Hou, T.Y., Huang, J., Kees, C.E., Kelley, C.T., Langtangen, H.P., 2013. Numerical simulation of water resources problems: Models, methods, and trends. Advances in Water Resources 51, 405-437. doi:10.1016/j.advwatres.2012.05.008

Nassar, M.K., Ginn, T.R., 2014. Impact of numerical artifact of the forward model in the inverse solution of density-dependent flow problem. Water Resources Research 50, 6322-6338. doi:10.1002/2013WR014672

Nick, H.M., Raoof, A., Centler, F., Thullner, M., Regnier, P., 2013. Reactive dispersive contaminant transport in coastal aquifers: Numerical simulation of a reactive Henry problem. Journal of Contaminant Hydrology 145, 90-104. doi:10.1016/j.jconhyd.2012.12.005

Nield, D.A., Kuznetsov, A.V., 2012. The Effect of Strong Heterogeneity and Strong
Throughflow on the Onset of Convection in a Porous Medium: Non-Periodic Global Variation. Transport in Porous Media 91, 927–938. doi:10.1007/s11242-011-9881-z

Njue, C.N., Cundy, A.B., Smith, M., Green, I.D., Tomlinson, N., 2012. Assessing the impact of historical coastal landfill sites on sensitive ecosystems: A case study from Dorset, Southern England. Estuarine, Coastal and Shelf Science 114, 166–174. doi:10.1016/j.ecss.2012.08.022

Oz, I., Shalev, E., Yechieli, Y., Gvirtzman, H., 2015. Saltwater circulation patterns within the freshwater–saltwater interface in coastal aquifers: Laboratory experiments and numerical modeling. Journal of Hydrology 530, 734–741. doi:10.1016/j.jhydrol.2015.10.033

Parker, J.C., Kim, U., 2015. An upscaled approach for transport in media with extended tailing due to back-diffusion using analytical and numerical solutions of the advection dispersion equation. Journal of Contaminant Hydrology 182, 157–172. doi:10.1016/j.jconhyd.2015.09.008

Rajabi, M.M., Ataie-Ashtiani, B., 2014. Sampling efficiency in Monte Carlo based uncertainty propagation strategies: Application in seawater intrusion simulations. Advances in Water Resources 67, 46–64. doi:10.1016/j.advwatres.2014.02.004

Riva, M., Guadagnini, A., Dell'Oca, A., 2015. Probabilistic assessment of seawater intrusion under multiple sources of uncertainty. Advances in Water Resources 75, 93–104. doi:10.1016/j.advwatres.2014.11.002

Sebben, M.L., Werner, A.D., Graf, T., 2015. Seawater intrusion in fractured coastal aquifers: A preliminary numerical investigation using a fractured Henry problem. Advances in Water Resources 85, 93–108. doi:10.1016/j.advwatres.2015.09.013

Ségol, G., 1994. Classic Groundwater Simulations: Proving and Improving Numerical Models. Prentice Hall.

Shahkarami, P., Liu, L., Moreno, L., Neretnieks, I., 2015. Radionuclide migration through fractured rock for arbitrary-length decay chain: Analytical solution and global sensitivity analysis. Journal of Hydrology 520, 448–460. doi:10.1016/j.jhydrol.2014.10.060

Shao, Q., Fahs, M., Younes, A., Makradi, A., Mara, T., 2016. A new benchmark reference solution for double-diffusive convection in a heterogeneous porous medium. Numerical Heat Transfer, Part B: Fundamentals 70, 373–392. doi:10.1080/10407790.2016.1215718

Simpson, M.J., Clement, T.P., 2004. Improving the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models. Water Resources Research 40, n/a-n/a. doi:10.1029/2003WR002199

Suk, H., 2017. Semi-analytical solution of land-derived solute transport under tidal fluctuation in a confined aquifer. Journal of Hydrology 554, 517–531. doi:10.1016/j.jhydrol.2017.09.033

Suk, H., 2016. Generalized semi-analytical solutions to multispecies transport equation coupled with sequential first-order reaction network with spatially or temporally variable transport and decay coefficients. Advances in Water Resources 94, 412–423. doi:10.1016/j.advwatres.2016.06.004

Tartakovsky, D.M., 2000. An analytical solution for two-dimensional contaminant transport during groundwater extraction. Journal of Contaminant Hydrology 42, 273–283. doi:10.1016/S0169-7722(99)00086-8

Tartakovsky, D.M., Di Federico, V., 1997. An Analytical Solution for Contaminant Transport in nonuniform Flow. Transport in Porous Media 27, 85–97. doi:10.1023/A:1006579817668 van Reeuwijk, M., Mathias, S.A., Simmons, C.T., Ward, J.D., 2009. Insights from a
pseudospectral approach to the Elder problem. Water Resources Research 45, n/a-n/a.
doi:10.1029/2008WR007421

Volker, R.E., Zhang, Q., Lockington, D.A., 2002. Numerical modelling of contaminant transport
in coastal aquifers. Mathematics and Computers in Simulation 59, 35–44.
doi:10.1016/S0378-4754(01)00391-3

Werner, A.D., Bakker, M., Post, V.E.A., Vandenbohede, A., Lu, C., Ataie-Ashtiani, B.,
Simmons, C.T., Barry, D.A., 2013. Seawater intrusion processes, investigation and
management: Recent advances and future challenges. Advances in Water Resources 51,
3–26. doi:10.1016/j.advwatres.2012.03.004

Woumeni, R.S., Vauclin, M., 2006. A field study of the coupled effects of aquifer stratification, fluid density, and groundwater fluctuations on dispersivity assessments. Advances in Water Resources 29, 1037–1055. doi:10.1016/j.advwatres.2005.09.002

Xie, Y., Simmons, C.T., Werner, A.D., Diersch, H.-J.G., 2012. Prediction and uncertainty of free
convection phenomena in porous media. Water Resources Research 48, n/a-n/a.
doi:10.1029/2011WR011346

Yang, J., Graf, T., Herold, M., Ptak, T., 2013. Modelling the effects of tides and storm surges on
coastal aquifers using a coupled surface–subsurface approach. Journal of Contaminant
Hydrology 149, 61–75. doi:10.1016/j.jconhyd.2013.03.002

Yeh, H.D., Huang, C.S., Chang, Y.C., Jeng, D.S., 2010. An analytical solution for tidal
fluctuations in unconfined aquifers with a vertical beach. Water Resources Research 46,
n/a-n/a. doi:10.1029/2009WR008746

Younes, A., Ackerer, P., 2008. Solving the advection-dispersion equation with discontinuous Galerkin and multipoint flux approximation methods on unstructured meshes. International Journal for Numerical Methods in Fluids 58, 687–708. doi:10.1002/fld.1783

Younes, A., Fahs, M., 2014. A semi-analytical solution for saltwater intrusion with a very
narrow transition zone. Hydrogeology Journal 22, 501–506. doi:10.1007/s10040-014-1102-8

Younes, A., Fahs, M., Ahmed, S., 2009. Solving density driven flow problems with efficient
spatial discretizations and higher-order time integration methods. Advances in Water
Resources 32, 340–352. doi:10.1016/j.advwatres.2008.11.003

Zhang, Q., Volker, R.E., Lockington, D.A., 2002. Experimental investigation of contaminant
transport in coastal groundwater. Advances in Environmental Research 6, 229–237.
doi:10.1016/S1093-0191(01)00054-5

Zidane, A., Younes, A., Huggenberger, P., Zechner, E., 2012. The Henry semianalytical solution
for saltwater intrusion with reduced dispersion. Water Resources Research 48.
doi:10.1029/2011WR011157
Appendix I. Coefficients of the final systems

$$ \sigma_i = \begin{cases} 2 & \text{if } i=0 \\ 1 & \text{if } i \neq 0 \end{cases} \qquad (I.1) $$

$$ \tilde{A}_{i,j} = \begin{cases} A_{i,j} & 0 \le i \le Nm \\ 0 & i > Nm \end{cases} \qquad (I.2) $$

$$ \tilde{B}_{i,j} = \begin{cases} B_{i,j} & 0 \le j \le Ns \\ 0 & j > Ns \end{cases} \qquad (I.3) $$

$$ \Lambda_{i,j} = \begin{cases} \frac{(-1)^{i+j}-1}{i+j} + \frac{(-1)^{i-j}-1}{i-j} & \text{if } i \neq j \\ 0 & \text{if } i = j \end{cases} \qquad (I.4) $$

$\delta_{i,j}$ is the kronecker symbol

$$ \eta_{i,j,k} = \delta_{j+k,i} + \delta_{k-j,i} + \delta_{j-k,i} \qquad (I.5) $$

$$ \hat{\eta}_{i,j,k} = \delta_{i,j+k} + \delta_{i,k-j} - \delta_{i,j-k} \qquad (I.6) $$

$$ \bar{\eta}_{i,j,k} = \delta_{i,j-k} + \delta_{i,k-j} + \delta_{i,j+k} \qquad (I.7) $$

$$ \kappa_{i,j,k} = -\delta_{j+k,i} + \delta_{k-j,i} + \delta_{j-k,i} \qquad (I.8) $$

$$ \hat{\kappa}_{g,m,k} = \delta_{g,m+k} + \delta_{g,m-k} - \delta_{g,k-m} \qquad (I.9) $$

$$ \bar{\kappa}_{g,m,k} = \delta_{g,m-k} + \delta_{g,k-m} - \delta_{g,k+m} \qquad (I.10) $$

$$ \theta_{i,j,k} = \frac{(-1)^{i+j+k}-1}{i+j+k} + \frac{(-1)^{i+j-k}-1}{i+j-k} + \frac{(-1)^{i-j+k}-1}{i-j+k} + \frac{(-1)^{i-j-k}-1}{i-j-k} \qquad (I.11) $$

$$ \hat{\theta}_{i,j,k} = \frac{(-1)^{i+j+k}-1}{i+j+k} - \frac{(-1)^{i-j+k}-1}{i-j+k} + \frac{(-1)^{i-j-k}-1}{i-j+k} - \frac{(-1)^{i-j-k}-1}{i-j-k} \qquad (I.12) $$

$$ \bar{\theta}_{i,j,k} = - \left( \frac{(-1)^{k+i-j}-1}{k+i-j} + \frac{(-1)^{k-i+j}-1}{k-i+j} + \frac{(-1)^{k+i-j}-1}{k+i+j} + \frac{(-1)^{k-i-j}-1}{k-i-j} \right) \qquad (I.13) $$

$$ \lambda_{i,j,k} = \frac{(-1)^{i+j-k}-1}{i+j-k} + \frac{(-1)^{i-j+k}-1}{i-j+k} - \frac{(-1)^{i+j+k}-1}{i+j+k} - \frac{(-1)^{i-j-k}-1}{i-j-k} \qquad (I.14) $$

$$ \hat{\lambda}_{i,j,k} = \frac{(-1)^{i+j+k}-1}{i+j+k} - \frac{(-1)^{i-j+k}-1}{i-j+k} + \frac{(-1)^{i-k+j}-1}{i-k+j} - \frac{(-1)^{i+k-j}-1}{i+k-j} \qquad (I.15) $$
$$ \bar{\lambda}_{i,j,k} = \left( \frac{(-1)^{k+i-j}-1}{k+i-j} + \frac{(-1)^{k-i+j}-1}{k-i+j} - \frac{(-1)^{k+i+j}-1}{k+i+j} - \frac{(-1)^{k-i-j}-1}{k-i-j} \right) \quad (I.16) $$

$$ \Upsilon_{i,j} = \begin{cases} 2 \left( \frac{(-1)^{j+i}}{j+i} + \frac{(-1)^{i-j}}{i-j} \right) & \text{if } j \neq i \\ \frac{1}{i} & \text{if } j = i \end{cases} \quad (I.17) $$

$$ \Sigma_{i,j} = \frac{1}{\pi} \left( \frac{(-1)^{i+j}-1}{(i+j)^2} - \frac{(-1)^{i-j}-1}{(i-j)^2} \right) \quad (I.18) $$

$$ I_1^{BC} = \int_0^\xi f_{Top}(X) \sin\left(n\pi\frac{X}{\xi}\right) \cos\left(\frac{h\pi X}{\xi}\right) dX \quad (I.19) $$

$$ I_2^{BC} = \int_0^\xi f'_{Top}(X) \cos\left(\frac{h\pi X}{\xi}\right) \cos\left(n\pi\frac{X}{\xi}\right) dX \quad (I.20) $$

$$ I_3^{BC} = \int_0^\xi (Pe^c f_{Top}''(X) - f'_{Top}(X)) \cos\left(\frac{h\pi X}{\xi}\right) dX \quad (I.21) $$

$$ J_1^{BC} = \int_0^1 f_{Land}(Z) \cos(m\pi Z) \cos(g\pi Z) dZ \quad (I.22) $$

$$ J_2^{BC} = \int_0^1 f'_{Land}(Z) \sin(m\pi Z) \cos(g\pi Z) dZ \quad (I.23) $$

$$ J_3^{BC} = \int_0^1 f_{Land}(Z) \cos(g\pi Z) dZ \quad (I.24) $$

$$ J_4^{BC} = \int_0^1 f''_{Land}(Z) \cos(g\pi Z) dZ \quad (I.25) $$

The integrals $I_1^{BC}$, $I_2^{BC}$, $I_3^{BC}$, $J_1^{BC}$, $J_2^{BC}$, $J_3^{BC}$ and $J_4^{BC}$ are evaluated numerically using QDAWO integration routine from the IMSL libraries.
## Appendix II. Coefficients for the regression models

<table><caption>**Table II.1.** Coefficients for the bi-exponential model used for the regression of $A_c$ and $M_c$ (surface contamination scenario).</caption><thead><tr><th><i>Pe</i></th><th><i>Y</i><sub>0</sub></th><th><i>A</i><sub>1</sub></th><th><i>t</i><sub>1</sub></th><th><i>A</i><sub>2</sub></th><th><i>t</i><sub>2</sub></th></tr></thead><tbody><tr><td colspan="6" style="text-align:center;"><i>A</i><sub>c</sub></td></tr><tr><td>10</td><td>-6273.80334</td><td>0.32654</td><td>2.31438</td><td>6274.1983</td><td>1.10×10<sup>6</sup></td></tr><tr><td>50</td><td>0.09419</td><td>0.12038</td><td>1.72974</td><td>0.08971</td><td>11.84595</td></tr><tr><td>200</td><td>0.03512</td><td>0.06188</td><td>1.90672</td><td>0.04985</td><td>22.71617</td></tr><tr><td colspan="6" style="text-align:center;"><i>M</i><sub>c</sub></td></tr><tr><td>10</td><td>-3678.1559</td><td>3678.28691</td><td>3.08962</td><td>0.08493</td><td>2.33×10<sup>6</sup></td></tr><tr><td>50</td><td>-111.85293</td><td>111.90627</td><td>2.67844</td><td>0.04164</td><td>151094.592</td></tr><tr><td>200</td><td>0.01589</td><td>0.01521</td><td>2.15698</td><td>0.01559</td><td>13.78005</td></tr></tbody></table>

<table><caption>**Table II.2.** Coefficients for the second degree polynomial model used for the regression of $F_c$ and $\overline{Sh}$ (surface contamination scenario).</caption><thead><tr><th><i>Pe</i></th><th><i>B</i><sub>0</sub></th><th><i>B</i><sub>1</sub></th><th><i>B</i><sub>2</sub></th></tr></thead><tbody><tr><td colspan="4" style="text-align:center;"><i>Sh</i></td></tr><tr><td>10</td><td>3.824225</td><td>0.307125</td><td>-0.006025</td></tr><tr><td>50</td><td>7.75745</td><td>0.7242</td><td>-0.01385</td></tr><tr><td>200</td><td>16.2514</td><td>1.435125</td><td>-0.02315</td></tr><tr><td colspan="4" style="text-align:center;"><i>F</i><sub>c</sub></td></tr><tr><td>10</td><td>0.07882</td><td>0.00475</td><td>-6.87×10<sup>-5</sup></td></tr><tr><td>50</td><td>0.03324</td><td>0.00285</td><td>-5.22×10<sup>-5</sup></td></tr><tr><td>200</td><td>0.01699</td><td>0.00143</td><td>-2.53×10<sup>-5</sup></td></tr></tbody></table>
<table><caption>**Table II.3.** Coefficients for the second degree polynomial model used for the regression of $A_c$ and $M_c$ (Landward contamination scenario).</caption><thead><tr><th>Pe</th><th>B<sub>0</sub></th><th>B<sub>1</sub></th><th>B<sub>2</sub></th></tr></thead><tbody><tr><td colspan="4">A<sub>c</sub></td></tr><tr><td>10</td><td>2.43973</td><td>-0.0613</td><td>0.00115</td></tr><tr><td>50</td><td>2.33241</td><td>-0.16018</td><td>0.00469</td></tr><tr><td>200</td><td>1.63068</td><td>-0.11858</td><td>0.00359</td></tr><tr><td colspan="4">M<sub>c</sub></td></tr><tr><td>10</td><td>0.75099</td><td>-0.02858</td><td>6.76E-04</td></tr><tr><td>50</td><td>0.78682</td><td>-0.04927</td><td>0.00133</td></tr><tr><td>200</td><td>0.78478</td><td>-0.05704</td><td>0.00159</td></tr><tr><td colspan="4">F<sub>c</sub></td></tr><tr><td>10</td><td>0.31201</td><td>-0.00366</td><td>3.60E-04</td></tr><tr><td>50</td><td>0.30187</td><td>-9.06E-04</td><td>1.98E-04</td></tr><tr><td>200</td><td>0.30528</td><td>-0.00107</td><td>1.22E-04</td></tr></tbody></table>
Figure captions

**Figure 1.** Domain of the studied problem and contamination scenarios: contamination source at aquifer top surface (a) and landward boundary (b).

**Figure 2.** Simultaneous depiction of contaminant plume, velocity field of and saltwater isochlors for surface contamination scenario: (a) test case1 (low Péclet number), (b) test case 2 (intermediate Péclet number) and (c) test case 3 (high Péclet number).

**Figure 3.** Simultaneous depiction of contaminant plume, velocity field of and saltwater isochlors for landward contamination scenario: (a) test case1 (low Péclet number), (b) test case 2 (intermediate Péclet number) and (c) test Case 3 (high Péclet number).

**Figure 4.** Comparison between semi-analytical and numerical solutions for both scenarios: (a) test case1, (b) test case 2 and (c) test case 3.

**Figure 5.** Effect of $N_g$ on the main contamination contours (20%, 50% and 80%) for surface contamination scenario: (a) test case 1, (b) test case 2 and (c) test case 3.

**Figure 6.** Variation of the metrics characterizing the contaminant transport versus the gravity number ($N_g$) and corresponding regression models (surface contamination scenario).

**Figure 7.** Effect of $N_g$ on the main contamination contours (20%, 50% and 80%) for landward contamination scenario: (a) test case 1, (b) test case 2 and (c) test case 3.
**Figure 8.** Variation of the metrics characterizing the contaminant transport versus the gravity number ($N_g$) and corresponding regression models (landward contamination scenario).

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_14.png" /></figure>
Surface contamination scenario

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

Landward contamination scenario
<figure><img src="image_18.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_20.png" /></figure>
<figure><img src="image_21.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(d)</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_27.png" /></figure>
<figure><img src="image_28.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(c)</figcaption></figure>
<table><caption>**Table 1.** Values of $A_c$, $M_c$, $F_c$ and $\overline{Sh}$ calculated analytically for both contamination scenarios</caption><thead><tr><th>Test Case</th><th>Nk</th><th>Nl</th><th>Number of Coefficients</th><th>Contaminated Zone Area (A<sub>c</sub>)</th><th>Mass of contaminant (M<sub>c</sub>)</th><th>Discharge Flux (F<sub>c</sub>)</th><th>Sherwood Number (Sh&#773;)</th></tr></thead><tbody><tr><td colspan="8" style="text-align: center;">Surface contamination scenario</td></tr><tr><td>1</td><td>50</td><td>30</td><td>1,550</td><td>0.43</td><td>0.14</td><td>0.09</td><td>4.95</td></tr><tr><td>2</td><td>50</td><td>90</td><td>4,550</td><td>0.17</td><td>0.06</td><td>0.04</td><td>10.15</td></tr><tr><td>3</td><td>70</td><td>110</td><td>7,770</td><td>0.09</td><td>0.04</td><td>0.02</td><td>20.12</td></tr><tr><td colspan="8" style="text-align: center;">Landward contamination scenario</td></tr><tr><td>1</td><td>30</td><td>90</td><td>2,790</td><td>2.22</td><td>0.65</td><td>0.30</td><td>NC*</td></tr><tr><td>2</td><td>40</td><td>150</td><td>6,150</td><td>1.79</td><td>0.62</td><td>0.30</td><td>NC*</td></tr><tr><td>3</td><td>50</td><td>300</td><td>15,300</td><td>1.20</td><td>0.59</td><td>0.30</td><td>NC*</td></tr></tbody></table>

<footnote>*NC: Not calculated</footnote>
<table><caption>**Table 2.** Values of *Ac*, *Mc*, *Fc* and *Sh* calculated numerically for both contamination scenarios</caption><thead><tr><th>Test Case</th><th>Contaminated Zone Area (A<sub>c</sub>)</th><th>Mass of contaminant (M<sub>c</sub>)</th><th>Discharge Flux (F<sub>c</sub>)</th><th>Sherwood Number (Sh&#773;)</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center;">Surface contamination scenario</td></tr><tr><td>1</td><td>0.43</td><td>0.14</td><td>0.09</td><td>4.94</td></tr><tr><td>2</td><td>0.17</td><td>0.06</td><td>0.04</td><td>10.19</td></tr><tr><td>3</td><td>0.08</td><td>0.03</td><td>0.02</td><td>20.51</td></tr><tr><td colspan="5" style="text-align: center;">Landward contamination scenario</td></tr><tr><td>1</td><td>2.22</td><td>0.65</td><td>0.30</td><td>NC*</td></tr><tr><td>2</td><td>1.80</td><td>0.62</td><td>0.30</td><td>NC*</td></tr><tr><td>3</td><td>1.21</td><td>0.59</td><td>0.30</td><td>NC*</td></tr></tbody></table>

<footnote>*NC: Not calculated</footnote>
Analytical solutions for contaminant transport in variable velocity field are scarce

Semianalytical solution is developed for transport in an aquifer subjected to seawater intrusion

Advantages of the semianalytical solution are evaluated by comparison against numerical
solution

Effect of seawater intrusion on contaminant transport is investigated analytically