A Generalized Semi-Analytical Solution for the Dispersive Henry Problem: Effect of Stratification and Anisotropy on Seawater Intrusion

Abstract

The Henry problem (HP) continues to play a useful role in theoretical and practical studies related to seawater intrusion (SWI) into coastal aquifers. The popularity of this problem is attributed to its simplicity and precision to the existence of semi-analytical (SA) solutions. The first SA solution has been developed for a high uniform diffusion coefficient. Several further studies have contributed more realistic solutions with lower diffusion coefficients or velocity-dependent dispersion. All the existing SA solutions are limited to homogenous and isotropic domains. This work attempts to improve the realism of the SA solution of the dispersive HP by extending it to heterogeneous and anisotropic coastal aquifers. The solution is obtained using the Fourier series method. A special hydraulic conductivity–depth model describing stratified heterogeneity is used for mathematical convenience. An efficient technique is developed to solve the flow and transport equations in the spectral space. With this technique, we show that the HP can be solved in the spectral space with the salt concentration as primary unknown. Several examples are generated, and the SA solutions are compared against an in-house finite element code. The results provide high-quality data assessed by quantitative indicators that can be effectively used for code verification in realistic configurations of heterogeneity and anisotropy. The SA solution is used to explain contradictory results stated in the previous works about the effect of anisotropy on the saltwater wedge. It is also used to investigate the combined influence of stratification and anisotropy on relevant metrics characterizing SWI. At a constant gravity number, anisotropy leads to landward migration of the saltwater wedge, more intense saltwater flux, a wider mixing zone and shallower groundwater discharge zone to the sea. The influence of stratified heterogeneity is more pronounced in highly anisotropic aquifers. The stratification rate and anisotropy have complementary effects on all SWI metrics, except for the depth of the discharge zone.

1. Introduction

The Henry problem (HP) [1] is widely used as a surrogate for the understanding of seawater intrusion (SWI) processes in coastal aquifers [2]. It is an abstraction of SWI in a vertical cross-section of a confined coastal aquifer perpendicular to the shoreline. In this aquifer, an inland freshwater flow is in equilibrium with the seawater intruded from the seaside, due to its higher density (Figure 1). The aquifer is assumed to be homogenous and isotropic.

Recently, a variation of the HP has been used to understand SWI in fractured coastal aquifers [4,5]. A reactive HP has been presented to investigate the interplay between dispersion and reactive processes in coastal aquifers [6] and to study the effect of calcite dissolution on SWI [7,8]. HP has been considered in several works investigating the effects of heterogeneity [9,10,11], anisotropy [12,13], mass dispersion [13], salinity dependent permeability [14], seafloor slope [15] and inland boundary conditions [16] on the SWI extent. Three variants of this problem were used in Post et al. [17] to interpret the groundwater ages in coastal aquifers. Javadi et al. [18] developed a multi-objective optimization algorithm and applied it to the HP in order to assess different management methods for controlling SWI. Hardyanto and Merkel [19], Herckenrath et al. [20], Rajabi and Ataie-Ashtiani [21], Rajabi et al. [22,23], Burrows and Doherty [24] and Riva et al. [25] investigated the HP and conducted uncertainty analyses to evaluate the effect of incomplete knowledge of the aquifer parameters on the predictions of SWI models. Sanz and Voss [26] developed an a posteriori parameter covariance matrix analysis of the HP and provided insight about parameter sensitivities that can help in improving the inverse modeling procedure. A parameter sensitivity analysis based on the HP was also reported in Carrera et al. [27], who presented a general analysis of the inverse problem for SWI.

This brief review points out the importance of the HP for theoretical and practical investigations of conceptual and modeling aspects of SWI. Nonetheless, the popularity of this problem stems precisely from the existence of a semi-analytical (SA) solution [1]. The mathematical formulation of the HP is based on the density-driven flow (DDF) model that incorporates a system of a variable density-flow equation, an advection–dispersion equation and a state relation, expressing the mixture density as function of the saltwater concentration. Because of the existence of this SA solution, HP has been accepted as one of the primary benchmarks for the assessment of DDF numerical codes, e.g., [6,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45].

The first SA solution of the HP was developed by Henry [1] using the Fourier series method. This solution underwent a revision by Segol [46], who pointed out some errors in Henry’s solution and proposed appropriate corrections. With this revision, Segol [46] was able to recreate the SA solution numerically. The worthiness of this SA solution for benchmarking DDF codes has been widely discussed in the literature [2,47,48,49,50,51,52,53]. These studies concluded that Segol’s solution is not sufficient to test DDF numerical models since the buoyancy effects are dominated by the influence of boundary forcing exerted by the freshwater inland flow and by the diffusion effects, as the solution was developed for exaggerated diffusion coefficients. Consequently, Simpson and Clement [48] suggested an improved version of the SA solution in which the freshwater recharge was decreased by half. Younes and Fahs [54] and Zidane et al. [55] developed new implementations of the SA method and obtained solutions with reduced diffusion coefficients. Moreover, another limitation of the Henry’s SA solution is its unrepresentativeness of the real-world SWI as only uniform diffusion and no velocity-dependent dispersion is considered. Dispersion processes are important factors that effectively influence the mixing between fresh and salt waters [12,13,56]. A dispersive HP has been suggested by Abarca et al. [13] based on numerical simulations. As claimed by Abarca et al. [13], the major drawback of this problem is that there is no reference solution with which the results can be compared. Recently Fahs et al. [57] addressed this issue and developed a new SA solution for the dispersive HP. While all the existing SA solutions are based on the positions of the main concentration isochlors, Fahs et al. [57] derived, analytically, several metrics, characterizing the SWI as the toe length, thickness of the mixing zone, saltwater flux to the aquifer and depth of the groundwater discharge zone at the sea boundary.

This paper attempts to generalize the SA solution of the dispersive HP to more realistic configurations. Indeed, anisotropy and heterogeneity are primary characteristics of real aquifers that the SA solutions of the HP do not account for. Several former studies indicated that anisotropy can affect the saltwater penetration, the saltwater flux and submarine groundwater discharge [12,13,58]. Heterogeneity has been also found to play an important role in SWI processes [9,10,59,60,61]. Thus, the aim of this paper is to extend the SA solution of the dispersive HP, developed by Fahs et al. [57], to anisotropic and heterogeneous coastal aquifers.

