Spectral approach to seawater intrusion in heterogeneous coastal aquifers

This work presents a two-dimensional horizontal plane model of seawater intrusion into a heterogeneous coastal unconfined aquifer under the spatially random recharge condition from a stochastic point of view. The stochastic response of the coastal unconfined aquifer is attributed to uncertainties related to the heterogeneity of aquifer properties and spatial variability of the recharge. Based on the sharp-interface assumption, analytical expressions for the variances of the seawater-freshwater interface elevation and specific discharge in the longitudinal and transverse directions are developed using the perturbation approximation and spectral Fourier-Stieltjes nonstationary representations. We focus on the study of the impact of the heterogeneity of aquifer properties and recharge on these results.


Introduction
Coastal aquifers constitute an important source for freshwater supply, especially in arid and semi-arid zones (Bear, 1999).Many coastal aquifers are nowadays facing the threat of seawater intrusion.Therefore, the assessment of the extent of seawater intrusion is essential for the planning and management of groundwater resources in coastal aquifers.The extent of seawater intrusion in coastal aquifers is affected by heterogeneity and parameter uncertainty (e.g., Dagan and Zeitoun, 1998;Naji et al., 1998;Ababou and Al-Bitar, 2004;Al-Bitar and Ababou, 2005;Al-Bitar, 2007).Motivated by this, the present study attempts to quantify the uncertainty of predication of the seawater-freshwater interface taking into account heterogeneity and parameter uncertainty.
Toward this goal, the steady problem of two-dimensional flow in a heterogeneous coastal unconfined aquifer under the Correspondence to: H.-D. Yeh (hdyeh@mail.nctu.edu.tw)spatially random recharge condition is considered for the investigation.To achieve simple results of an analytical nature, it is assumed that the seawater and freshwater are separated by a sharp interface (e.g., Dagan and Zeitoun, 1998;Naji et al., 1998;Al-Bitar and Ababou, 2005).In addition, following Al-Bitar and Ababou (2005), a quadratic transform will be introduced to linearize the seawater-freshwater interface equation, which is the key to our investigation.
The introduction of spatial random recharge causes nonuniformity in the mean head gradient and results in nonstationarity in hydraulic head and velocity fields.Therefore, a nonstationary spectral approach (Li andMcLaughlin, 1991, 1995) based on Fourier-Stieltjes representations for the perturbed quantities is adopted to account for the spatial variability of nonstrationary head fields.
To the best of our knowledge, the closed-form expressions developed by the application of the nonstationary spectral approach for quantifying the uncertainty of predication of the seawater-freshwater interface with consideration of heterogeneity and parameter uncertainty have never before been presented.Note that Al-Bitar and Ababou (2005) have done similar analyses numerically in the case of zero recharge.Coastal areas often have moist climates and therefore may receive a large amount of rainfall.It is therefore of nature to characterize the interaction between freshwater and seawater with the consideration of recharge effect.It is hoped that the present results serve as estimates of the uncertainty of prediction of saltwater intrusion in helping the planning and management of the coastal groundwater resources.

