Modelling and optimisation of water management in sloping coastal aquifers with seepage, extraction and recharge

TL;DR

This paper develops and validates a mathematical model for managing sloping coastal aquifers, optimizing recharge strategies, and assessing pollutant transport and seawater intrusion under various scenarios.

Contribution

It introduces a comprehensive modeling framework combining Darcy flow, pollutant transport, and seawater intrusion, validated with real data and used for optimization and scenario analysis.

Findings

Optimized recharge strategies reduce water usage compared to current practices.

Pollutants can reach the aquifer in about three years under worst-case leakage.

Double recharge rates accelerate pollutant removal from the aquifer.

Abstract

We consider the management of sloping, long and thin, coastal aquifers. We first develop a simple mathematical model, based on Darcy flow for porous media, which gives the water table height, and the flow velocity as a function of the underground seepage rate, the recharge rates and the extraction rates, neglecting sea water intrusion. We then validate the model with recent data from the Germasogeia aquifer which caters for most of the water demand in the area of Limassol, Cyprus. The data is provided over a three year period by the Cyprus Water Development Department (WDD), the governmental department managing the aquifer. Furthermore, based on our model, we subsequently develop an optimised recharge strategy and identify the optimal recharge rates for a desired extracted water volume while the water table height is maintained at the acceptable level. We study several scenarios of practical interest and we find that we can achieve considerable water savings, compared to the current empirical strategy followed by WDD. Additionally, we model the transport of pollutants in the aquifer in the case of accidental leakage, using an advection-diffusion equation. We find that in the case of an undetected and unhindered contamination (worst case scenario) the aquifer would get polluted in about three years. Also, we find that double recharge rates flush the pollutant out of the aquifer faster. Finally, to incorporate the possibility of sea water intrusion, which can render aquifers unusable, we develop a new, transient two-dimensional model of groundwater flow based on the Darcy-Brinkman equations, and determine the position of the water table and the seawater-freshwater interface for conditions of drought, moderate rainfall and flooding. The validation of the new seawater intrusion modelling approach has been carried out via comparison with a widely-accepted code.

Equations

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}}$$

$$u$$

$$w$$

$$w=\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial x}$$

$$p=p_a$$

$$w=0\mathrm{on}z=0.$$

$${\int }_0^{h}udz=Q(t)/\ell \mathrm{at}x=0,$$

$$h=H_b\mathrm{at}x=L/\cos \alpha ,$$

$$\begin{array}{rl}& x=L\hat{x},z=H\hat{z},t=\frac{L^2}{KH}\hat{t},p=p_a+\rho gH\hat{p},u=\frac{KH}{L}\hat{u},\\ & w=\frac{K{H}^2}{L^2}\hat{w},s=\frac{KH}{L^2}\hat{s},h=H\hat{h},k=\kappa \hat{k},\end{array}$$

$$ϵ\approx \alpha .$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\partial\hat{u}}{\partial\hat{x}}+\frac{\partial\hat{w}}{\partial\hat{z}}}$$

$$\hat{u}$$

$${ϵ}^2\hat{w}$$

$$\hat{p}=\hat{h}-\hat{z}.$$

$$\hat{u}=\frac{\hat{k}}{\varphi }\left(1-\frac{\partial \hat{h}}{\partial \hat{x}}\right).$$

$$\varphi \frac{\partial \hat{h}}{\partial \hat{t}}+\frac{\partial }{\partial \hat{x}}\left(\hat{k}\hat{h}\left(1-\frac{\partial \hat{h}}{\partial \hat{x}}\right)\right)=\varphi \hat{s}\hat{h}.$$

$$\hat{h}\hat{k}\left(1-\frac{\partial \hat{h}}{\partial \hat{x}}\right)=\hat{q}:\mathrm{at}\hat{x}=0,$$

$$\hat{h}=\hat{H}:\mathrm{at}\hat{x}=1,$$

$$\hat{s}=\hat{R}-\hat{E}.$$

$$\hat{R}=\sum _{i=1}^4{\hat{R}}_i(\hat{t})\delta \left(\hat{x}-{\hat{x}}_i\right).$$

$$\hat{R}=\sum _{i=1}^4{\hat{R}}_i(\hat{t})\frac{1}{\sqrt{2\pi {\sigma }_i^2}}\exp \left[-\frac{(\hat{x}-{\hat{x}}_i{)}^2}{2{\sigma }_i^2}\right],$$

$$\hat{R}=\sum _{i=1}^4{\hat{R}}_i(\hat{t})\frac{1}{2^{k/2}\Gamma (k/2)}{\left(\hat{x}-{\hat{x}}_i\right)}^{\frac{k}{2}-1}\exp \left[-\frac{(\hat{x}-{\hat{x}}_i)}{2}\right],$$

$$\hat{E}=\sum _{i=1}^{19}{\hat{E}}_i(\hat{t})\delta \left(\hat{x}-{\hat{x}}_i\right),$$

$$Err=\sum _{T_N}\sum _{D_N}{\left[\frac{\hat{h}-{\hat{h}}_f}{{\hat{h}}_f}\right]}^2$$

$$\underset{\{{\hat{R}}_i,{\hat{x}}_i,\hat{h}\}}{\min }\sum _{i=1:n_r}{\hat{R}}_i,$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\partial}{\partial\hat{x}}\left(\hat{h}\left(1-\frac{\partial\hat{h}}{\partial\hat{x}}\right)\right)=\phi\hat{s}\hat{h}}\quad$$

$$\hat{h}\left(1-\frac{\partial \hat{h}}{\partial \hat{x}}\right)=\hat{q}$$

$$\hat{h}=\hat{H}$$

$$\hat{h}-{\hat{h}}_{min}\ge 0$$

$$0\le {\hat{R}}_i\le {\hat{R}}_{max}$$

$$0\le {\hat{x}}_i\le 1$$

