Fourier series solution for an anisotropic and layered configuration of the dispersive Henry Problem

Henry Problem (HP) still plays an important role in benchmarking numerical models of seawater intrusion (SWI) as well as being applied to practical and managerial purposes. The popularity of this problem is due to having a closed-form semi-analytical (SA) solution. The early SA solutions obtained for HP were limited to extensive assumptions that restrict its application in practical works. Several further studies expended the generality of the solution by assuming lower diffusion coefficients or including velocity-dependent dispersion in the results. However, all these studies are limited to homogeneous and isotropic domains. The present work made an attempt to improve the reality of the SA solution obtained for dispersive HP by considering anisotropic and stratified heterogeneous coastal aquifers. The solution is obtained by defining Fourier series for both stream function and salt concentration, applying a Galerkin treatment using the Fourier modes as trial functions and solving the flow and the salt transport equations simultaneously in the spectral space. In order to include stratified heterogeneity, a special depth-hydraulic conductivity model is applied that can be solved analytically without significant mathematical complexity. Several examples are proposed and studied. The results show excellent agreement between the SA and numerical solutions obtained with an in-house advanced finite element code.


INTRODUCTION
Henry Problem has always been important in theoretical and practical investigation of general concepts and characteristics of seawater intrusion (SWI).It is an abstraction of SWI in a vertical cross-section of a confined coastal aquifer perpendicular to the shoreline.A schematic representation of the HP is given in Figure 1.The most important aspect of HP stems from the existence of a semi-analytical (SA) solution (Henry 1964).This solution is obtained using the Fourier series method applied to the stream function form of the flow and salt concentration 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.
In the scope of previous works concerning the SA solution of the HP, several simplifying assumptions were applied to make the calculation of SA solution feasible.While all the studies deal with uniform diffusion coefficient, recently, Fahs et al. (2016) developed a new SA solution with velocity dependent dispersion.Anisotropy and heterogeneity are primary characteristics of real aquifers that the SA solutions of the HP do not account for.Thus the aim of this work is to address this gap by extending the SA of the dispersive Henry problem to heterogeneous and anisotropic domain.In order to address the heterogeneity and to model it analytically, the Gardner exponential permeability function (Gardner 1958) was considered.This function is widely used for layered aquifers.It makes the handling of the space derivatives of permeability emerging in the process, much more applicable.The approach developed to get the SA solution is described in the next section.Then, the results of the SA solution are compared against an in-house finite element code throughout different test cases investigating the effect of anisotropy and heterogeneity.

The semi analytical solution
The mathematical model of SWI is based on coupled variable-density flow and salt transport partial differential equations.Under Oberbeck-Boussinesq approximation and steady-state conditions the flow system is as followed: (1) where q is the Darcy's velocity [LT -1 ]; the freshwater density [ML -3 ]; K the freshwater hydraulic conductivity tensor [LT -1 ] (diagonal with components K x and K z ); h the equivalent freshwater head [L]; the density of mixture fluid [ML -3 ] and z is the elevation [L].The salt transport process can be described by the advection-dispersion equation: Where c is the dimensionless solute concentration [-]; D m the molecular diffusion coefficient [L 2 T -1 ]; the porosity [-], I the identity matrix And the dispersion tensor D is defined by: where and are the longitudinal and transverse dispersion coefficient respectively [L].Flow and transport equations are coupled via the linear mixture density equation: (5) where is the difference between the freshwater ( ) and saltwater ( ) density.The geometry of the aquifer is simplified to a rectangle of length ( ) and depth ( ).The freshwater recharge flux per unit of width imposed on the inland side is noted . By differentiating Darcy's law with respect to x and z direction and eliminating pressure and applying stream function theory the system is results as followed: (6) where and and is the anisotropy coefficient that corresponds to the anisotropic nature of the domain that is being studied in this work.The non-dimensional form of the Eq(6) is defined as follows: where: is the local gravity number which compares the buoyancy flux to the inland freshwater flux ( ).
In order to address the layered heterogeneity of the aquifer, the depth-hydraulic conductivity model is defined by the following equation: where and are the hydraulic conductivity at the bottom of the aquifer in x-and zdirections respectively.is the rate of hydraulic conductivity change with depth.The non-dimensional form of the salt transport equation is presented as follows: ( ) where ⁄ and ⁄ .
The unknowns are expanded into infinite Fourier series that satisfy the boundary conditions: where A m,n (resp.B r,s ) are the Fourier series coefficients for the stream function (resp.concentration).Nm and Nn are the truncation orders for the stream function in the X-and Zdirections, respectively.Nr and Ns are the ones for salt concentration.The Fourier series expansions are then appropriately substituted in Eqs. ( 7) and (10).Then, a Galerkin treatment is applied with the Fourier modes as trial functions.This leads to the following system of equations (for flow and salt transport respectively): where (resp. ) is the residual vector for the flow equation (resp.transport equation) and is the Kronecker delta function.Np refers to the number of integration points used to evaluate the dispersion terms.
, and are the integration weight and the coordinates of the integration points, respectively.
( ) is the local gravity number at the aquifer bottom surface.
In this work a new technique is developed to solve the system of equations in the spectral space (Eqs 13 and 14).With this technique the Fourier coefficients of the stream function are analytically calculated in terms of the salt concentration and the spectral transport equation is solved the Newton's method, with only the Fourier series coefficients of the salt concentration as primary unknowns.

RESULTS
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 (Younes et al. 2009).All the simulations of the numerical model are developed under the transient regime for a long duration to reach the steady-state solution.
Two dispersive cases are considered using the same parameters as in Abarca et al. (2007).The first case deals with a homogenous domain ( ) while in the second one, the aquifer is assumed to be stratified ( ).Anisotropy is acknowledged with ( ).The non-dimensional parameters and the corresponding physical parameters used for these cases are given in Table 1.For the homogenous cases, oscillation-free solution has been obtained using 4,725 Fourier modes (Nm=15; Nn=90; Nr=20 and Ns=160).This is equal to the number of Fourier modes used in Fahs et al. [2016] for the isotropic domain.This indicates that the anisotropy does not affect the stability of the Fourier series solution.For the heterogeneous case, the same truncation modes.For the heterogeneous case, the same number of Fourier modes (4,725) leads to unstable solution with some unphysical oscillations at the aquifer top.In fact, in this upper zone, the local permeability is 5 times more important than at the aquifer bottom leading to stronger buoyancy effects.An oscillation-free solution has been obtained with 6,405 coefficients (Nm=15, Nn=90, Nr=20 and Ns=240).Several runs confirm that heterogeneity affects only the truncation order of the concentration Fourier series in the x-direction.The SA results for both cases (homogenous and heterogeneous) are depicted in Figures 2 and 3, respectively.
As mention previously, both test cases are simulated using an in-house numerical code.The physical parameters used in the numerical simulations are given in Table 1.The numerical isochlors for both test cases (homogenous and heterogeneous) are depicted in Figure 2 and 3, respectively.These Figures highlight the excellent agreement between the analytical and numerical solutions.It should be noted that, with the new technique developed in this work for solving the HP in the spectral space, the solution can be obtained with a reduced number of unknowns.For instance, in the heterogeneous case, the final system is solved using 4,800 unknowns instead of 6,405.

Figure 1 .
Figure 1.Schematic configuration of Henry Problem.