The developed SA solution is based on the Fourier series method, applied to the stream function form of the governing equations. This method combines the exactness of the analytical methods with an important extent of generality in describing the geometry and boundary conditions of the numerical methods [62]. It proceeds by expending the unknowns (stream function and salt concentration) into infinite Fourier series and applying a Galerkin treatment using the Fourier modes as trial functions. We derived the stream function formulation for flow in heterogeneous and anisotropic domain. Heterogeneity introduces terms involving the space derivative of the permeability that require specific permeability–depth relationships to be addressed analytically. This was achieved by considering a stratified domain characterized via the Gardner exponential permeability function [63], which is widely used in layered aquifers, e.g., [64,65,66,67,68,69,70]. Solving the flow and transport equations in the spectral space is a challenging task, especially when large numbers of Fourier modes should be used to obtain oscillation-free solutions. Heterogeneity compounds these challenges as it involves high local gravity numbers that require a large number of Fourier modes to be represented adequately [71]. A new technique is developed here for solving the spectral flow and transport equations with a reduced number of unknowns. With this technique, we show that the HP can be solved in the spectral space with only the concentration as the primary unknown.

Various illustrative test cases are generated, and the SA solutions are compared against an in-house finite element code. The SWI metrics are evaluated analytically based on the Fourier series. These metrics represent quantitative indicators more suitable for code validation and benchmarking than the isochlors’ positions. They are also used to investigate the effects of anisotropy and stratified heterogeneity on SWI. Indeed, few studies have reported the effect of anisotropy on the isochlors’ positions [12,13,58]. The contradictory results of these studies call for further investigations [12]. In addition, the effect of heterogeneity on SWI has been considered experimentally and numerically in [10,61,72,73,74]. To the best of our knowledge, coupled effects of anisotropy and stratified heterogeneity have never been investigated. Here, we take advantage of the developed SA solution to address this gap. Moreover, while previous numerical or experimental studies investigated, visually, the uncoupled effects of anisotropy and stratification on the isochlors, this work aims to analytically provide a deeper understanding of their coupled effects on several relevant metrics characterizing SWI.

2. The Mathematical Model and Boundary Conditions

The DDF mathematical model for SWI is based upon the coupled variable-density flow and mass transport equation. Under the Oberbeck–Boussinesq approximation [75] and steady-state conditions, the flow system is given by the following continuity equation and Darcy’s law [76]:

$$\nabla \cdot q=0,$$

(1)

$$q=-K\left(\nabla h+\frac{\rho -{\rho }_0}{{\rho }_0}\nabla z\right),$$

(2)

where $q$ is the Darcy’s velocity $\left[L{T}^{-1}\right]$; ${\rho }_0$ is the freshwater density $\left[M{L}^{-3}\right]$; $K$ is the freshwater hydraulic conductivity tensor $\left[L{T}^{-1}\right]$ (diagonal with components $K_x$ and $K_z$); $h$ is the equivalent freshwater head $\left[L\right]$; $\rho$ is the density of mixture fluid $\left[M{L}^{-3}\right]$ and $z$ is the elevation $\left[L\right]$.

The salt transport processes can be described by the advection–dispersion equation:

$$\epsilon \frac{\partial c}{\partial t}+q\nabla c-\nabla \cdot \left(\epsilon D_{m}I+D\right)\nabla c=0$$

(3)

where $c$ is the dimensionless solute concentration $[-]$; $D_m$ is the molecular diffusion coefficient $\left[L^2{T}^{-1}\right]$; $\epsilon$ is the porosity $[-]$, $I$ is the identity matrix and $D$ is the dispersion tensor defined by:

$$D=\left({\alpha }_L-{\alpha }_T\right)\frac{q\times q}{\left|q\right|}+{\alpha }_T\left|q\right|I$$

(4)

where ${\alpha }_L$ and ${\alpha }_T$ are the longitudinal and transverse dispersion coefficients, respectively $\left[L\right]$.

The flow and transport equations are coupled via the linear mixture density equation:

$$\rho ={\rho }_0+\Delta \rho .c,$$

(5)

where $\Delta \rho ={\rho }_1-{\rho }_0$ is the density difference between seawater (${\rho }_1$) and freshwater (${\rho }_0$).

The geometry of the aquifer is simplified to a rectangle of length $\left(\ell \right)$ and depth $\left(d\right)$. The freshwater recharge flux per unit of width imposed on the inland side is noted as $q_d$ $\left[L^2{T}^{-1}\right]$. At the seaside, the seawater’s equivalent freshwater head is specified and constant salt concentration (c = 1) is imposed.

3. The Semi-Analytical Solution

In isotropic domains, the freshwater head (h) can be eliminated from the flow equation by applying the curl operator to Darcy’s law. For the anisotropic domain, an equivalent procedure is to differentiate the x-component (respectively, z-component) of Darcy’s law with respect to z (respectively, x) and divide it by $K_x$ (respectively, $K_z$), while accounting for spatial hydraulic conductivity variation, as the domain is heterogeneous. Subtracting the resulting equations, one from the other, gives:

$$r_K\frac{\partial q_x}{\partial z}-\frac{\partial q_z}{\partial x}=\frac{K_z}{{\rho }_0}\frac{\partial \rho }{\partial x}+r_K{q}_x\frac{1}{K_x}\frac{\partial K_x}{\partial z}-q_z\frac{1}{K_z}\frac{\partial K_z}{\partial x},$$

(6)

where $q_x$ and $q_z$ are the velocity components and $r_K=K_z/K_x$ is the hydraulic conductivity anisotropy ratio, as defined in Abarca et al. [13].

The stream function ($\psi$), defined by $q_x=\partial \psi /\partial z$ and $q_z=-\partial \psi /\partial x$, can satisfy the continuity equation. Using Equation (5) and the stream function, Equation (6) is written as:

$$r_K\frac{{\partial }^2\psi }{\partial z^2}+\frac{{\partial }^2\psi }{\partial x^2}-r_K\frac{\partial \psi }{\partial z}\frac{1}{K_x}\frac{\partial K_x}{\partial z}-\frac{\partial \psi }{\partial x}\frac{1}{K_z}\frac{\partial K_z}{\partial x}-\frac{K_z\Delta \rho }{{\rho }_0}\frac{\partial c}{\partial x}=0.$$