Mathematical formulation of the problem
We consider the steady state flow of groundwater in a heterogeneous unconfined coastal aquifer under the spatially random recharge condition.A stochastic frame of reference is adopted to account for the spatial variability of the log Published by Copernicus Publications on behalf of the European Geosciences Union.
C.-M. Chang and H.-D. Yeh: Spectral approach to seawater intrusion hydraulic conductivity and recharge.The main assumptions used to develop analytical expressions, namely the first two statistical moments of seawater-freshwater interface elevation fluctuations, are: 1.The zone of transition from freshwater to seawater is relatively small compared to the aquifer extent and thickness.
2. The seawater wedge in the unconfined aquifer is quasistatic.
3. The freshwater flow is quasi-horizontal with a vertically hydrostatic pressure distribution.
As such, the transition zone separating the freshwater and seawater may be approximated as a sharp interface and the Ghyben-Herzberg approximation can be used to relate the seawater-freshwater interface elevation and the freshwater table.Note that as approaching the sea, the vertically hydrostatic freshwater is not possible and the assumption of essentially horizontal flow is no longer valid.As a consequence, the steady-state mass conservation of the incompressible of fluid in the freshwater zone of the unconfined aquifer takes the form (e.g., Bear, 1979) subject to the boundary conditions where X 1 and X 2 are horizontal Cartesian coordinates, η is the elevation of the seawater-freshwater interface, h f is the freshwater head, K is the hydraulic conductivity in the freshwater zone, W is the recharge rate, X 1 = X G is the location of the intersection of the interface with the substratum, and H s is the elevation of the sea level.Here and subsequently, bold face letters denote vector quantities.X L is the domain size in the X 1 -direction, between the two randomly prescribed head boundaries h f = H s and h f = H L and H L is the depth of the freshwater level above the substratum at X 1 = X L .
To simplify the analyses, it is assumed that H L and H s in boundary conditions (2) and (3), respectively, may be approximated as H L =< H L > = h L and H s =< H s >= h s by neglecting zero mean perturbations, where < > stands for the expected value operator.Note that Eq. ( 1) is the Dupuit-Boussinesq flow approximation for steady freshwater flow.A schematic representation of seawater intrusion into an unconfined coastal aquifer is illustrated in Fig. 1, where W denotes the source term in Eq. (1).Based on the Ghyben-Herzberg approximation, we can relate the elevation of the freshwater head to that of the interface as In Eq. ( 4), ρ f is the freshwater density, ρ s the saltwater density, and δ=ρ f /(ρ s −ρ f ).For the analyses below and throughout the paper, we will use ρ f D=1000 Kg/m 3 and ρ s =1025 Kg/m 3 .Thus, substituting Eq. ( 4) into Eq.( 1) leads to the following equation: where υ=δ 2 /(1+δ).Equation ( 5) can be simplified by using the φ-transform (Al-Bitar and Ababou, 2005) The idea behind it is that it is easier to compute the statistical moments of interface elevation fluctuations in the φ-field.This allows us to rewrite Eq. ( 5) as Note that Eq. ( 8) is obtained by substituting Eqs. ( 4) and (6) into Eq.( 2).The partial differential Eq. ( 7) with boundary conditions ( 8) and ( 9) constitutes the start point for developing solutions characterizing the variability of the interface elevation fluctuations in the sequel.

Stochastic equations
To develop the first two statistical moments of fluctuations in φ, the local parameters in Eq. ( 7), such as the log-hydraulic conductivity and recharge rate, are regarded as realizations of second-order stationary random fields typically represented by the sum of a mean and a small zero-mean perturbation: where < > stands for the expected value operator.The spatial persistence of the random fields, lnK and W , can be fully characterized in terms of the covariance between two locations.The spatially correlated random heterogeneity in the lnK and W parameters results in spatially correlated random fluctuations in the φ-field where ϕ(X) is the zero-mean perturbation.
A stationary process is one whose joint probability distribution is invariant when shifted in space.This implies that there is some kind of repetition of a similar structure to the variability at different locations.Thus, there is the prospect that the statistical properties of the flow field can be inferred from the analysis of one large realization.The validity of one large realization assumption requires that the correlation length of the random fields is much smaller than the domain size.
A first-order approximation of the mean boundary value problem results from substituting Eqs. ( 10)-( 12) into Eqs.( 7)-( 9) and taking the expected value as follows: subject to deterministic boundary conditions Subtracting this mean boundary value problem from Eqs. ( 7)-( 9) and dropping all products of perturbations leads to a first-order stochastic boundary value problem For convenience, the X 1 direction is aligned with the opposite direction of the mean flow so that Eqs. ( 13) and ( 16) reduce to where N = /e F and J = −∂ /∂X 1 .The first Eq.( 19) is a deterministic equation that can be solved with the boundary conditions ( 14) and ( 15).The second is a random partial differential equation for output fluctuations in φ in terms of variations in inputs lnK and W . Equations ( 19) and ( 20) provide the framework required to develop the first two moments of ϕ fluctuations in terms of the spatial covariances of the input hydraulic parameters.

