Modelling Interaction Between Unconfined Groundwater and Surface Water Bodies in Aquifers Subjected to Periodic Forcing  

1. Introduction and Objective
Cyclic responses to
periodic sinusoidal forcing are well documented in the study of
groundwater systems, however there are few results in 2D or 3D.
The objective of this research is to understand the dynamics of
groundwater flow in a 2D vertical section beneath a wide shallow
lake in a regional aquifer. Full details are presented in Smith
(1999).
2. Background
Nield et al. [1994] identified 39 fundamentally different groundwater flow patterns as a function of water body geometry and the relative magnitude of aquifer flows and net recharge to the water table. Steady flow patterns were classified as being recharge, discharge or flowthrough.
Townley [1995] presented analytical solutions for the response of 1D aquifers to periodic forcing.
Two further papers by Townley [in press] present analytical solutions for a 1D aquiferlakeaquifer system, and a periodic finite element model called AQUIFEMP for computing periodic fluctuations in heads and velocities in a 2D region.
This presentation
introduces ongoing research that uses AQUIFEMP, with special
mixed boundary conditions developed from Townley [1995], to
simulate 2D flow patterns beneath a lake in a regional setting.
It generalises the 1D aquiferlakeaquifer results to 2D in the
near field of a surface water body, and allows us to determine
the extent to which the steady results obtained by Nield et
al. [1994] can be applied to dynamic systems.
3. Notation
For reasons of convenience, complex notation is used to represent sinusoidal variation. In general
f(t) = f_{s} + Re{f_{p} exp(i w t)}
where f_{s}
is the steadystate or mean value of f(t) and f_{p}
= f_{r} + i f_{i}
is the complex amplitude of fluctuation. The amplitude and phase
are the modulus f_{p} and argument
arg(f_{p}), respectively.
"Re{...}" is an operator that returns the real part of
a complex number, i = sqrt(1) and w = 2p / P is the angular frequency of period P.
4. Conceptual Model and Solution Technique
Consider a crosssection through an aquifer, with groundwater flowing from a groundwater divide towards an ocean boundary (Figure 1a).
Figure 1  Problem decomposition 
Because flow in the far field is essentially horizontal, we assume that these regions can be represented by 1D models with a transmissivity T that is spatially varying but constant in time (Figure 1b). The near field is approximated as a rectangular region, following Nield et al. [1994]. The fact that the upper surface is horizontal and fixed in time is a geometric approximation which still allows many features of the true solution to be studied.
Because the equations to be solved in the 1D and 2D regions are linear, we separate the response of the system into two parts: a steady solution, and a harmonic (purely sinusoidal) solution (Figure 1c). We use the analytical solutions of Townley [1995] to develop 3rd Type (mixed) boundary conditions at the ends of the 2D near field, such that these boundary or matching conditions encapsulate the steady and harmonic behaviour of the far field (Figure 1d). We solve both steady and harmonic problems in 2D using AQUIFEMP, and then use visualisation techniques to show pathlines and streaklines beneath the surface water body. All results are presented in terms of key nondimensional ratios (Appendix 1). The lakeaquifer geometry and boundary fluxes are depicted in Figure 2; aquifer properties include the horizontal and vertical components of hydraulic conductivity (K_{x}, K_{z}), specific yield (S_{y}), specific storativity (S_{0}), and the period of fluctuations is P.
Figure 2  Dimensional Analysis: (a) geometry, (b) boundary fluxes 
5. Results
Because heads have steady and harmonic components, specific discharge and velocity also have both steady and harmonic components, in both x and z directions. Harmonic variations in velocities result in velocity and displacement ellipses (Figure 3), with a displacement ellipse being the same shape and orientation as the velocity ellipse but scaled in magnitude by a factor of P/2p or 1/w.. When superimposed on steady flow, these result in pathlines of three kinds: smooth wavy pathlines, cuspy pathlines (with particles passing through instantaneous stagnation points) and loopy pathlines (Figure 4).
Figure 3  Velocity and displacement ellipses in harmonic flow 
Figure 4  Types of pathlines in periodic flow 
We have been keen to find out under what circumstances we would find complex cuspy or loopy flowlines in the real world. Systematic searching in the nondimensional parameter space has started to provide some indications.
Figures 5 and 6 show ellipses for lakeaquifer systems that are aquiferdriven and lakedriven, respectively. In Figure 5, the dominant fluctuation is backwards and forwards in the aquifer, with a flowthrough fluctuation through the lake. In Figure 6, the dominant fluctuation is into and out of the aquifer, with the lake alternating between recharge and discharge behaviour. In both cases, ellipses are open only for intermediate values of L_{r}^{2}S_{y}/K_{x}BP (i.e. for aquifer response times which are of the same order as the period of fluctuations), and for large values of R_{p}P/S_{y}B (i.e. when the fluctuation in water table due to recharge fluctuations is significant relative to aquifer thickness). While these results may be intuitively reasonable, they have not previously been described in quantitative terms.
Figure 5  Aquiferdriven harmonic fluctuations in groundwater flow (for enlarged view click on individual pictures) 
Figure 6  Lakedriven harmonic fluctuations in groundwater flow (for enlarged view click on individual pictures) 
Consider a steady flow pattern with groundwater discharge toward the lake. In Nield et al.'s [1994] classification this system would be described as a discharge regime of type D4 (Figure 7a). With different sinusoidal recharge to the lake and aquifer, AQUIFEMP predicts velocity and/or displacement ellipses as shown in Figure 7b. Superimposing the steady and fluctuating velocity fields and starting particles at a number of times through a period shows that particles can be completely captured by the lake (Figure 7c), or they can migrate from the lake into the lake's release zone (Figure 7d), or they can appear to be released by the lake, only to be recaptured at a later time (Figure 6e).
Figure 7  Pathline beneath a
lake with sinusoidal fluctuation in recharge (click here to see these results animated) 
6. Discussion and Conclusions
Groundwater flow patterns beneath surface water bodies during periods of dynamic changes in water levels may be much more complex than previously realised. The complex advective flow patterns beneath surface water bodies may be responsible for effects that would commonly be described as mixing or dispersion.
A sinusoidal dynamic is one of the simplest to deal with mathematically, thus AQUIFEMP provides a versatile tool for studying the details of flow patterns in vertical section. Special mixed boundary conditions at the ends of the near field, allow AQUIFEMP to account for regional behaviour.
Appendix 1  Nondimensional Ratios
Variable Ratios  
2a / B  Geometric ratio relating the size of the lake to the saturated thickness of aquifer. 
2a / L_{r}  Geometric ratio relating the size of the lake to the length of the regional aquifer. 
L_{L }/ L_{r}  Geometric ratio describing the position of the lake within the regional setting. 
R_{Ls} / R_{s}  Recharge flux ratio; relates steady state lake recharge to steady aquifer recharge. 
R_{Lp}S_{y} / R_{p}  Recharge amplitude ratio; relates the amplitude of fluctuation in lake recharge to the amplitude of fluctuation in aquifer recharge. 
L_{r}^{2}S_{y} / K_{x}BP  Non dimensional response time; describes the ability of the lakeaquifer system to respond to, or propagate, periodic forcing applied at the boundaries. Equal to the response time of the aquifer divided by the period of fluctuation. 
R_{s}L_{r}^{2 }/ K_{x}B^{2}  Non dimensional 'mound' height; a measure of the height of the regional groundwater mound. Equal to the total aquifer recharge divided by the ability of the aquifer to 'carry' it away (i.e. transmissivity K_{x}B, times a gradient B/L_{r}  allowing for a head B drop over a distance L_{r}). 
R_{p}P_{ }/ S_{y}B  Non dimensional scale ratio; describes the magnitude of the harmonic response R_{p}P_{ }/ S_{y} , relative to the depth of aquifer B. 
arg(R_{Lp} / R_{p})  Recharge phase lag; describes the phase difference between fluctuations in lake recharge and fluctuations in aquifer recharge. 
S_{y}  Aquifer specific yield. 
Constrained Ratios  
L / B = 2  Geometric ratio relating the length of the near field to the aquifer depth; chosen so that vertical flow effects induced by the lake are negligible at the near field boundaries. This implies that and are approximately uniform with depth. 
K_{x} / K_{z} = 1  Anisotropy ratio; equal to one for isotropic conditions in this work. 
S_{o}B = 0  Aquifer elastic storage coefficient; assumed negligibly small in a phreatic aquifer. 
arg(R_{p}) = 0  Phase of aquifer recharge; provides a time datum for measuring the phase of fluctuations. 
References
Nield SP. Townley LR. Barr AD. 1994. A framework for quantitative analysis of surface water  groundwater interaction: Flow geometry in a vertical section. Water Resources Research 30(8), 24612475.
Smith AJ. 1999. Periodic forcing of surface watergroundwater interaction; modelling in vertical section. PhD Thesis, Murdoch University, Western Australia, 259 pp.
Townley LR. 1995. The response of aquifers to periodic forcing. Advances in Water Resources 18(3), 125146.
Copyright © 2018 by GWSW 