(7)

As in the standard HP, the following changes in variables are considered to obtain periodic boundary conditions (essential for the Fourier series method) and to derive the non-dimensional form of the flow equation (the complete mathematical development can be found in [46,55,57]):

$$X=\frac{x}{d},Z=\frac{z}{d},\Psi =\frac{\psi }{q_d}-Z,C=c-\frac{X}{\xi },$$

(8)

where $\xi =\ell /d$ is the aspect ratio.

$$r_K\frac{{\partial }^2\Psi }{\partial Z^2}+\frac{{\partial }^2\Psi }{\partial X^2}-r_K\frac{1}{K_x}\frac{\partial K_x}{\partial Z}\left(\frac{\partial \Psi }{\partial Z}+1\right)+\frac{1}{K_z}\frac{\partial K_z}{\partial X}\frac{\partial \Psi }{\partial X}-NG\left(\frac{\partial C}{\partial X}+\frac{1}{\xi }\right)=0,$$

(9)

where $NG=K_{z}d\Delta \rho /{\rho }_0{q}_d$ is the local gravity number which compares the buoyancy flux to the inland freshwater flux.

Due to heterogeneity, the stream function formulation of the flow equation incorporates the spatial derivatives of the hydraulic conductivity. Consequently, the SA solution cannot be obtained under discontinuous depth–hydraulic conductivity relationship. To overcome this limitation, we considered heterogeneities corresponding to vertical stratification. This picture is widespread and agrees with several coastal sedimentary aquifers [61,77]. We assume that heterogeneity is perfectly layering with an exponential trend of hydraulic conductivity with depth. The exponential trend of hydraulic conductivity has been suggested as a simplified model to describe stratified heterogeneity by Gardner [63]. It appears to have considerable generality and field evidence [64,65,66,67,68,69,70]. In anisotropic aquifers, the exponential depth–hydraulic conductivity model was suggested by Zlotnik et al. [68]:

$$K_x(Z)=K_{x,0}{e}^{\mathsf{\Upsilon }Z};K_z(Z)=K_{z,0}{e}^{\mathsf{\Upsilon }Z},$$

(10)

where $K_{x,0}$ and $K_{z,0}$ are the local horizontal and vertical hydraulic conductivities at the aquifer bottom surface, respectively. $\mathsf{\Upsilon }$ is the rate of hydraulic conductivity change with depth.

Usually, the hydraulic conductivity is decreasing with depth. Thus, $\mathsf{\Upsilon }$ is mostly positive. In certain cases, hydraulic conductivity can increase with depth ($\mathsf{\Upsilon }$ negative) due to hydro-fracturing. The developed SA method can be applied to both negative and positive values of $\mathsf{\Upsilon }$, but results are presented only for the common case of decreasing hydraulic conductivity with depth. The considered model (Equation (10)) assumes constant local (or small) scale anisotropy, as the anisotropy ratio is the same for all layers ($r_K=K_z/K_x=K_{z,0}/K_{x,0}$). Based on the adopted depth–hydraulic conductivity model, Equation (9) becomes:

$$r_K\frac{{\partial }^2\Psi }{\partial Z^2}+\frac{{\partial }^2\Psi }{\partial X^2}-\mathsf{\Upsilon }{r}_K\left(\frac{\partial \Psi }{\partial Z}+1\right)-N{G}_0{e}^{\mathsf{\Upsilon }Z}\left(\frac{\partial C}{\partial X}+\frac{1}{\xi }\right)=0,$$

(11)

where $N{G}_0$ ($=K_{z,0}d\Delta \rho /{\rho }_0{q}_d$) is the local gravity number at the aquifer bottom surface.

The boundary conditions for Equation (11) are given by (see Fahs et al. [57]):

$$\begin{array}{l}\frac{\partial \Psi }{\partial X}\left(0,Z\right)=\frac{\partial \Psi }{\partial X}\left(\xi ,Z\right)=C\left(0,Z\right)=C\left(\xi ,Z\right)=0,\\ \Psi \left(X,0\right)=\Psi \left(X,1\right)=\frac{\partial C}{\partial Z}\left(X,0\right)=\frac{\partial C}{\partial Z}\left(X,1\right)=0,\end{array}$$

(12)

The unknowns are expanded into infinite Fourier series that satisfy the boundary conditions:

$$\Psi ={\displaystyle \sum _{m=1}^{Nm}{\displaystyle \sum _{n=0}^{Nn}{A}_{m,n}\sin \left(m\pi Z\right)\cos \left(n\pi X/\xi \right)}},$$

(13)

$$C={\displaystyle \sum _{r=0}^{Nr}{\displaystyle \sum _{s=1}^{Ns}{B}_{r,s}\cos \left(r\pi Z\right)\sin \left(s\pi X/\xi \right)}},$$

(14)

where $A_{m,n}$ (respectively, $B_{r,s}$) are the Fourier series coefficients for the stream function (respectively, concentration). $Nm$ and $Nn$ are the truncation orders for the stream function in the X and Z directions, respectively. $Nr$ and $Ns$ are the ones for the salt concentration.

The Fourier series expansions are appropriately substituted into Equation (11). Then, a Galerkin treatment is applied, with the Fourier modes as trial functions. This leads to the following system of equations:

$$\begin{array}{l}{R}_{g,h}^F={\varpi }_h^1{\pi }^2\xi A_{g,h}\left(r_k{g}^2+\frac{h^2}{{\xi }^2}\right)+{\varpi }_h^1{r}_K\mathsf{\Upsilon }\xi {\displaystyle \sum _{m=1}^{Nm}m{A}_{m,h}{\Pi }_{g,m}}+{\pi }^{2}N{G}_{0}h{\displaystyle \sum _{r=0}^{nr}{\tilde{B}}_{r,h}{\Gamma }_{g,r}}\\ +2\pi N{G}_0{\delta }_{h,0}{\Gamma }_{g,0}+\frac{2r_K\mathsf{\Upsilon }\xi }{\pi }{\Pi }_{g,0}{\delta }_{h,0}=0\left(g=1..Nm,h=0..Nn\right),\end{array}$$

(15)

where $R_{}^F$ is the residual vector of the flow equation and ${\delta }_{i,j}$ is the Kronecker symbol.

The above procedure is similarly applied to the transport equation. This gives:

$$\begin{array}{l}{b}_m\left(\frac{{\partial }^{2}C}{\partial X^2}+\frac{{\partial }^{2}C}{\partial Z^2}\right)-\frac{\partial \Psi }{\partial Z}\frac{\partial C}{\partial X}+\frac{\partial \Psi }{\partial X}\frac{\partial C}{\partial Z}-\frac{1}{\xi }\frac{\partial \Psi }{\partial Z}-\frac{\partial C}{\partial X}-\frac{1}{\xi }+{\Delta }_{1,1}^{}\frac{{\partial }^{2}C}{\partial X^2}\\ +2{\Delta }_{1,2}^{}\frac{{\partial }^{2}C}{\partial X\partial Z}+{\Delta }_{2,2}^{}\frac{{\partial }^{2}C}{\partial Z^2}+\left(\frac{\partial C}{\partial X}+\frac{1}{\xi }\right)\left(\frac{\partial {\Delta }_{1,1}^{}}{\partial X}+\frac{\partial {\Delta }_{1,2}^{}}{\partial Z}\right)+\frac{\partial C}{\partial Z}\left(\frac{\partial {\Delta }_{1,2}^{}}{\partial X}+\frac{\partial {\Delta }_{2,2}^{}}{\partial Z}\right)=0,\end{array}$$

(16)

where $b_m=\epsilon D_m/q_d$ and ${\Delta }_{i,j}=D_{i,j}/q_d$.

The Fourier series method applied to Equation (16) gives the following system of equations (see Fahs et al. [57]):

$$\begin{array}{l}{R}_{g,h}^T={\varpi }_g^2{b}_m{\pi }^2{B}_{g,h}\left(\frac{h^2}{{\xi }^2}+g^2\right)\xi -\frac{\pi }{4}{\displaystyle \sum _{m=1}^{Nm}{\displaystyle \sum _{n=0}^{Nn}{\displaystyle \sum _{r=0}^{Nr}{\displaystyle \sum _{s=1}^{Ns}{B}_{r,s}{A}_{m,n}\left(s.m.{\eta }_{g,m,r}{\theta }_{h,n,s}-r.n.{\kappa }_{g,m,r}{\lambda }_{h,n,s}\right)}}}}\\ -g{\displaystyle \sum _{n=0}^{Nn}{\tilde{A}}_{g,n}{\Lambda }_{h,n}}-{\varpi }_g^2{\displaystyle \sum _{s=1}^{Ns}s{B}_{g,s}}{\Lambda }_{h,s}-\frac{2}{\pi }{\Lambda }_{h,0}{\delta }_{g,0}\\ -4{\displaystyle \sum _{i=1}^{Np}W{p}_i{F}_{}^{Disp}\left(X{p}_i,Z{p}_i\right)\cos \left(g\pi Z{p}_i\right)\sin \left(h\pi \frac{X{p}_i}{\xi }\right)}=0\left(g=0..Nr;h=1..Ns\right),\end{array}$$

(17)

where $R_{}^T$ is the residual vector corresponding to the transport equation. The coefficients and matrices of Equations (15) and (17) are defined in Appendix A. $Np$ is the number of integration points used to evaluate the dispersion terms. $Wp$, $Xp$ and $Zp$ are the integration weight and the coordinates of the integration points, respectively.

We should mention that, as in Fahs et al. [57], the Fourier series are used to evaluate analytically the four dimensionless metrics characterizing the SWI:

• The length of the toe ($L_{toe}$): The dimensionless distance between the seaside boundary and the point where the 50% isochlor intersects the aquifer bottom.
• The average vertical width of the mixing zone ($W_{MZ}$): Defined as the average of the vertical dimensionless distances between the 10% and 90% isochlors. The mixing zone is defined by the interval $0.3\times L_{toe}$ to $0.7\times L_{toe}$.
• The total dimensionless salt flux ($Q_s$): Represents the advective, diffusive and dispersive salt flux that enters the domain from the seaside boundary, normalized by the freshwater inland flux.
• The depth of the zone of groundwater discharge to the sea ($d_{disch}$): Equal to the distance from the aquifer top surface to the point separating the discharge zone and seawater inland flow zone at the sea boundary.

4. The New Technique for Solving the System of Equations in the Spectral Space

Solving the coupled system of flow and salt transport of the HP in the spectral space (Equations (15) and (17)) is not an easy task, particularly for cases involving a large numbers of Fourier modes. Henry [1], Segol [46] and Simpson and Clement [48] employed the Picard method and solved the system by a sequential procedure between Equations (15) and(17). This approach is computationally efficient because it splits the system into two smaller systems. However, it cannot converge for high nonlinear cases. This is why the first solutions of the HP were limited to high diffusion coefficients with relatively small number of Fourier modes. Zidane et al. [55] solved the system simultaneously using the Newton method with a numerical approximation of the Jacobian matrix. This method provides better convergence properties, but it is limited in computational efficiency, as both flow and transport equations are solved simultaneously in one large system. Younes and Fahs [54] tried to circumvent this drawback by using the Powell algorithm (a variation of Newton’s method) with the analytical Jacobian matrix. In the present work, we develop a new technique that preserves the same convergence property as the one developed by Younes and Fahs [54], but without scarifying computational efficiency of the sequential resolution (i.e., solving small systems). We show that the flow (stream function) can be analytically calculated in terms of the salt concentration, and the spectral transport equation can be solved via Newton’s method, with only the Fourier series coefficients of the salt concentration as primary unknowns.

The proof is firstly developed for homogenous aquifer. In this case, the second and fifth terms of Equation (15) drop out. This allows explicit calculation of the Fourier series coefficients of the stream function in terms of the ones for salt concentration. This gives:

$$A_{i,j}=-\frac{1}{\pi \xi }\frac{N{G}_0}{{\varpi }_j^1}\frac{\pi j{\displaystyle \sum _{r=0}^{Nr}{\tilde{B}}_{r,j}{\Gamma }_{i,r}}+2{\delta }_{j,0}{\Gamma }_{i,0}}{r_k{i}^2+j^2/{\xi }^2}.$$

(18)

Equation (18) is then substituted into Equation (17). This eliminates the coefficient $A_{g,h}$:

$$\begin{array}{l}{R}_{g,h}^T={\varpi }_g^2{b}_m{\pi }^2{B}_{g,h}\left(\frac{h^2}{{\xi }^2}+g^2\right)\xi \\ +\frac{N{G}_0}{4\xi }{\displaystyle \sum _{m=1}^{Nm}{\displaystyle \sum _{n=0}^{Nn}{\displaystyle \sum _{r=0}^{Nr}{\displaystyle \sum _{s=1}^{Ns}{B}_{r,s}\frac{1}{{\varpi }_n^1}\frac{\pi n{\displaystyle \sum _{r=0}^{Nr}{\tilde{B}}_{r,n}{\Gamma }_{m,r}}+2{\delta }_{n,0}{\Gamma }_{m,0}}{r_k{m}^2+n^2/{\xi }^2}\left(s.m.{\eta }_{g,m,r}{\theta }_{h,n,s}-r.n.{\kappa }_{g,m,r}{\lambda }_{h,n,s}\right)}}}}\\ +\frac{gN{G}_0}{\pi \xi }{\displaystyle \sum _{n=0}^{Nn}\frac{1}{{\varpi }_n^1}\frac{\pi n{\displaystyle \sum _{r=0}^{Nr}{\tilde{B}}_{r,n}{\Gamma }_{g,r}}+2{\delta }_{n,0}{\Gamma }_{g,0}}{r_k{g}^2+n^2/{\xi }^2}{\Lambda }_{h,n}}-{\varpi }_g^2{\displaystyle \sum _{s=1}^{Ns}s{B}_{g,s}}{\Lambda }_{h,s}-\frac{2}{\pi }{\Lambda }_{h,0}{\delta }_{g,0}\\ -4{\displaystyle \sum _{i=1}^{Np}W{p}_i{F}_{}^{Disp}\left(X{p}_i,Z{p}_i\right)\cos \left(g\pi Z{p}_i\right)\sin \left(h\pi \frac{X{p}_i}{\xi }\right)}=0.\end{array}$$

(19)

In the case of a heterogeneous domain, expression of $A_{i,j}$ in Equation (18) can be obtained explicitly by considering the matrix form of Equation (18):

$$\left[M\right].A=\left[N\right].B+V,$$

(20)

where $M$ is the matrix coefficients of $A_{g,h}$ (first two terms in Equation (15)), $N$ is the matrix coefficients of $B_{g,h}$ (third term in Equation (15)), $A$ and $B$ are two vectors representing, respectively, the Fourier series coefficients, $A_{g,h}$ and $B_{g,h}$, and $V$ is a constant vector (last two terms in Equation(15)).

Equation (19) represents a closed system of $\left(Nr+1\right)\times Ns$ equations with the coefficients $B_{g,h}$ as unknowns. This system can be solved alone without any coupling with the flow equation. In this work, we solved it using the nonlinear solver of the IMSL library. The Jacobian matrix is provided analytically and the term with four nested summations is simplified to three summations, as in Fahs et al. [57]. Similarly, the three resolving steps (i.e., residual vector, Jacobian matrix and linear system) are implemented in parallel on shared memory architecture.

5. Verification of the Fourier Series Solution: Stability and Comparison against Numerical Solution

A FORTRAN code has been developed to solve the final system of nonlinear equations resulting from the Fourier series method. To examine the correctness of this code, we compared it against a full numerical solution obtained using an advanced in-house model [32]. The numerical model is based on the combination of the method of lines for higher order time integration and advanced methods for spatial discretization. The mixed hybrid finite element method, which is appropriate for anisotropic and heterogeneous domains, is used to solve the flow equation. The transport equation is discretized based on the discontinuous Galerkin finite element method that can accurately handle sharp solutions with relatively coarse meshes [32]. The numerical simulations are performed with a fine mesh involving 10,000 triangular elements, obtained by subdividing 2500 squares into four equal triangles. This regular mesh is used to avoid numerical artifacts related to irregular meshes. All the simulations are developed under a transient regime for a long duration to reach the steady-state solution.

Besides validation, this section aims also to investigate the effects of heterogeneity and anisotropy on the stability of the Fourier series solution. As mentioned previously, the anisotropic layered aquifers considered in our study are described by the exponential depth–hydraulic conductivity model. Four test cases, based on the standard pure diffusive [1] and dispersive HP [13], are generated. The results are compared in terms of isochlors’ positions and SWI metrics ($L_{toe}$,$W_{MZ}$, $Q_s$ and $d_{disch}$). For an adequate comparison, we assume the same effective large-scale gravity number for all the cases. The large-scale gravity number is calculated based on the average hydraulic conductivity ($\overline{K_z}$):

$$\overline{NG}=\frac{\overline{K_z}d\Delta \rho }{{\rho }_0{q}_d}.$$

(21)

The average hydraulic conductivity (in the vertical direction) for layered aquifers with depth-dependent hydraulic conductivity is given by [68]:

$$\overline{K_z}=\frac{1}{{\displaystyle {\int }_0^1\frac{dZ}{K_z^{}(Z)}}}=\frac{K_{z,0}^{}\mathsf{\Upsilon }}{1-e^{-\mathsf{\Upsilon }}}.$$

(22)

Thus, the large-scale gravity number is calculated as follows:

$$\overline{NG}=\frac{\mathsf{\Upsilon }}{1-e^{-\mathsf{\Upsilon }}}N{G}_0^{}.$$

(23)

For all cases, the domain aspect ratio ($\xi$) is increased to 4. This is to ensure perfect horizontal streamlines at the land boundary, essential to satisfy the freshwater recharge boundary conditions with the Fourier series solution [55].

5.1. The Pure Diffusive Henry Problem

Two pure diffusive cases are considered using the same parameters as in the standard HP [1]. The first case deals with a homogenous domain ($\mathsf{\Upsilon }=0$) while, in the second one, the aquifer is assumed to be stratified ($\mathsf{\Upsilon }=1.5$). Anisotropy is acknowledged with $r_k=0.66$, as in Abarca et al. [13]. The variations of $K_x$ and $K_z$ with depth are given in Figure 2. The non-dimensional parameters used for these cases are given in

Table 1. Non-dimensional parameters used in the semi-analytical solution for the verification test cases.