Spectral solutions of the perturbation equations
The approach followed is to solve the perturbations Eq. ( 20) to fully characterize the first two moments of ϕ fluctuations.In attempting to arrive at a solution of Eq. ( 20), one must know the spatial behavior of the mean gradient J (X 1 ) in Eq. ( 20).Thus, Eq. ( 19) must be solved first in order to develop an expression for J (X 1 ).
The solution of Eq. ( 19) is given by where ξ 1 = X 1 /h s and ξ L = X L /h s .In the above, we have used the boundary conditions ( 14) and (15).Equation ( 21) holds for 0 < ξ L ≤ ξ G (where ξ G = X G /h s ).Using the above expression, we can easily obtain With Eq. ( 22), the solution of Eq. ( 20) can now be developed using the principle of superposition (e.g., Farlow, 1993) and the spectral representation (e.g., Bakr et al., 1978;Li andMcLaughlin, 1991, 1995) based on Fourier-Stieltjes representation for the perturbed quantities in wave number domain.Superposition is a useful tool in analyzing linear groundwater problems (e.g., Townley, 1995;Trefry, 1999).The principle of superposition states that a complex equation can be divided into sub-equations and the solution to the original equation is then obtained by summing the individual solution to each of the sub-equations.Based on this principle, the perturbation ϕ in Eq. ( 20) can be separated conveniently into two components: where the first term on the right-hand side (RHS) of Eq. ( 23) reflects the effect of the variation of the log-hydraulic conductivity, while the second term reflects the effect of the variation of the recharge.This leads to separate differential equations for the two components of ϕ C.-M. Chang and H.-D. Yeh: Spectral approach to seawater intrusion The solutions of Eqs. ( 24) and ( 25) can be determined using Fourier-Stieltjes representations for the perturbed quantities.By using this approach, the random perturbations in Eqs. ( 24) and ( 25) are assumed to be second-order stationary random processes and represented by the following two-dimensional wave number integral: where dZ f (R), dZ ω (R) and dZ ϕ (R) are the complex Fourier amplitudes of lnK, W and ϕ processes, respectively, and R = (R 1 , R 2 ) is the wave number vector.Here and subsequently, the integration is over two-dimensional wave number space.The mean gradient of φ is dependent of X 1 as indicated in Eq. ( 22).This space-dependent mean gradient causes the random output fluctuations in Eq. ( 24) to be nonstationary.However, these perturbed quantities in Eq. ( 24) can be presented using the nonstationary spectral representation (Li andMcLaughlin, 1991, 1995) as where ϕf (X, R) is a transfer function to be given.Thus, substituting Eqs. ( 26)-( 29) into Eqs.( 24) and ( 25) and invoking the uniqueness of the spectral representation gives the following equations The solution of ( 30) is found to be Taking the advantage of a closed-form expression, the perturbation-boundary effect on ϕ is assumed negligible (e.g., Li and McLaughlin, 1995) in obtaining the solution of Eq. ( 30).Since the perturbation approach used to quantify the perturbed fluctuations in variables is based on the assumption that the variations in parameters are relatively small, it is expected that the perturbation-boundary effect is largely limited to a small zone next to the medium boundary.Numerical simulations of two-dimensional seawater intrusion in randomly heterogeneous unconfined aquifers by Al-Bitar and Ababou (2005) illustrate that the theoretical standard deviation of φ deduced from infinite domain spectral solutions matches with the numerical standard deviation of φ approximately away from the sea and the salt-wedge tip.This implies that the assumption of negligible perturbationboundary effects is applicable, only far enough from the sea and the salt-wedge tip.

Variance of interface elevation
The expected value of the product of ϕ and its complex conjugate gives the variance of fluctuations in φ where S ff (R) is the spectrum of lnK and S ωω (R) is the spectrum of recharge rate.In Eq. ( 34), the log-hydraulic conductivity and recharge fields have been assumed to be uncorrelated.
With Eq. ( 6), the first two moments of η fluctuations can be related to those of φ by (Al-Bitar and Ababou, 2005) Thus, the equation for the variance of the interface elevation can be determined from Eqs. ( 34) and (36).

Variance of specific discharge
The vertically integrated flow equation has the form where Q i is the specific discharge in the X i direction.With the application of Eqs. ( 4) and ( 6), Eq. ( 37) becomes The first-order approximation of mean equation for the integrated specific discharge is obtained by substituting Eqs. ( 10) and ( 12) into Eq.( 38) and taking the expected value (e.g., Gelhar and Axness, 1983) Dropping products of perturbed quantities, the meanremoved form of Eq. ( 40) is Substituting Fourier-Stieltjes representations for the specific discharge perturbations, and the representations in Eq. ( 33) with its gradients into Eq.( 40) and using the uniqueness of the spectrum presentations yields where dZ Qi (R) is the complex Fourier-Stieltjes amplitude of q i .Using the representation theorem with Eqs. ( 44) and ( 45) leads to the following specific discharge spectra in the longitudinal and transverse directions, respectively, where S ff (R) and S ωω (R) are the spectra of lnK and recharge rate, respectively.Finally, the variance of specific discharge in the different directions can now be developed using Eqs.( 46) and ( 47) through the following In the remainder of this paper, we will focus on developing closed-form analytical solutions for the variances of φ fluctuations and integrated specific discharge.