$$\hat{\rho }=\{\begin{array}{ll}1& \mathrm{in}\mathrm{fresh}\mathrm{water},\\ 1+\beta & \mathrm{in}\mathrm{sea}\mathrm{water},\end{array}$$

$$\hat{\mu }=\{\begin{array}{ll}1& \mathrm{in}\mathrm{fresh}\mathrm{water},\\ 1+\gamma & \mathrm{in}\mathrm{sea}\mathrm{water},\end{array}$$

$$\varphi \hat{\rho }\left(\frac{\partial \hat{u}}{\partial \hat{x}}+\frac{\partial \hat{w}}{\partial \hat{z}}\right)$$

$$\frac{{\mathrm{Fr}}^2}{ϵ}\hat{\rho }\frac{\partial \hat{u}}{\partial \hat{t}}+\hat{\mu }\varphi \hat{u}$$

$$\frac{{\mathrm{Fr}}^2}{ϵ}\hat{\rho }\frac{\partial \hat{w}}{\partial \hat{t}}+\hat{\mu }\varphi \hat{w}$$

$$\text{continuity of mass:}{\hat{\rho }}_1{\hat{\mathbf{u}}}_1$$

Full text

††journal: arXiv.org

1 Introduction

Climate change, characterised by unpredictable seasonal fluctuations, a rise in annual temperatures and a growing population make natural water conservation one of the most important environmental challenges of the 21st century, especially in locations with limited fresh water resources.

In island countries and coastal regions (for example, Cyprus [1], Texas, USA [2], Turkey [3], Morocco [4], Taiwan [5] and Greece [6]) where a significant part of the water supplies rely on aquifers, drought and over-extraction can lead to a drop in the aquifer water table. This drop increases the risk of sea water intrusion, which can render the aquifers unusable for a long period of time. In the case of drought, the water table lowers due to a decrease in the natural underground seepage flow in the aquifer. In order to maintain the water table and mitigate the risk of sea water intrusion, many coastal aquifers are replenished by pumping water back into the aquifer. Therefore, the water table height and consequently the management of aquifers requires a balance between the natural underground seepage, the extraction of water for consumption and the artificial recharge of water.

By improving aquifer management strategies, it is possible to prevent the risk of sea water contamination and, hence, reduce the dependence on alternative desalination technologies, which are energy intensive, expensive and often not feasible for small economies or in the long term.

In the past three decades, there have been many studies (both laboratory experiments and mathematical modelling) on water management in coastal aquifers. The most notable mathematical strategy for modelling the aquifer water table is the Dupuit-Forchheimer approximation [7, 8, 9, 10, 11, 12], which neglects the pressure variation in the vertical direction and considers the aquifer to be long and thin. In this approximation the Darcy equations for flow in porous media are integrated across the aquifer depth to obtain a nonlinear partial differential equation for the water table height. However, in the studies mentioned above the interplay between seepage flow rates, recharge rates and extraction rates have not been studied and in only a few cases (e.g. [10, 11]) was the developed model quantitatively compared to field data in the long term.

Moreover, there have been several reports on the controlled optimisation of the water management in the last three decades [1, 2, 3, 4, 5, 6, 13, 14, 15, 16]. Several optimisation techniques for different objectives were used. For example, Casola et al. [16] performed a time-dependent optimal control of pumping costs, Makinde et al. [3] used feedback control of the water table height with respect to a critical level, and Shamir et al. [17] used a network-based linear programming approach for multi-objective optimisation. Some of these studies modelled sea water intrusion [3, 4, 5, 6, 15] and some of them used field data from various locations, such as Cyprus [1], Texas, USA [2], Turkey [3], Morocco [4], Taiwan [5] and Greece [6].

The most comprehensive work on optimisation was by Mantoglou et al. [6], and it is very relevant to the work we present here. The model in [6] for the water table height incorporates sea water intrusion using the single potential formulation of Strack [18], and a non-linear constrained optimisation method is carried out to find the optimum pumping rates for a confined, non-sloping aquifer. However, the optimisation protocol does not incorporate the balance between underground seepage, recharge and extraction and a point source/sink formulation is assumed for the recharge/extraction process respectively. Here we develop an optimisation protocol for sloping aquifers that also incorporates the balance between seepage, recharge, and extraction. Unlike the study of Mantoglou et al. we do not incorporate a sea water intrusion model into our optimisation framework, but instead we impose a constraint that the water table must stay above a minimum critical level. This constraint serves as a proxy for avoiding sea water intrusion. Furthermore, considering that the recharge water disperses longitudinally as it enters the aquifer from the injection points at the surface, we model this with several realistic distribution functions, finding close agreement with the field data. There is nothing of probability or statistics being referred to here.

In the case of sea water intrusion, the two main modelling approaches can be divided into those in which the interface between sea water and fresh water is assumed to be sharp (i.e. no diffusion) and those in which the interface is considered diffuse. A simple and popular sharp-interface model, based on the assumption that the fresh water potential lines are vertical, is the Ghyben-Herzberg (GH) relation [19, 20]. The single potential formulation of Strack [18] is a sharp-interface model based on the GH relation (see also [6, 11]). Strack uses the method of images to find an analytical solution, assuming the water is recharged and extracted at point sources and sinks, respectively. In [21] a streamfunction-vortex formulation is used instead which models the sharp interface with a distribution of vortices [9]. Other sharp-interface models use conformal mapping techniques in the hodograph plane [22, 23, 24]. However, the above models, besides being complex, do not allow for more realistic distributions of recharge and extraction locations.

In diffuse-interface models, Henry [25] used a Boussinesq assumption and an advection-diffusion equation to model the concentration of salt in water in the aquifer. Many others have followed Henry’s diffuse-interface model [26, 27, 25, 28, 29] but these models, although more realistic than sharp-interface models, are much more computationally intensive, usually without a significant increase in accuracy.

Additional to the risk of sea water contamination, in aquifers near residential and agricultural areas there is also risk of contamination due to leaking of surface chemicals (e.g. pesticides, fertilisers, construction chemicals and human waste) [26, 30]. However, since the aquifer water quality cannot be monitored everywhere and continuously, it is important to model and study contamination transport in the aquifer and be able to make predictions about the spread and rate of transport of the contaminants. In most cases hydrodynamic dispersion is quite slow so it is quite likely that by the time contamination is detected it has affected a large part of the aquifer, and this makes modelling even more essential. Models of contaminant transport also provide an estimate of the time required to flush the contaminant naturally, before it can be reused [30, 31].

In the current study, we develop a model for the aquifer, based on the Dupuit-Forchheimer approximation, and we predict the underground water table height as a function of the recharge rates (from injection at specific points), extraction rates and the natural underground seepage rate. The natural underground seepage originates from water entering the aquifer at the bottom of the upstream dam. In Section 2 the predicted water table is validated with field measurements of the Germasogeia aquifer over a three-year period444The current work emanated from the 125th European Study Group with Industry (www.esgi-cy.org), a problem-solving workshop which took place in Cyprus, in December 2016.. The aquifer, which leads to the Mediterranean Sea, has a small inclination angle, there is a dam upstream (the source of underground seepage) and it serves most of the demand in the town of Limassol, Cyprus. The aquifer is managed by the Water Development Department (WDD), the governmental department responsible for water management in Cyprus. We explore the balance between the extraction rates, the recharge rates, and the seepage rate, and how this balance affects the water table height in the aquifer. Initially, we model the recharge and extraction sources as point sources/sinks, respectively, and then extend our formulation to line sources, where the total discharged amount is the same as in through the point injection.

In Section 3 we use a numerical optimisation method to identify the optimal recharge rates for given extraction and seepage rates. The objective of the optimisation is to minimise the total recharge whilst maintaining the water table above a critical level (so that no sea water intrusion takes place). The optimisation method not only allows for variation of the recharge rate, but also changing the number and location of the recharge sources. We can, thus, identify potential water savings, compared to the empirical protocol currently followed by WDD and make recommendations on how to improve the existing procedure by relocating the recharge sources or adding new ones.

Finally, in Section 4 we model sea water intrusion using the Darcy-Brinkman equations [32, 33, 34, 35] combined with a sharp-interface model based on the GH relation. We use this model to investigate how recharge, extraction and seepage affect sea water intrusion. Furthermore, the results of this approach are compared to those of a popular variable density groundwater tool - SUTRA, to establish reliability of the proposed model. It may be emphasized that the present approach is computationally simpler, as we only solve for the momentum equations at the interface. The mass transport equations are not solved for, unlike the case with the Henry’s approach. Generally, the available sharp interface codes are based on the physics of solving the species transport equation (creating density gradients) in the high-Peclet number case, using an optimized numerical algorithm or through the Ghyben–Herzberg approach (which inherently involves some critical assumptions) [36]. However, here we solve the fluid momentum equations for an interface, with the idea of a thin mass transfer boundary layer which is the case with high Peclet number, applicable for any practical scenario of coastal sea water intrusion. Additionally, in the case of contamination of the aquifer, we predict the contaminant spread and speed in an aquifer using an advection-diffusion equation, where the water velocity is given by the model developed in Section 2. We also make appropriate recommendations for minimising the contamination impact. This is detailed in the Appendix C of the paper.

2 Modelling the water table height in a sloping, long and thin, coastal aquifer

We consider a gravity-driven flow in a sloping, porous aquifer. We assume that the aquifer is wide enough so that width effects are negligible and the problem can be approximated in two dimensions. This is a good approximation for the Germasogeia aquifer, in the district of Limassol, Cyprus. For other aquifers, where width effects are more prominent, a full three-dimensional model may be more suitable.

This alluvial aquifer lies in the Germasogeia river valley and extends from the Germasogeia dam to the coast. It is $5.5$km long, and its depth ranges from 35m near the dam to about 55m near the sea while its width ranges from 100m to 800m (wider close to the sea)–see Fig. 1. The WDD (Limassol district) manage this aquifer and the groundwater, which is naturally filtered in the aquifer, is extracted through boreholes at several locations and used to fulfill water needs in the greater Limassol area. The WDD also artificially recharge the aquifer at four locations in order to maintain the water table at an acceptable level. Beneath the aquifer there is an impermeable bedrock, which is assumed to be flat, and the top boundary is the ground surface (see Fig. 2). There is a saturated water table level at some height below ground level. Below the water table the aquifer is fully saturated and approximately dry above it. Upstream of the aquifer a dam holds a large body of water behind a concrete barrier, and there is water seepage through the barrier with volumetric flow rate $Q$. Rainfall contribution is accounted in the model through the dam seepage. Downstream, the water table meets the sea level at height $H_{b}$ above the bedrock.

We choose a rotated coordinate system such that the $x$ direction is parallel to the bedrock level, inclined at a constant angle $\alpha$ to the horizontal, and the $z$ direction is perpendicular to the bedrock. $L$ and $H$ are the length and elevation of the aquifer, respectively, and $\tan\alpha=H/L$ (see Fig. 2). We denote by $(u,w)$ the velocity components in the $(x,z)$ directions, and by $p$ the pressure. The incompressible flow of a Newtonian fluid in a porous medium is governed by the continuity equation and the Darcy equations as follows:

[TABLE]

where $\rho$ is the density of the fluid (fresh water), $g$ is the gravitational acceleration constant, $\mu$ is the viscosity of the fluid and $\phi$ is the aquifer porosity, which we assume to be constants. We assume that the permeability of the porous medium $k(x)$ varies along the length of the aquifer. The rate of extraction (sinks) and recharge (sources) are accounted in $s(x,t)$ in the right hand side of (1). Later on, we discuss the form of $s(x,t)$ and different possible recharge and extraction distributions. The water table is denoted by $z=h(x,t)$ (blue dashed curve in Fig. 2). The kinematic and dynamic boundary conditions are

[TABLE]

where $p_{a}$ is the atmospheric pressure. At the bedrock, we have the impermeability condition

[TABLE]

The dam upstream the aquifer is a source of underground seepage and we have the following boundary condition at the dam-aquifer barrier:

[TABLE]

where $\ell$ is the typical dam width. Downstream of the aquifer the water table is at sea level

[TABLE]

which we assume constant (ignoring tidal fluctuations). We non-dimensionalise all variables using the scalings

[TABLE]

where $K=\kappa\rho g/\mu$ is the hydraulic conductivity and $\kappa$ is a typical permeability value. The inclination angle $\alpha$ is considered small ($\alpha\approx 0.01$ radians for the Germasogeia aquifer) and the aquifer is assumed long and thin such that the aspect ratio $\epsilon=\tan\alpha$ is approximately

[TABLE]

Considering (10), equations (1)$-$(3) become

[TABLE]

If we ignore terms which are of order $\mathcal{O}(\epsilon)$ and smaller, equation (13) is reduced to the classical Dupuit-Forchheimer approximation [7, 8, 9, 10, 11, 12], which we integrate (using (5)) to get

[TABLE]

Inserting (14) into (12) we find the horizontal velocity

[TABLE]

Integrating equation (11) across the water table $0\leq\hat{z}\leq\hat{h}$ and using the boundary conditions (4), (5) and (6) (and the Leibniz integral rule), we get

[TABLE]

The boundary conditions (7) and (8) become

[TABLE]

where $\hat{q}(\hat{t})=\phi QL/(KH^{2}\ell)$ and $\hat{H}=H_{b}/H$. Equation (16)$-$(18) and an appropriate initial condition, can be solved for the water table height, assuming knowledge of the source/sink terms $\hat{s}$, the seepage rate $\hat{q}$ and the permeability $\hat{k}$. We represent the contribution to $\hat{s}$ from recharge sources by $\hat{R}$ and extraction sinks by $\hat{E}$, such that

[TABLE]

Below, we will discuss various distributions for $\hat{R}$ and $\hat{E}$, and compare them when validating the model with field data from the Germasogeia aquifer, over a three-year period (Oct 2013–Nov 2016). Currently, there are four recharge locations and nineteen extraction points in operation, distributed over the length of the aquifer. The dam seepage rate varies seasonally over this period (it depends on the overground water level). The extraction points are monitored every fortnight and the recharge rate is altered by WDD every fortnight depending on the demographic demand (inferred from the extraction rates). The time series data of the dam seepage, extraction and recharge rates is presented in A (Fig. 9). The location of the four recharge and nineteen extraction points are given in the figure legends (Fig. 9). The density and viscosity of fresh water are $\rho=998~{}\,\mathrm{kg/m^{3}}$ and $\mu=0.9\times 10^{-3}~{}\mathrm{Pa.s}$, respectively. The porosity of the aquifer soil material is $\phi=0.4$, as reported by WDD. The permeability of the soil is spatially varying, but by not more than $\pm 20\%$ of its mean value, which is given in terms of the conductivity as $K=1.5\times 10^{-3}~{}\,\mathrm{m/s}=130~{}\,\mathrm{m/day}$. In our calculations, we have considered both a constant conductivity (mean value), and a spatially varying conductivity and found only minor differences in the predicted water levels. Therefore we proceeded with the assumption of constant conductivity (see Fig. 4(a)).

The recharge sources are initially considered to be point sources, on the ground. After the recharge water enters the porous medium, it percolates underground before reaching the water table. Since there are no observations of the percolation, in our model we incorporate this effect by considering three different recharge distribution patterns, sourced into the water table. We consider point sources, symmetric Gaussian distributions and asymmetric $\chi^{2}$ distributions of the water profile fed into the water table. In the case of the point source distribution,

[TABLE]

In the case of the Gaussian distribution,

[TABLE]

where $\sigma_{i}$, the standard deviation, is the parameter which controls the spatial width of the flow rate distribution, fed into the aquifer water table from the surface at single point injections (recharge). Finally, in the case of the $\chi^{2}$ distribution

[TABLE]

where $k$ is the degree of freedom of the $\chi^{2}$ distribution. Unlike the recharge process, extraction is performed using pumping wells which are drilled deep underground beneath the water table. Therefore, their effect is not distributed and we can model extraction sinks with points, as follows

[TABLE]

where the values of $\hat{E}_{i}$ are taken from the Germasogeia field data. The different distributions considered in this study are illustrated qualitatively in Fig. 3(e). The initial condition for $\hat{h}$ is chosen according to the first measurement of the water table height (Oct 2013). Since the field measurements are obtained at discrete points throughout the aquifer, we have fitted a fourth order polynomial to generate a smooth curve for the initial table height.

In Fig. 3 we display the measured water table height and the result of our model for the three different water distribution functions. Here the $\sigma_{i}$ for the Gaussian distribution and the parameter $k$ in the $\chi^{2}$-distribution are fitting parameters for the water table field data. For the $\sigma_{i}$ the values are chosen using a least-squares error (residual) minimization for the field data values and the observed model results for all the observation times recorded in the field study. The objective function in the error minimization is

[TABLE]

where $\hat{h}_{f}$ represents the observed water table height, $D_{N}$ is the number of field measurement data points and $T_{N}$ is the number of time points in the three year period of Oct 2013 - Nov 2016. The values of the optimisation variables, which are the $\sigma_{i}$ in this case, are constrained within the range of 0.0002 to 0.1. We have found $\sigma_{1}=0.01$, $\sigma_{2}=0.02$, $\sigma_{3}=0.001$, and $\sigma_{4}=0.01$ using this minimization routine. The value of $k$ in the $\chi^{2}$ distribution is obtained from a similar least-squares error minimization. However, in this case only integer values constrained in the range of 3 to 6 are allowed for $k$. We have found $k=3$ to be the optimum value in this case. There is also another reason behind the preference for choice of a $\chi^{2}$ distribution, as it asymmetric in nature and therefore realistically connects to the water flow pattern in the sloping aquifer (since water cannot travel uphill).

The comparison between the model and the data for the first 200 days is good since this data is close in time to the initial polynomial fit. Therefore, we only plot later comparisons at $285$, $570$, $855$ and $1140$ days after the initial time (Fig. 3(a)-(d), respectively). It is evident from Fig. 3 that considering the recharge sources as point sources correlates poorly with the data, compared with the Gaussian and $\chi^{2}$ distributions. We also observe that the results of the Gaussian and the $\chi^{2}$ distributions are similar. However, the Gaussian distribution is found to be better overall since agreement to data is better at larger times. The close agreement of the Gaussian distribution approach gives confidence on the accuracy of the model predictions. Note that there are some measurements of the water table height (for example in Figs. 3(c) and (d)) which are clear outliers. Regarding figures 3c and 3d, one can observe a “small cluster” of outliers at a certain water table level and distance from the dam (same at both cases), not appearing in figures 3a and 3b, which are probably due to various reasons: drought in that region, geology – impermeability of the bottom rock, changing characteristics of the aquifer soil in that region with time, which is not included in the model formulation.

In the Germasogeia aquifer, the permeability of the porous medium varies spatially along the aquifer (solid red curve in Fig. 4(a)). In Fig. 4(a) we plot the water table height $\hat{h}$ and in Fig. 4(b) we plot $\hat{u}$, calculated using our model for $\hat{t}=0.01$ (one month after the first measurement), for constant permeability (spatially averaged value) and for the varying measured permeability.

We find that the discrepancy in the water table height and in the velocity is not significant (maximum deviation of 7.4%), and hence the constant permeability assumption, in Fig. 3, is well justified. We note that in addition to incorporating a spatially varying permeability, the model is valid also for spatially varying porosity, as well as for the temporal variation of these quantities, although this is not analysed in this work.

3 Water management and optimal recharge strategy

The aquifer management issues we address here are sea water intrusion and water wastage. The types of aquifers we consider are managed with recharge points which pump water into the aquifer, thereby maintaining the water table. If the water table in the aquifer falls below a certain level, often due to excessive extraction or insufficient recharge, sea water can intrude, making it unsuitable for human consumption. The minimum water table height was determined by WDD using measurements from the Germasogeia aquifer. On the other hand, excessive recharge of the aquifer is a waste of the clean water and resources. In this study, using the recharge points (varying the recharge amount, location and number of points) to control the water table, we seek to minimise the amount of water needed to maintain the water table above the critical level. We use the optimisation procedure and data from the Germasogeia aquifer to investigate optimal recharge strategies but the methodology is valid for any sloping, sufficiently wide, long and thin aquifer.

3.1 Formulating the optimisation problem

There are a number of recharge points and we shall denote the recharge strength associated with each point by $R_{i}\,\mathrm{m^{3}/day}$. The objective of our water management is to minimise the total recharge $\sum\limits_{i}R_{i}\,\mathrm{(m^{3}/day)}$, balancing against the extraction and maintaining the minimum water level to prevent sea-water intrusion. We consider scenarios in which the strength of the recharge points is held constant for a long period of time such that the water table reaches its steady state. Hence, our optimal recharge strategy applies to the equilibrium water table, rather than an instantaneous water table. In other situations where the recharge points vary their strength much more rapidly, a time-dependent optimisation approach is more suitable.

The recharge sources are considered as Gaussian distributions, as they have been shown to have the best agreement with data (see Section 2). In addition to the input recharge rates $0\leq\hat{R}_{i}\leq\hat{R}_{max}$ (in non-dimensional terms), which are optimisation variables, we also allow the positions of the recharge points $0\leq\hat{x}_{i}\leq 1$ to vary, as this may yield further water savings. In reality, this would require installing new recharge points, but the savings in water could potentially outweigh the cost. We also consider adding new recharge points. The water table $\hat{h}$ satisfies the governing equation (16) in the steady state (no time derivative term), and boundary conditions (17) and (18). Furthermore, for the purpose of this optimisation, we restrict our attention to the case where the permeability is constant, $\hat{k}=1$. Finally, we impose a constraint on the water table height $\hat{h}\geq\hat{h}_{min}$, where $\hat{h}_{min}$ represents the critical (minimum) water level.

Therefore, to summarise, the mathematical description of our optimisation problem (in non-dimensional form) is

[TABLE]

subject to the constraints

[TABLE]

where $n_{r}$ is the total number of recharge points. We solve the optimisation problem numerically by discretising $\hat{h}$ into $N$ points. To solve equation (26) we use a second-order finite differences central scheme in space. For the boundary condition (27), we use a second-order forward scheme. We use the interior point method [37] (with the IpOpt library [38]) to solve the non-linear constrained optimisation problem. Gradients are calculated using automatic differentiation in the JuMP package [39] of the Julia programming language [40].

3.2 Application to the Germasogeia aquifer with fixed recharge locations

Having outlined the optimisation approach above, we now discuss how to use this to identify improved management protocols for the Germasogeia aquifer. For the Germasogeia aquifer, there are $n_{r}=4$ recharge points and the maximum recharge rates (due to tap capacity) are approximately $R_{max}=15000\,\mathrm{m^{3}/day}$ from each recharge point (giving a total of $60000\,\mathrm{m^{3}/day}$ maximum recharge). Firstly, we take the locations of the recharge points as their field locations, as in Section 2. Later on, we also consider variable locations. As discussed in Section 2, extraction is from 19 different locations along the aquifer. These locations lie between $x=1864\,\mathrm{m}$ and $3889\,\mathrm{m}$ (or $\hat{x}=0.33,\,0.68$). We consider each of them as point sinks, as described in Section 2. We choose extraction data from two periods which are characterised by high and low extraction rates, respectively. The first set of data is from November 2013 with a total extraction rate of $30035\,\mathrm{m^{3}/day}$, and the second is from February 2015 with a total extraction rate of $14795\,\mathrm{m^{3}/day}$. Using these two examples we develop an understanding of how to manage the aquifer when there is large or small water demand. Finally, the measured dam seepage rates are $2000-6000\,\mathrm{m^{3}/day}$ (see A). We consider rates between $500\,\mathrm{m^{3}/day}$ and $9000\,\mathrm{m^{3}/day}$ to account for and study extreme situations (e.g. drought or flooding).

Firstly, as an illustrative example, we optimise the recharge protocol for the dam seepage rate in November $2013$, which was $5076\,\mathrm{m^{3}/day}$ (as interpreted from the dam water level data, refer to Fig. 9a). Using this value, together with the extraction data measured in that period, we solve equations (25)$-$(31) to find the optimal recharge rates. In Fig. 5(a) we show the non-rotated optimised water table (solid blue line), for the optimised recharge rates, together with the minimum allowed water table (dashed red line), and the measured water table data for that period (red circles). We also plot the impermeable bedrock (dash-dotted black line), approximated with a linear slope. The optimal recharge rates are shown as black stars (refer to the values on the right hand axis). We see that the optimal recharge rates decrease with $x$, suggesting that pumping upstream is more important than downstream, and the final two points are not used at all. The total recharge amount in this example is $\sum\limits_{i=1:4}R_{i}=20699\,\mathrm{m^{3}/day}$, which is $69\%$ of the total extraction rate. Therefore, if we choose a water management strategy which recharges exactly the same amount of water as is extracted, then we would be wasting $31\%$ of the equivalent total extraction rate. We can clearly see that the optimised water table stays above the critical level, as it is supposed to, and is slightly higher than the measured water table near the start of the aquifer. This is due to large recharge rates from $R_{1}$ and $R_{2}$ which are sufficient to maintain the water table further downstream above the critical level.

The next step is to see how variation in the dam seepage and extraction rates affect the optimal recharge strategy. In Fig. 5(b) we plot the optimal recharge rates for each individual recharge point, shown with four different marker styles, whilst varying both the dam seepage (horizontal axis) and extraction rates. High extraction rates (solid lines) are taken from November $2013$ data and low extraction rates (dashed lines) are taken from February $2015$ data. For both sets of extraction data, there is a common trend in that the optimal recharge rates $R_{i}$ decrease with increasing dam seepage. This is expected as, with more seepage from the dam, there is less need to recharge in order to maintain the water table above the critical level. For small extraction rates, each recharge point has a lower optimal value than for the larger extraction rate. This is expected because, if more water is being extracted from the dam, we must account for the loss of water with extra recharge. Another observation is that there is a hierarchy of importance for the recharge points, with $R_{1}$ (first) having the largest flow rate and $R_{4}$ (last) having the smallest. Therefore, above some critical dam seepage values, recharge points $R_{2}$, $R_{3}$, and $R_{4}$ can be turned off completely. A line is drawn in the figure to indicate a dam seepage rate of $5076\,\mathrm{m^{3}/day}$, which corresponds to Fig. 5(a).

3.3 Application to the Germasogeia aquifer with variable recharge locations

If we now allow the location of the recharge points, $x_{i}$, to vary, it is interesting to see whether we can achieve more water savings. In Fig. 6 we plot the total recharge rate for when $x_{i}$ are fixed at their real values (red line with asterisks), and when $x_{i}$ are free variables (blue line with circles). We also consider the four recharge locations fixed and introduce an extra recharge point $R_{5}$ as a free variable (green line with asterisks). For both (a) high extraction rates and for (b) low extraction rates we explore the effect of variable recharge location(s). Both cases, either having low dam seepage rates and allowing variable recharge point locations, or adding a fifth point, result in a significant contribution to water savings. Water quantities of up to $9478\,\mathrm{m^{3}/day}$ and $21860\,\mathrm{m^{3}/day}$ can be saved for the low extraction regime and high extraction rates, respectively. Interestingly, the flexibility of an additional recharge point has almost the same contribution to water savings as allowing the remaining other recharge points to change location. For high seepage rates, the effect of variable recharge point locations or a new recharge point on the water saving is minimal. This is because the only active recharge point is $R_{1}$, which is located very close to the dam at $x=0$. When we allow the locations to vary freely, the locations tend to cluster towards small values of $x$ because, as we have seen, these are the most important. Thus, at large dam seepage rates there is a little advantage in relocating the existing recharge points, or installing a new one. However, for the Germasogeia aquifer, dam seepage rates are low and an optimised recharge protocol becomes more important.

In Fig. 6, we also indicate the total extraction rate (horizontal line), corresponding to both the low and high extraction cases. For high extraction rates ($30035\,\mathrm{m^{3}/day}$), if we do not constrain the recharge point locations, or add a new recharge point, then a total recharge rate equal to the total extraction rate is sufficient to maintain the water table above its critical value (except for extremely small dam seepage $\sim 500\,\mathrm{m^{3}/day}$). However, when the recharge points are fixed at their current locations, dam seepage rates below $\sim 3000\,\mathrm{m^{3}/day}$ require recharge rates much larger than the total extraction rate. For the low extraction case ($14795\,\mathrm{m^{3}/day}$), it is clear that for dam seepage rates smaller than $\sim 4500\,\mathrm{m^{3}/day}$, the total recharge must be larger than the total extraction, irrespective of the recharge point locations. Therefore, though it may seem reasonable that a total recharge rate equivalent to the total extraction rate should be sufficient to maintain the critical water level, we have seen here that this is not always the case. In particular, for small dam seepage rates it is often necessary to recharge as much as double the total extraction rate (see the left hand points in Fig. 6(b)). However, by relocating the recharge points, or adding a new point, it is possible to greatly reduce the necessary recharge rates in these situations, thereby saving water.

Although we have restricted our attention here to the optimisation of the equilibrium water table, it would be interesting to investigate how our findings would be different in the case of a time-dependent optimisation of the dynamic water table. However, we leave this to the list of future work.

4 Modelling sea water intrusion

In the case of drought or excessive extraction from the aquifer, we expect the water table to drop, thereby drawing heavier salty water from the sea into the aquifer. The interface between the fresh water and the sea water is thin because the relevant Peclet number is large, and the transport of salt is advection-dominated. The Peclet number Pe${}_{s}=KH/D_{s}$ is given in terms of the hindered diffusivity $D_{s}=\phi D_{s_{0}}$, where $D_{s_{0}}$ is the bulk diffusivity (which we assume isotropic). Therefore, if we take the bulk diffusivity of salt in water as $D_{s_{0}}\approx 10^{-9}\,\mathrm{m^{2}/s}$ [41] and, considering the characteristics of the Germasogeia aquifer $\phi=0.4$, $K=1.5\times 10^{-3}\,\mathrm{m/s}=130~{}\,\mathrm{m/day}$ (spatially averaged) and $H=80\,\mathrm{m}$, then Pe${}_{s}\approx 3\times 10^{8}$. Since, the transition region between sea water and fresh water is thin, we model it as a sharp interface between two fluids (with different densities). This reduces the computation load significantly as we don’t have to solve an additional transport equation. We introduce the non-dimensional density function

[TABLE]

and the viscosity function

[TABLE]

where $\beta$ and $\gamma$ are non-dimensional constants relating the densities and viscosities of sea water (and which can be related to salinity). We model the flow in the aquifer using the two-dimensional, unsteady Darcy-Brinkman equations [32, 33, 34, 35], which include viscous effects. Unlike the Dupuit-Forchheimer model in Section 2, this allows us to impose appropriate stress conditions on the sharp interface. As before, for the sake of simplicity, we consider constant permeability $\hat{k}=1$. In non-dimensional form, the continuity and the Darcy-Brinkman equations are

[TABLE]

where Re=$\rho KH/\mu$ is the Reynolds number and Fr${}^{2}=K^{2}H^{2}/(gL^{3})$ is the Froude number squared. Unlike the Dupuit-Forchheimer approximation in Section 2, here we do not neglect terms of order $\mathcal{O}(\epsilon)$ or higher.

We impose the boundary conditions (17)$-$(18) and additional conditions at the sharp interface between sea water and fresh water, which is denoted by $\hat{z}=\hat{h}_{i}(\hat{x},\hat{t})$, as follows:

[TABLE]

where $\mathbf{u}=(u,w)$. The index 1,2 represents the sea and fresh water characteristics and $\mathbf{I}$ represents the identity matrix. On the aquifer base $\hat{z}=0$, in addition to the impermeability boundary condition (6), we should also impose the no-slip condition

[TABLE]

As initial condition for the interface $\hat{h}_{i}$, we consider a vertical interface at the sea front ($\hat{x}=1$).

In addition, we impose the GH relationship, which relates the water table height above sea level to the depth of the sharp interface below sea level. If we denote the sea level height (in non-dimensional rotated coordinates) as $\hat{h}_{s}$, then the GH relationship is

[TABLE]

We solve the problem (34)–(40) to analyse the effect of sea water intrusion, choosing two cases to investigate. We discuss the numerical techniques used to solve this system in B. The first scenario corresponds to severe sea water intrusion. This is achieved by setting an initially very low water table, below the sea water level. Furthermore, we consider no recharge or extraction from the aquifer, and the dam seepage is set to minimal, $Q=100~{}\,\mathrm{m^{3}/day}$ . We have chosen this scenario in order to illustrate the isolated effect of having an excessively low water table.

In the second scenario, we use three different values for the dam seepage rate, $Q=500$, $2500$ and $12500\,\mathrm{m^{3}/day}$, respectively for a drought period, a moderate rainfall period, and a flooding period. We use constant extraction rates, given by the Germasogeia aquifer management data for October 2013, whereas we keep the recharge sources switched off to illustrate the effect of excessive extraction. In this case, the initial water table height is set as a polynomial fit with the observation data from October 2013 (as discussed in Section 2). In both scenarios, we assume constant permeability in the aquifer. Additionally, the sharp interface between sea water and fresh water is initially vertical, set at $\hat{x}=1$.

We plot the results of the first scenario in Fig. 7 at different times $t=0$, $1.6$, $3$ and $6.9$ years after the initial condition. The background colour represents the magnitude of the fluid velocity in the horizontal direction $\hat{u}$. We can see that there are significant regions of negative velocity. In particular, the velocity is negative at the bottom of the sharp interface, indicating that sea water is being drawn into the fresh water region. At the top of the water table the velocity is positive, indicating that the fresh water is still draining into the sea, despite the backward flow beneath. It can be observed that both the water table $\hat{h}$ in the aquifer region ($\hat{x}<1$), and the sharp interface $\hat{h}_{i}$ change position significantly over $6.9$ years.

Since sea water is denser than fresh water (with a density ratio of $1.03$), the base of the interface progresses faster than the top, as is evident from Figs. 7(c) and (d). The difference of the interface location from sea level to the aquifer base is more than 1 km after $6.9$ years.

Whilst the first scenario illustrates the effect of an excessively low water table, the second scenario focuses on the role of dam seepage in sea water intrusion. In Fig. 8a we display the position of the sharp interface between fresh water and sea water for the different dam seepage rates at two different time intervals, corresponding to 1.6 ($\hat{t}=0.2$) and 5.6 years ($\hat{t}=0.7$). In each case the sharp interface is initially at the same location and as time progresses the interface intrudes into the aquifer. As expected, with smaller dam seepage rates the intrusion is much larger. For a dam seepage rate of $500\,\mathrm{m^{3}/day}$, after around 5.6 years, sea water encroaches almost $30\%$ of the aquifer (at the base).

In order to have credibility and validate the results of the developed model, we have set up a similar model using the widely accepted saturated-unsaturated transport (SUTRA) variable-density ground-water flow package [42]. The comparison of the results using these two approaches are presented in Fig. 8b for a specific value of the dam seepage rate of 500 $m^{3}$/day. It can be seen here that the results from the SUTRA code are close to the developed model calculations (within reasonable accuracy), corresponding to $\beta=0.005$ which is equivalent to 16.7% seawater isochlor. The details of the model set up in the SUTRA package are described in B.

5 Conclusions and future work

In this work we have developed a holistic framework for the management of sloping, long and thin, coastal aquifers. We have taken as a test case the Germasogeia aquifer in Limassol, Cyprus, and used data provided over a three-year period by the Cyprus Water Development Department.

We first developed a simple two-dimensional mathematical model, based on Darcy flow for porous media, for the water table height and the water velocity, neglecting sea water intrusion (Dupuit-Forchheimer approximation). The model was validated with the data provided by WDD. Spatial permeability variation and its effect on the water table has been studied and it has been shown that small variations do not have a strong effect on the water table. Since this is the case for the Germasogeia aquifer, we have adopted in all parts of this study a constant permeability. However, in future work the model can easily be extended to study spatially varying and time-dependent permeability as may be needed in other aquifer studies.

Using the steady version of our model, we subsequently developed an optimised recharge strategy and identified the optimal recharge rates for a desired extracted water volume while the water table height was maintained at an acceptable level. We considered various practical scenarios where we varied the rates and locations of the recharge sites (while keeping the extraction sites fixed) and showed that considerable water savings can be achieved as compared to the current empirical strategy followed by WDD. In future work, it would be of further interest to develop a time-dependent optimisation protocol for a dynamically varying water table, and to additionally explore the effect of variation of the extraction locations and rates.

We then considered the case of an accidental leakage and predicted the transport of pollutants in the aquifer for various scenarios of practical interest. We used an advection-diffusion equation for the pollutant and assumed that its velocity is given by the water velocity, as this was determined in Section 2. Three scenarios were considered for the duration of the leakage: non-stop leakage, and “pulse” leakage lasting one day or 30 days. In the case of non-stop, unhindered and undetected leakage (worst case scenario) we found that the aquifer would be completely polluted in approximately three years. In the two cases of pulse leakage, the impact is less severe. Exploring also variation of the recharge rates we found that contamination can be flushed out faster as the recharge rates increase. As the Peclet number is very large (diffusion effects are very small) the contaminant propagation front is almost vertical and this fact suggests that a monitoring device would be expected to detect sudden changes in contaminant concentration levels.

Finally, we developed a time-dependent, two-dimensional model of saturated-unsaturated groundwater flow based on the Darcy-Brinkman equations, and we showed that the problem can be case as a two-fluids problem with a sharp interface since diffusion effects are very small. This framework allowed us to easily study sea water intrusion by tracking the position of the water table and of the seawater-freshwater interface. We applied the model in four scenarios of interest; initially a very low water table with no extraction or recharge and then for three potential cases-drought, moderate rainfall and flooding. We found that even in the case of drought (dam seepage rate of $500\,\mathrm{m^{3}/day}$) we expect sea water to encroach almost $30\%$ of the aquifer (at the base) in less than 6 years. The analysis done here offers an alternative, less computationally expensive methodology to solving the fully coupled density-concentration diffusion problem. This part of the work is validated using the popular SUTRA code for the variable density groundwater flow and the comparison of the results is also shown here. Comparison of the computational resource requirement for this new approach is a subject of further investigation, which is outside the scope of the present study. However, the current work is aimed towards providing a foundation in the development of new alternative efficient technique.

In the present study the width of the aquifer has been assumed infinite (no edge effects) and hence we could treat the problem as two-dimensional. As we have seen this is a good approximation for the Germasogeia aquifer for a significant part of the aquifer but we expect a closer agreement with data if we take the width variation into account. It is also expected that also in other aquifers we would not be able to neglect the effects of width variation.

Acknowledgements

The authors would like to acknowledge the 125th European Study Group with Industry (ESGI125, www.esgi-cy.org) which was held in Limassol, Cyprus, 5-9 December 2016. Dr K. Kaouri (Chair of the ESGI125 Organising COmmittee) would especially like to thank the Mathematics for Industry Network (MI-NET, COST Action TD1409), KPMG and the Cyprus University of Technology for supporting this workshop. The authors are also thankful to the Cyprus Water Development Department (WDD), Limassol District Office, for proposing the problem, fruitfully collaborating with the Study Group participants and for financial support. The authors also acknowledge support from WDD that aided them to continue the work after the Study Group. The authors are also thankful to all the members of the ESGI125 team for stimulating discussions: Prof D. Papageorgiou (Imperial College London, UK), Dr Hilary Ockendon (University of Oxford, UK), Dr C. Nikolopoulos (University of the Aegean, Greece), Dr T. Ivanov (Sofia University, Bulgaria), Prof. P. Frolkovic (Slovak University of Technology, Slovakia), Mr A.N. Riseth (University of Oxford, UK), Prof. A. Lacey (Heriot-Watt University, UK) and Dr V. Rottschäfer (Leiden University, Netherlands). We would also like to thank Prof G. Sander (University of Loughborough/University of Oxford) for helpful discussions.

G. Benham acknowledges the partial financial support from WDD and the University of Oxford in the second phase of the work, after the Study Group. Dr R. Mondal acknowledges financial support from MI-NET through a Short Term Scientific Mission (STSM) grant. Dr S. Mondal also acknowledges the partial support from WDD and the Royal Society towards travel and subsistence.

Appendix A Time-series data of the recharge, extraction and dam seepage

Appendix B Numerical details for sea water intrusion

The numerical computations are carried out using the finite-element-based commercial package COMSOL v5.2 [43, 44]. The discretization of the variational form of the governing equations follows the standard Galerkin formalism [45]. To save computational cost, the Jacobian matrix is updated once in several iterations. Within each Newton iteration, the sparse linear system is solved by preconditioned Krylov methods [46] such as the generalized minimum residual (GMRES) method [47, 48] and the bi-conjugate gradient stabilized (BCGSTAB) method [49]. The segregated solver [43] is used with a relative tolerance of 0.001. Two segregated steps decoupling the hydrodynamics and the director equations are set up.

The time-dependent solver is utilized with a relative tolerance of $10^{-4}$ along with the moving mesh method based on Arbitrary Lagrangian–Eulerian (ALE) formulation [45, 50, 51], which enables the moving boundaries without the need for mesh movement. This allows for mesh deformation with the translation of the physical boundary. Smoothing the deformation (because of the movement of the interior nodes) throughout the domain is achieved by solving the hyperelastic [52, 53] smoothing partial differential equation with the boundary conditions specified. Second-order shape elements are used to mesh the geometry in the case of free deformation with particle movement. The spatial derivatives of the dependent variables are now transformed with respect to the new transformed coordinates. When the mesh deformation produces a low-element quality (defined as the ratio of the inscribed to circumscribed circles’ radii for each simplex element), where the minimum element quality is less than 0.3 (it is 0 for the degenerated element and 1 for the best possible element), remeshing of the entire domain is done and the simulation progresses.

Regarding the SUTRA simulation, the present situation is similar to the Henry problem. Initially a 2D rectangular box was created, representing the calculation domain, as shown in Fig. 2, except we extend the length of the aquifer to be $2L$ so as to account for the downstream sea region, as described in Section 4. The fresh water is on the left side and the sea water is through the right boundary, similar to the configuration described in Fig. 2. We have used an irregular mesh generation method [54] with an average mesh size of 5m, as illustrated in Fig. 10. The transport scheme selected is “solute using pressure” and the flow conditions are considered to be saturated. The boundary conditions as presented in the current work are incorporated in this simulation. For the physical properties of the aquifer, we have considered the hydraulic conductivity of the domain to be constant at 130 m/day and the porosity was set to 0.4. The ratio of the density of sea water to freshwater is considered to be 1.03 ($\beta=0.03$). For the operating conditions we have used constant extraction rates equivalent to the values of October 2013 at all the 19 locations, and the recharge rates are switched off to clearly illustrate the effect of the sea water intrusion. The dam seepage rate is assumed to be constant at the rate of 500 $m^{3}$/day. The model is set up using the ModelMuse v3.10 GUI package [55] which is the subsequently fed in to the SUTRA code. First the simulation is run for steady state to generate the initial data and then it is “restarted” using the result (of the steady state run) as an initial condition for the transient calculations over a 6 year period. The contour lines for the concentration are for $\beta=0.005$ (a measure of the density ratio) plotted at different times, as shown in Fig. 8b.

Appendix C Contaminant transport

In the event of an accidental spillage or seepage of chemical contamination in an aquifer, it is important to understand and forecast the impact of contamination. For simplicity, we restrict our attention to the case where the contamination is uniform in the vertical direction $\hat{z}$. Therefore, the contaminant undergoes both horizontal advection by the fluid velocity $\hat{u}$ (see equation (15)), and molecular diffusion. In practice, usually the contamination volume is small relative to the volume of the aquifer, so that it does not affect the flow field significantly, and we proceed with this assumption. In this way, the transport of contaminant (species) concentration $c$ (in the dilute range) can be modelled with an one-dimensional advection-diffusion equation. Here we define the non-dimensional concentration $\hat{c}=c/c_{0}$, where $c_{0}$ is a typical concentration value, such that

[TABLE]

where the Peclet number Pe${}_{c}=KH/D_{c}$ is given in terms of the hindered diffusivity of the contaminant in the porous medium, $D_{c}\approx\phi D_{c0}$ (where $D_{c0}$ is the bulk diffusivity). Here we also do not consider the hydrodynamic dispersion effects, so effectively $D_{c}\not\approx f(u)$. We also assume that the spillage occurs at an isolated location for a finite time $\hat{x}_{c}$ for $0\leq\hat{t}=[\leq\hat{t}_{c}]$ (“pulse” spillage). In equation (41), $G$ represents the spillage and can be approximated by a Dirac-delta function $\delta$, such that

[TABLE]

We consider a very diffusive solute with $D_{c_{0}}=2\times 10^{-9}\,\mathrm{m^{2}/s}$, which leads to a Peclet number Pe${}_{c}\approx 2\times 10^{8}$. Even with a very diffusive solute, we still find Pe${}_{c}\gg 1$, and hence the contaminant transport is advection dominated.

Therefore, it is insightful to analyse the worst case scenario, when the spill happens at the dam ($\hat{x}_{c}=0$). The boundary condition at the far end of the aquifer ($\hat{x}=1$) is the outlet condition that there is no change in the spatial concentration there, that is

[TABLE]

The boundary condition at the dam $(\hat{x}=0)$ is also that no concentration flux exists there, since we do not expect the contaminant to diffuse across the dam boundary. Furthermore, we assume that at the time of the spillage, there is no contaminant anywhere in the aquifer, so we have the initial condition

[TABLE]

We solve equations (41)$-$(44) for the concentration of a chemical contaminant $\hat{c}$, given the spillage function (42), where the advection velocity $\hat{u}$ is calculated using the Dupuit-Forchheimer model (Equations (16)$-$(18)) with recharge, extraction and dam seepage rates taken from the Germasogeia aquifer data with Gaussian distributions for the recharge sources (as described in Section 2). The permeability of the aquifer is considered constant, $\hat{k}=1$. Since the contaminant travels downstream due to fluid advection, the spillage affects more if it occurs upstream.

We also consider different spillage durations $\hat{t}_{c}$ . In the first case we choose a continuous (uninterrupted) spillage ($\hat{t}_{c}\rightarrow\infty$), in the second case a spillage time at $\hat{t}_{c}=3.4\times 10^{-4}$ (1 day), and in the third case at $\hat{t}_{c}=0.01$ (30 days). The second and third cases correspond to situations in which the spillage is detected and contained successfully. In Figs. 11 and 12, we show how the contaminant is transported across the aquifer as time progresses for the three cases considered.

In the first case (Fig. 11) we see that the contaminant propagation front is almost vertical, which is explained by the fact that the Peclet number is very large. This suggests that a monitoring device would detect a very sudden change in contaminant concentration levels (zero to unity within a couple of weeks as seen in Fig. 11b). We can see that if the contamination spillage is undetected and unhindered, the aquifer gets completely polluted within 3 years ($\hat{t}\approx 0.37$).

In the second and third cases, where the contaminant spillage is detected, the spillage source is stopped in a finite time. The contaminant, which enters the aquifer, gets advected and eventually the concentration levels throughout the aquifer fall to an acceptable level $(\varepsilon_{c})$. The concentration profile at different times is displayed for the second case in Fig. 12(a) and for the third case in Fig. 12(b). In each case, the contaminant transportation behaviour is similar, with the concentration isolated to a pulse which spreads out over time. In the case of the longer spillage, the magnitude of the concentration is higher overall. The results also provide information about the time dependent position of the contaminant and one can infer the region of safety ($\hat{c}<\varepsilon_{c}$).

To illustrate the effect of recharge on contaminant transport, we display the concentration profiles for double recharge rates as dotted lines in Fig. 12. We clearly see that by increasing recharge, the contaminant is flushed more quickly. Furthermore, in the case of the briefer spillage, increasing the recharge has a greater flushing effect than for the longer spillage (note the profiles for $\hat{t}=0.0125$ in Figs. 12a and b).