Dimensionless Parameters	Value	Cases
$\overline{NG}=\frac{\overline{K_z}d\Delta \rho }{{\rho }_0{q}_d}.$	3.11	All cases
$b_m=\epsilon D_m/q_d$	0.1 5 × 10−4	Diffusive cases Dispersive cases
Non-dimensional longitudinal dispersion $b_L={\alpha }_L/d$	0 0.1	Diffusive cases Dispersive cases
Transverse to longitudinal dispersion coefficients ratio $r_{\alpha }={\alpha }_T/{\alpha }_L$	0 0.1	Diffusive cases Dispersive cases
$r_K=K_z/K_x$	0.66	All cases
The rate of heterogeneity $\mathsf{\Upsilon }$	0 1.5	Homogenous cases Heterogeneous cases

. The corresponding physical parameters used in the numerical code are summarized in

Table 2. Physical parameters used in the numerical model for the verification test cases.

Parameters	Value	Cases
$\Delta \rho$ [kg/m3]	25	All cases
${\rho }_0$ [kg/m3]	1000	All cases
$q_d$ [m2/s]	6.6 × 10−5	All cases
$d$ [m]	1	All cases
$\ell$ [m]	4	All cases
$\overline{K_z}$ [m/s]	8.213 × 10−3	All cases
$r_k$ [-]	0.66	All cases
$\epsilon$ [-]	0.35	All cases
$D_m$ [m2/s]	53.88 × 10−6 3.300 × 10−8	Diffusive cases Dispersive cases
${\alpha }_L$ [m]	0 0.1	Diffusive cases Dispersive cases
${\alpha }_T$ [m]	0 0.01	Diffusive cases Dispersive cases
$\mathsf{\Upsilon }$ [-]	0 1.5	Homogenous cases Heterogeneous cases

.

For the homogenous case, the SA solution for the anisotropic domain is sought using 1868 Fourier modes ($Nm=8$, $Nn=40$, $Nr=10$ and $Ns=140$), similar to the isotropic domain in Fahs et al. [57]. This indicates that anisotropy does not affect the stability of the Fourier series solution. For the heterogeneous case, the same number of Fourier modes leads to an unstable SA solution with some unphysical oscillations at the aquifer top (Figure 3). In fact, in this upper zone, the local permeability is five times more important than at the aquifer bottom. This high permeability leads to stronger buoyancy effects than in the homogenous case. The buoyancy forces act on a smaller space scale than convection–dispersion processes and require, consequently, a higher number of Fourier modes to be accurately simulated [71]. An oscillation-free solution has been obtained with 2968 coefficients ($Nm=8$, $Nn=40$, $Nr=10$ and $Ns=240$). Several runs confirm that heterogeneity affects only the truncation order of the concentration Fourier series in the x-direction. The SA and numerical isochlors are plotted in Figure 4. The corresponding SWI metrics are given in

Table 3. Physical parameters used in the numerical model for the verification test cases (L_toe: length of the toe, W_MZ: average vertical width of the mixing zone, Q_S: total dimensionless salt flux, d_disch: depth of the zone of groundwater discharge to the sea).

	Semi-Analytical Solution				Numerical Solution
Metrics	$L_{toe}$	$W_{MZ}$	$Q_s$	$d_{disch}$	$L_{toe}$	$W_{MZ}$	$Q_s$	$d_{disch}$
Diffusive homogenous	0.74	0.78	1.06	0.57	0.74	0.79	1.09	0.56
Diffusive heterogeneous	0.95	0.83	1.09	0.35	0.95	0.84	1.01	0.36
Dispersive homogenous	1.54	0.29	1.07	0.46	1.53	0.29	1.09	0.47
Dispersive heterogeneous	2.30	0.59	1.09	0.31	2.29	0.59	1.12	0.31

. Results show an excellent agreement between the SA and numerical solutions. This gives more confidence in the correctness of both numerical and SA codes. The agreement between the numerical and SA values of $Q_s$ is less than for the other metrics. This can be explained by the fact that $Q_s$ involves the x-derivative of the concentration that is very sensitive to the Fourier modes (in the SA solution) and the mesh size (in the numerical solution).

5.2. The Dispersive Henry Problem

For the dispersive HP, all parameters are set as in Abarca et al. [13], except for the molecular diffusion coefficient ($b_m$) which is, for convenience, assumed to be very small (5 × 10^-4) [57]. The anisotropy ratio ($r_k$) and the rate of stratification ($\mathsf{\Upsilon }$) are kept the same as for the pure diffusive cases. Non-dimensional and corresponding physical parameters are given in Table 1 and Table 2. For homogenous and heterogeneous cases, oscillation-free solutions have been obtained using 4725 ($Nm=15$, $Nn=90$, $Nr=20$ and $Ns=160$) and 6,405 ($Nm=15$, $Nn=90$, $Nr=20$ and $Ns=240$) Fourier modes, respectively. This confirms that only the concentration Fourier modes in the x-direction are mainly affected by heterogeneity. The SA and numerical isochlors are depicted in Figure 5. SWI metrics are given in Table 3. As for the homogenous case, Figure 5 and Table 3 highlight the excellent agreement between the analytical and numerical solutions. The same remark, as in the homogenous case, can be made regarding the agreement between SA and numerical values of $Q_s$. Table 3 provides quantitative indicators that can be useful for benchmarking DDF codes in realistic configurations of anisotropy and heterogeneity. First an observation from the comparison between the Figure 4 and Figure 5 indicates that the effect of heterogeneity on SWI is more important in the dispersive case. Figure 5 shows that the solution of the dispersive HP is sharper and more realistic than for the diffusive case. This sharpness explains why the dispersive test case requires more Fourier modes than the diffusive one.

It should be noted that, with the new technique developed here for solving the HP in the spectral space, the solution can be obtained with a reduced number of unknowns. For instance, in the pure diffusive homogenous case, the final system to be solved involves about 1400 coefficients, instead of 1868. For more complex cases, such as the heterogeneous dispersive case, the final system is solved with 4800 unknowns, instead of 6405.

6. Effects of Anisotropy and Heterogeneity on SWI

6.1. Effect of Anisotropy on SWI in a Homogenous Aquifer