Closed-form solutions
In order to evaluate the variances of φ fluctuations and specific discharge explicitly, the spectra S ff (R) and S ωω (R) must be specified.For this analysis, the random lnK perturbation field under consideration is characterized by the following spectral density function (e.g., Mizell et al., 1982;Li and McLaughlin, 1995;Li and Graham, 1999) where α f = 3π/(16λ f ) and σ 2 f and λ f are the variance and correlation scale of lnK, respectively.Note that Eq. ( 49) is found by taking the Fourier transform of an exponential correlation function of the lnK process.
In addition, the spectral density function to characterize the amplitude of recharge rate fluctuation is assumed to be described by (e.g., Li and Graham, 1999) where α ω = 3π/(16λ ω ) and σ 2 ω and λ ω are the variance and correlation scale of the amplitude of recharge rate fluctuation, respectively.

First two moments of interface fluctuations
Using Eqs. ( 49) and ( 50), integration of Eq. ( 34) over the wave number domain yields where µ f =λ f /h s and µ ω =λ ω /h s .Thus, from Eq. ( 36), it yields where is defined previously by Eq. ( 21).Equation ( 52) suggests that the variation of the interface increases linearly with the heterogeneity of the medium.
It follows from Eq. ( 4) that the head variance is related to the variance of seawater-freshwater interface elevation by The head variance is simply found from Eqs. ( 52) and ( 53) to be in the form In addition, the mean elevation of the interface can be explicitly determined after substituting Eqs. ( 21) and (51) into Eq.( 35).The behavior of the variance of seawater-freshwater interface elevation in Eq. ( 52) as a function of position for various values of µ f is illustrated in Fig. 2. It indicates that the interface elevation variance decreases with position, while it increases with µ f at a fixed location.The increase in the interface elevation variance with µ f at a fixed location can be explained by the fact that an increase in µ f produces more persistence of interface elevation fluctuations, which leads to larger deviations from the mean.It is evident from the last term on the RHS of Eq. ( 52) that the interface elevation variance always increases with µ ω at a fixed location.However, analysis of Eq. ( 52) from the field parameters shows that the increasing rate in σ 2 η with respect to the change in µ ω , is relatively small such that the impact of µ ω on the interface elevation variance is negligible compared to the impact of µ f on the same parameter.
The mean elevation of the seawater-freshwater interface resulting from Eqs. ( 21), ( 35) and ( 51) is presented graphically as a function of position in Fig. 3, which shows that larger parameter N results in reduced mean elevation of the interface.This agrees with our physical intuition that large recharge in the coastal unconfined aquifer reduces the seawater intrusion.In this paper, the position of the wedge tip is independent of the heterogeneity of the medium due to the neglect of perturbation-boundary effect.Note that the simulation results of Albitar and Ababou (2005) showed a significant increase in the mean of X G with heterogeneity. Figure 4 shows the reliability of the mean model subject to spatial heterogeneity.
The analysis leading to the results is restricted to the case of relative small hydraulic conductivity variations (weak heterogeneity) so that the second-order terms may be neglected in the flow perturbation equation.In addition, the perturbation-boundary effect is neglected in the development of analytical results.Therefore, there arises a need for comparing analytical results with numerical simulations to give some indication of the range of applicability of these results.Note that Al-Bitar and Ababou (2005) have conducted the numerical validations for stronger heterogeneity (σ 2 f = 1.6) and higher (Al-Bitar, 2007).