The effect of anisotropy on SWI has been often investigated using the sharp interface model [77], which assumes that freshwater and saltwater are immiscible [78,79,80,81]. Abarca et al. [13] were the first to investigate the influence of anisotropy on the saltwater wedge using the DDF model. They demonstrated that a decrease in the anisotropy ratio ($r_K=K_z/K_x$) pushes the saltwater wedge landward. Michael et al. [58] studied numerically (based on the DDF model) the effect of anisotropy on the position of the salinity transition zone. They showed that for a lower value of the vertical hydraulic conductivity relative to a constant horizontal conductivity (i.e., lower value of $r_K=K_z/K_x$), the salinity transition zone moves seaward and can occur farther offshore. Qu et al. [12] investigated the effect of anisotropy on salinity distribution in an unconfined aquifer. Their results show that the toe recedes when the ratio of horizontal to vertical hydraulic conductivity increases (i.e., $r_K=K_z/K_x$ decreases). The results obtained by Qu et al. [12] and Michael et al. [58] are coherent but they are in contradiction with Abarca et al. [13]. This contradiction was reported in Qu et al. [12]. The authors tried to explain it by simulating their SWI problem with vertical sea–land boundaries (similar to [13]) and showed that the freshwater–seawater interface also moves seaward as $r_K$ decreases. In this section, we aim to take advantage of the developed SA solution to provide a better understanding of the effect of anisotropy on SWI. Our goal is to (i) explain the contradictory results of the previously mentioned studies and (ii) investigate the effect of anisotropy, not only on isochlors’ positions, but also on several relevant metrics characterizing the SWI.

Thus, we considered again the dispersive homogenous configuration described above in Table 2. We derived the SA solution for different cases dealing with an increasing anisotropy ratio from 0.1 to 1. Specific Fourier modes are used to obtain a stable solution for each case. Figure 6 represents the salinity distribution for three selected values of $r_K$ (0.2, 0.5 and 0.8). This figure shows that the saltwater wedge moves landward with a decrease in $r_K$. It confirms, analytically, the results of Abarca et al. [13]. The contradiction with the results of Michael et al. [58] and Qu et al. [12] can be explained based on the non-dimensional analysis presented in Section 2. In fact, Michael et al. [58] and Qu et al. [12] changed the anisotropy ratio by changing the vertical hydraulic conductivity ($K_z$) while keeping the horizontal conductivity ($K_x$) constant. This means that the isochlors’ positions are investigated at different gravity numbers as the latter is based on the vertical hydraulic conductivity ($NG=K_{z}d\Delta \rho /{\rho }_0{q}_d$). In this work, as well as in Abarca et al. [13], $K_z$ is kept constant and $r_K$ is varied based on $K_x$. This means that the isochlors’ positions are compared at a constant $NG$. In Michael et al. [58] and Qu et al. [12] the decrease in $K_z$ will decrease both $r_K$ and $NG$. The decrease in $NG$ can be alternatively seen as an increase in $q_d$ which expresses the intensification of the advection effects related to the inland freshwater recharge against the buoyancy effects (responsible for SWI). This explains why the seawater–freshwater interface is pushed back seaward when $r_K$ is decreased in Michael et al. [58] and Qu et al. [12]. In Abarca et al. [13], as well as in this work, increasing $K_x$ is equivalent to decreasing $r_K$ at a constant gravity number. At a fixed balance between the advection and buoyancy effects, the increase in horizontal permeability ($K_x$) intensifies the horizontal flow of seawater. Consequently, the saltwater wedge moves landward. To verify this conclusion, we re-evaluated the SA solution of the homogenous dispersive case by decreasing vertical hydraulic conductivity. The results are illustrated in Figure 7, which confirms that, when the decrease in $r_K$ is obtained by decreasing the vertical hydraulic conductivity, the saltwater wedge moves seaward.

Figure 8 presents the variations of $L_{toe}$ ,$W_{MZ}$, $Q_s$ and $d_{disch}$ with respect to $r_K$. Results for the homogenous case are plotted in black solid lines. Figure 8a shows that $L_{toe}$ decreases with the increase of $r_K$, which is coherent with the results in Figure 6. Figure 8b indicates that the average width of the mixing zone decreases also with $r_K$. This means that anisotropy reduces the mixing between freshwater and intruded saltwater.

Figure 8c shows that the saltwater flux to the aquifer is also controlled by the anisotropy. It decreases with the increase of $r_K$. Qu et al. [12] studied the influence of anisotropy on the submarine groundwater discharge flow. They demonstrated that the freshwater discharge is independent of anisotropy, while the seawater recirculation rate decreases as $K_x/K_z$ increases. This is in contradiction with our results which is understandable as in Qu et al. [12] $r_K$ is changed based on $K_z$ which means that the gravity number is not constant. To explain the variation of $Q_s$ against $r_K$, let us interpret the increase of $r_K$ at constant gravity number as a decrease of the horizontal hydraulic conductivity ($K_x$). The horizontal flow near the sea–aquifer interface, which mainly depends on $K_x$, is attenuated and consequently reduces the advective saltwater flux. A closer look at the saltwater flux confirms that the advective saltwater flux is mostly influenced by $r_K$. It decreases when $r_K$ is increased. From Figure 8d we can observe that the depth of the zone of groundwater discharge to the sea increases with $r_K$. This is also related to the drop of horizontal velocity near the sea boundary which pushes the inflection point of the velocity downward.

6.2. Coupled Effect of Anisotropy and Stratified Heterogeneity on SWI

Significant research has been devoted to the evaluation of the influence of heterogeneity on SWI [9,10,60,72]. Particular attention has been paid to the stratified heterogeneity using the sharp interface model [77,82,83,84]. The influence of stratification on SWI has been also investigated experimentally and numerically (using the DDF model) [61,77,85]. All these studies focused on the effect of stratified heterogeneity on the isochlors’ positions. Very limited works have investigated the influence of stratification on the thickness of the mixing zone and saltwater flux to the aquifer [61]. On the other hand, combined effects of anisotropy and heterogeneity on solute transport mechanisms have been demonstrated in Qin et al. [86] and Woumeni and Vauclin [87]. Few studies have considered the combined effects of these on SWI [72]. To the best of our knowledge, the combined effects of stratification and anisotropy on SWI have been never discussed in previous studies. Here, we used the SA solution to address these gaps. To do so, we evaluated the SA solutions for two cases dealing with increased rates of heterogeneity ($\mathsf{\Upsilon }=1.5$ and $3$) and a constant large-scale gravity number ($\overline{NG}=3.11$). We should mention that exponent $\mathsf{\Upsilon }$ is assumed to be positive. This corresponds to a decreasing hydraulic conductivity with depth which is naturally prevalent in coastal aquifers.

For each case of heterogeneity, the SA is evaluated for increasing anisotropy ratio from 0.1 to 1. The main isochlors (10%, 50% and 90%) for three selected values of $r_K$ (0.2, 0.5 and 0.8) in the case of homogenous and heterogeneous (both $\mathsf{\Upsilon }=1.5$ and $3$) aquifers are depicted in Figure 9. This figure shows that the anisotropy and rate of stratification ($\mathsf{\Upsilon }$) have complementary effects on the position of the isochlors. As for the homogenous case, increasing anisotropy pushes the saltwater wedge landward. At a constant large-scale gravity number, the increase in the heterogeneity rate leads also to a landward migration of the transition zone. This is not in agreement with the results of Lu et al. [61] (see Figure 8 in their paper). This disagreement is related to the gravity number. The results presented by Lu et al. [61] were obtained using two different large-scale gravity numbers. In Lu et al. [61], the average vertical hydraulic conductivities were ${10}^{-4}$ and ${2.7\times 10}^{-5}$ m/s for the homogenous and increasing stratified permeability, respectively. Thus, the large-scale gravity number in the homogenous case is greater in than the stratified case. This leads to more intruded saltwater wedge in the homogenous case, as a higher gravity number corresponds to a stronger buoyancy effect against freshwater recharge. To confirm this conclusion, we have compared the homogenous case of the dispersive HP (Figure 5a) to a heterogeneous case dealing with a smaller large-scale gravity number. This latter is obtained by using the same parameters as the heterogeneous dispersive case (as in Table 1), except the large-scale gravity number which is fixed at 0.84 (to get the same ratio as in Lu et al. [61]). Figure 10 illustrates the 50% isochlor for the homogenous case with $\overline{NG}=3.11$, the heterogeneous case ($\mathsf{\Upsilon }=1.5$) with the same gravity number $\overline{NG}=3.11$ and the heterogeneous case ($\mathsf{\Upsilon }=1.5$) with a reduced gravity number ($\overline{NG}=0.84$). This figure shows a significant receding of the saltwater wedge in the case of a heterogeneous aquifer with $\overline{NG}=0.84$, which is in agreement with the results of Lu et al. [61].

Figure 8a shows that $L_{Toe}$ increases with $\mathsf{\Upsilon }$ whatever the level of anisotropy. This agrees well with the results in Figure 9. However, for highly anisotropic cases (smallest value of $r_K$), $L_{Toe}$ becomes slightly sensitive to $\mathsf{\Upsilon }$. This is related to the domain length as, in this case, the saltwater wedge reaches the inland boundary. Figure 8b shows that, in general, heterogeneity intensifies the mixing between freshwater and saltwater. Except for the highest anisotropic case, heterogeneity increases the average width of the mixing zone. However, it can be observed that for a moderate hydraulic conductivity increase ($\mathsf{\Upsilon }=1.5$), the influence of heterogeneity is negligible, except in the high anisotropic cases. Figure 8c shows that the heterogeneity intensifies the saltwater flux into the aquifer. The increase in saltwater flux becomes more pronounced in highly anisotropic domains. From Figure 8d, it can be observed that the depth of the zone of groundwater discharge increases with heterogeneity, whatever the degree of anisotropy.

7. Conclusions

In this paper, the SA solution of the dispersive Henry problem was generalized to anisotropic and stratified heterogonous porous media. The solution was developed based on the Fourier–Galerkin method. This method requires a continuous depth–hydraulic conductivity relationship. The stratified heterogeneity described with Gardner’s exponential model was adopted here because it allows analytical evaluation of the Galerkin integrals and it has field evidence. The paper makes the following major contributions:

• It derives the first SA solution of SWI with the DDF model in an anisotropic and heterogeneous domain with velocity-dependent dispersion. The SA solution is useful for testing and validating DDF numerical models in realistic configurations of anisotropy and stratification. In this context, we derived, analytically (using the Fourier series), quantitative indicators (i.e., seawater intrusion metrics $L_{toe}$ ,$W_{MZ}$, $d_{disch}$ and $Q_s$) that can be effectively used for code verification.
• From a numerical point of view, an efficient technique is presented for solving the HP in the spectral space. With this technique, we showed that the governing equations in the spectral space can be solved with only the concentration as a primary unknown. The spectral velocity field can be analytically expressed in terms of concentration. This technique improves the practicality of the Henry problem’s SA solution and renders it more suitable for further studies requiring repetitive evaluations as in inverse modeling or sensitivity analysis.
• The developed SA solution is used to investigate the effects of anisotropy and stratification on SWI. This is the first time that these effects have been investigated analytically with the DDF model. In previous works, analytical studies on this issue have been limited to the sharp interface model. While in most of the existing studies, the effects of anisotropy and heterogeneity have been mainly discussed in regard to the position of the saltwater wedge, we provided here a deeper understanding of these effects on several metrics characterizing SWI.
• Taking advantage of the SA solution, we explained the contradictory results in regard to the effect of anisotropy on the position of the saltwater wedge. We showed that at a constant gravity number, the decrease in the anisotropy ratio leads to landward migration of the saltwater wedge. Contradictions observed in the previous studies are related to the way in which the anisotropy ratio is changed (whether by varying horizontal or vertical hydraulic conductivity). The SA solution shows also that anisotropy leads to a wider mixing zone and intensifies the saltwater flux to the aquifer. It leads to a shallower zone of groundwater discharge to the sea.
• The combined effects of anisotropy and stratification on SWI have been investigated. We showed that the width of the mixing zone is slightly sensitive to the rate of stratification. This sensitivity is more significant in highly anisotropic aquifers. Complementary effects of anisotropy and heterogeneity are observed on the saltwater wedge and toe position as well as on the saltwater flux, while opposite effects are observed on the depth of the groundwater discharge zone.

The developed SA solution is limited to vertical seafloor. The salt concentration along the seaside is assumed to be constant. This neglects the effect of freshwater discharge on the saltwater concentration in the discharge zone near the aquifer top surface. Future extensions of this work could include how the SA solution can be used to deal with inclined seafloor and variable saltwater concentration at the sea boundary.