Variance of specific discharge
The variation of the specific discharge in the longitudinal direction is obtained from substitution of Eqs. ( 46), (49), and (50) into Eq.( 48) The variation of the specific discharge in the transverse direction following from Eq. ( 48) through the application of Eqs. ( 47), (49), and ( 50) is of the form Figure 5a and b show how the components of specific discharge variance in the longitudinal and transverse directions, respectively, vary with position for various values of N .
The variances of the longitudinal and transverse specific discharge increase with the mean recharge rate in the downstream of the flow domain, while they decrease with the increasing recharge rate in the upstream of the flow domain.This behavior is not different from the case of groundwater flow in a heterogeneous semiconfined aquifer (Lu and Zhang, 2002).It can be concluded from the examination of the last two terms on the RHS of Eqs. ( 56) or (57) that variations in µ f and µ ω do not appear to alter the longitudinal and transverse specific discharge values significantly.

A note on the mean discharge
The last term on the RHS of the mean discharge expression (39) is evaluated by using the representation theorem with Eqs. ( 26) and ( 42) From Eqs. ( 22) and ( 58), a second-order estimate of the mean discharge in Eq. ( 39) is then given by As expected, the mean specific discharge is not dependent of the form of the covariance of log-hydraulic conductivity.

A note on results as a function of the flux boundary
Using the condition, Q = −Q L at X 1 =X L (where the negative specific discharge is taken toward the coast), in Eq. ( 59), we obtain the following expression where From Eqs. ( 60), ( 21) and ( 51) can be rewritten, respectively, as The mean interface elevation as a function of the flux boundary is then obtained through applying Eqs. ( 61) and ( 62)  into Eq.( 35).In addition, with Eq. ( 60), we can easily express the variances of seawater-freshwater interface elevation in Eq. ( 52) and the longitudinal and transverse specific discharges in Eqs. ( 56) and (57), respectively, as For the case of flow through a homogeneous medium where the domain is uniformly recharged at a constant rate , we can set σ 2 f = 0, e F =K and σ 2 ω = 0.As such, coupling Eq. ( 21) with Eq. ( 60) results in Let X be a coastal directed horizontal coordinate and the origin be located distance X L from the coast.Then, the X coordinate can be related to the X 1 coordinate by X = X L −X 1 .Using this relationship, Eq. ( 66) implies which is equivalent to (9.7.7) in Bear (1972) for the homogeneous medium case.
The presence of the seawater intrusion in coastal aquifer constitutes a continuous threat to the water use.There is clear a need to view the extent of seawater intrusion and the discharge of freshwater to the sea as constraints for the selection of the optimal groundwater management scheme of a coastal aquifer.In reality, the aquifer properties such as hydraulic conductivity and the forcing functions such as recharge are uncertainty quantities.Therefore, it is important to quantify the levels of uncertainty associated with uncertainty in aquifer properties and forcing functions when designing the exploitation of fresh groundwater.Our results provide the quantification of the influence of aquifer and recharge heterogeneities on the extent of seawater intrusion and the freshwater specific discharge, which should be useful in designing constraints for the groundwater management of a coastal aquifer.

Conclusions
The goal of this study is to assess the impact of the heterogeneity of aquifer properties and recharge on the variation of the seawater-freshwater interface elevation and the specific discharge in the longitudinal and transverse directions in two-dimensional heterogeneous coastal unconfined aquifers.Toward this goal we have developed closed-form solutions for the fist two statistical moments of the interface elevation fluctuations and the specific discharge based on the perturbation approximation and spectral Fourier-Stieltjes nonstationary representations.To achieve simple results of an analytical nature, a sharp-interface assumption and the Ghyben-Herzberg approximation have been adopted.
It is found that the correlation scale of log-hydraulic conductivity plays an important role in increasing the variability of seawater-freshwater interface.On the other hand, in the presence of recharge variations in correlation scale of recharge do not appear to affect the variance of interface elevation appreciable.The mean elevation of the seawaterfreshwater interface is reduced as the recharge rate increases.
The longitudinal and transverse specific discharge variances increase with the mean recharge rate in the downstream of the flow domain, while they decrease with the increasing recharge rate in the upstream of the flow domain.In addition, they are not significantly affected by variations in correlation scales of log-hydraulic conductivity and recharge.

Fig. 2 .
Fig. 2. Dimensionless variance of seawater-freshwater interface elevation as a function of dimensionless position for various values of µ f .

Fig. 3 .
Fig. 3. Dimensionless mean elevation of seawater-freshwater interface as a function of dimensionless position for various values of N .

Fig. 5 .
Fig. 5. Dimensionless variance of the (a) longitudinal and of the (b) transverse specific discharge as a function of dimensionless position for various values of N .