Hypoaigic influences on groundwater flux to a seasonally saline river

Journal of Hydrology (2007) 335, 330– 353

available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/jhydrol

Hypoaigic influences on groundwater flux to a seasonally saline river M.G. Trefry *, T.J.A. Svensson 1, G.B. Davis CSIRO Land and Water, Private Bag No. 5, Wembley, WA 6913, Australia Received 9 August 2006; received in revised form 29 November 2006; accepted 1 December 2006

KEYWORDS Groundwater–surface water interaction; Saline; Seasonal; Seepage; Convection

Hypoaigic zones are aquifer volumes close to and beneath the shores of saline surface water bodies, and are characterized by the presence of time-dependent natural convection and chemical stratification. When transient and cyclic processes are involved there is significant potential for complex flow and reaction in the near-shore aquifer, presenting a unique challenge to pollutant risk assessment methodologies. This work considers the nature of some hypoaigic processes generated by the seasonally saline Canning River of Western Australia near a site contaminated by petroleum hydrocarbons. A dissolved hydrocarbon plume migrates within the shallow superficial aquifer to the nearby bank of the Canning River. Beneath the river bank a zone of complex fluid mixing is established by seasonal and tidal influences. Understanding this complexity and the subsequent ramifications for local biogeochemical conditions is critical to inferring the potential for degradation of advecting contaminants. A range of modelling approaches throws light on the overall topographic controls of discharge to the river, on the saline convection processes operating under the river bank, on the potential for fluid mixing, and on the various important time scales in the system. Saline distributions simulated within the aquifer hypoaigic zone are in at least qualitative agreement with previous field measurements at the site and are strongly affected by seasonal influences. Groundwater seepage velocities at the shoreline are found to be positively correlated with river salinity. Calculations of fluid age distributions throughout the system show sensitivity to dispersivity values; however, maximum fluid ages under the river appear to be diffusion limited to a few decades. The saline convection cell in the aquifer defines a zone of strong dispersive dilution of aged (many decades) deep aquifer fluids with relatively young (several months) riverine fluids. Seasonal recharge and river salinity cycles induce regular perturbations to the convection cell, yielding intra-annual variations of 50% in seepage velocity and almost 30% in wedge penetration distance at the plume location. ª 2006 Elsevier B.V. All rights reserved.

Summary

* Corresponding author. Tel.: +61 8 9333 6286; fax: +61 8 9333 6211. E-mail address: [email protected] (M.G. Trefry). 1 Present address: GEOSIGMA AB, Vattholmav. 8, Box 894, 751 08 Uppsala, Sweden. 0022-1694/$ - see front matter ª 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2006.12.001

Hypoaigic influences on groundwater flux to a seasonally saline river

331

Nomenclature $ vector gradient operator D difference operator aL, aT longitudinal and transverse dispersivities (m)  a density factor (–) h aquifer porosity (–) q, qmin, qmax density parameters of the saline fluids (g cm3) A unit thickness for a vertical section model (m) As van Genuchten capillary saturation scaling coefficient (m1) a,b,c empirical dispersion functional parameters (–) AHD Australian height datum AMG Australian map grid LNAPL light non-aqueous phase liquid B aquifer thickness (m) BTEX benzene, toluene, ethylbenzene, xylene (hydrocarbons) Cmax maximum saline concentration (g L1) D aquifer diffusivity (m2 d1) Dmol molecular (free) diffusion coefficient in water (m2 d1) fW empirical dispersion functional (–) h hydraulic head (m) K aquifer hydraulic conductivity (m d1) Kx, Ky, Kz vector components of hydraulic conductivity (m d1)

Introduction Recent advances in the study of physical and biological interactions in near-shore subsurface zones surrounding surface water bodies, the so-called hyporheic zones (from the ancient Greek words hypo meaning ‘‘beneath’’ and rheos meaning ‘‘stream’’), have stimulated interest in the dynamics, fates and risks of hyporheic contamination. Particular issues surround the natural biogeochemical abilities of near-shore zones to endure and even remediate groundwater pollutants as they follow local flow vectors. General concepts in the science of hyporheic environments have been established, mainly within the context of fresh surface waters, e.g. upland rivers and streams (Stanford and Ward, 1988), linking hyporheic controls to sediment respiration rates (Naegeli and Uehlinger, 1997), bacterial attenuation of nutrients (e.g. Clavero et al., 1999), riparian vegetation regimes (Lambs, 2004) and, more broadly, to ecological health (Young and Huryn, 1996). Seasonal reversals of seepage have been shown to induce significant shifts in freshwater aquifer hydrochemistry (Helena et al., 2000), and complex transient relationships between aquifer and stream chemistry have been reported in contaminated mining locales (Winde and van der Walt, 2004). In this paper, we draw a distinction between hyporheic zone processes involving fresh waters, and the very different sub-shore processes involving density contrasts between surface waters and groundwaters commonly found along

K1 m

lunisolar diurnal constituent, period 23.93 h van Genuchten capillary saturation exponent quantity (–) M2 principal lunar semidiurnal constituent, period 12.42 h n = 1/(1  m) van Genuchten capillary saturation exponent (–) O1 lunar diurnal constituent, period 25.82 h PEST parameter estimation module, implemented in Feflow PSC percentage recirculation ratio (%) q Darcy groundwater flux (m2 d1) Ra equivalent cellular Rayleigh number (–) Ra 0 modified Rayleigh number (–) S aquifer storativity (–) Sr soil residual water saturation (–) S2 principal solar semidiurnal constituent, period 12.00 h si arc length (m) T aquifer transmissivity (m2 d1) t time (s, days or years) v groundwater (pore fluid) speed (m d1) Vamb ambient regional groundwater horizontal speed (m d1) W ratio of aquifer depth to transverse dispersivity (–) x, y, z spatial coordinates (m)

coastlines and estuaries. We describe the latter convective interfacial environments as hypoaigic, from the ancient Greek word aigialos, meaning ‘‘of the sea-shore’’ (see also the term ‘subterranean estuary’’ introduced by Moore (1999) and the earlier discussion by Fetter (1994, p. 370)). Key hypoaigic concepts are the mixing of chemically disparate and variable density fluids (e.g. sourced from groundwater and surface waters), the effect of periodic fluid motions (e.g. induced by tidal or seasonal forcings) and the concentration and convective cycling of nutrients, contaminants and other variables in the groundwater seepage zone. Hypoaigic environments encompass the same rich array of physical, chemical and biological processes characteristic of hyporheic zones, but with the added complexity of transient and density-driven flow. Dense saline surface waters penetrate into the aquifer, displacing less dense aquifer fluids at the bottom of the aquifer, forming a saline wedge. Regional groundwater flux disperses with the wedge, creating saline concentration and density gradients which, in turn, allow denser surface waters to intrude from above. The result is a so-called convection cell, a zone of continuously recycling and dispersing fluid beneath the shoreline. There is evidence that hypoaigic environments can be associated with significant stratification in hydrochemistry. For example, Ullman et al. (2003) described the effects of mixing between surface waters and groundwaters on nutrient cycling and diagenesis in a (saline) near shore environ-

332 ment, concluding that nitrification is dominant higher in the zone and ammonification deeper in the zone. Duncan and Shaw (2003) reported large variations in rare earth element concentrations across the hypoaigic zone for a coastal salt marsh environment. The mobility and fate of metals advecting with groundwater in hypoaigic sediments has also been quantified (Simpson et al., 2004), with the bioavailability of the metals thought to be much lower in less permeable silts than in coarser sediments. At the same time, spatial and temporal influences can have profound impacts on fluid migration. Fluid exchanges between aquifer and surface waters have been shown to depend strongly on topographic and hydrological settings. C ¸ elik (2002) highlighted recirculation between river water and groundwater in a saline inland setting. Linderfelt and Turner (2001) examined seasonal and tidal nutrient discharge dynamics near river meanders, showing how ammonium may recirculate between fresh river waters and pore fluids through annual cycles. This work was extended to steady saline convection by Smith and Turner (2001). NeilsonWelch and Smith (2001) emphasize the complexity of groundwater flow vectors in the hypoaigic zone, especially in the presence of seasonal salinity variation. Hughes et al. (1998) discussed the importance of diurnal tidal processes in governing small-scale pore fluid exchanges in salt marsh environments and showed that evapotranspiration may be important in controlling saline concentrations in intertidal zones. Other variables may also be important in the hypoaigic zone. Kim et al. (2005) presented data of strong diurnal correlations between groundwater quality measures and sea tides. Novel thermo-convective mechanisms for mass exchange across the sediment–water interface have been examined (Golosov and Ignatieva, 1999) and thermal signatures may be used to identify local seepage flux (Taniguchi et al., 2003; Becker et al., 2004). Scale-dependent physical and biological processes are also important in determining mass transport across the sediment–water interface, e.g. Huettel et al. (2003), while saline fluid displacement by fresh groundwaters can induce significant changes to geochemical evolution (Dai and Samper, 2006). Assessing potential impacts to environmental and human health from these pollutants is also of interest. Field data on the attenuation of industrial pollutants in hypoaigic zones is sparse, making quantification of associated risk and exposure measures difficult, especially in the context of site-specific investigation. As yet, there are no convenient and broadly applicable measures of the natural degradative capacity and robustness of hypoaigic (nor hyporheic) sediments in the presence of contamination. Before these broad measures can be attained it is necessary to develop a conceptual framework for hypoaigic hydrology, the important boundary conditions and regimes that are active in convective surface water– groundwater interactions. In a recent field study of subsurface dissolved hydrocarbon migration near a seasonally saline river, Westbrook et al. (2005) observed that the saline wedge was a significant control on groundwater seepage patterns and had the potential to distort flow patterns deduced from freshwater theories. The present work seeks to build on the results of Westbrook et al. (2005) by clarifying the hydrodynamic mechanisms influencing the hypo-

M.G. Trefry et al. aigic flow regime at the study site. We hope that the outcome of the work is a clearer understanding of the length and time scales involved in the near-shore surface water–groundwater interaction and of the central role played by the saline convection in controlling discharge flow and chemistry. Here, where possible, computer simulations of groundwater flow processes and regimes are compared to measured data from the site, as provided in Westbrook et al. (2005). Of most interest is the regional hydrologic setting of the aquifer and plume, and the basic fluid dynamic processes generated by its interaction with the bounding Canning River. This paper does not consider the development, migration and transformation of the plume arising from the contaminated source zone. Instead the focus is on the hydrological and the fluid environment encountered by the plume on its journey from the source zone to the river. In particular, we note that degradative potential depends on both local hydrochemical and geochemical conditions and the resulting microbial distributions. In this paper, we elucidate some aspects of the richly complex fluid environment in the hypoaigic zone and how this environment fluctuates according to seasonal influences. Natural degradative capacity must necessarily depend on natural cycles of recharge, river condition and regional groundwater flow. In the following sections, we briefly discuss the site characteristics and establish some simple calibrated models of the site. Density-driven simulations are then carried out and the results compared with field measurements. As is appropriate for dissolved-phase BTEX plumes, buoyancy effects associated with the low-density hydrocarbon plume will be neglected in order to concentrate on the much denser saline species. We note, however, that dense contaminant plumes do exist and can exhibit interesting discharge phenomena in their own right (Zhang et al., 2002).

Site description and hydrogeology In this section, we provide a brief description of the site, its surrounds and important hydrological processes, before describing a sequence of simplified hydrological models used to define boundary conditions for a local site area. The reader is referred to Westbrook et al. (2005) for a more comprehensive discussion of site parameters.

Seasonal effects The Perth metropolitan area (see Fig. 1) enjoys a Mediterranean climate with warm dry summers and cool wet winters. Although there is some evidence of a recent reduction in rainfall, possibly due to global warming effects (Pitman et al., 2004), the long-term average annual rainfall is 870 mm yr1 with monthly averages as shown in Fig. 2. Annual evaporation potential is approximately twice the annual precipitation and outweighs monthly precipitation in all months except the May–August winter period (Davidson, 1995, p. 5). Perth itself is situated on the Swan River, which meets the Indian Ocean at Fremantle Harbour, as shown in Fig. 1. The Swan River has several tributaries, including the Canning River, which joins the Swan 15 km upstream from Fremantle Harbour near Waylen Bay at the south-east-

Hypoaigic influences on groundwater flux to a seasonally saline river

333

two measures (see Hamilton et al. (2001) for a useful summary of the Swan–Canning estuarine hydrology).

Basic hydrogeology

Figure 1 The Swan–Canning river system, Perth, Western Australia. The location of the contaminated site is indicated by a black dot (LNAPL source).

Figure 2 Perth mean monthly rainfall (bars, values in mm), estimated salinity of Canning River (open squares) and estimated recharge (filled squares) to the aquifer local to the study site. Over a full annual cycle the recharge function integrates to 25% of total annual rainfall. Rainfall data from the Australian Bureau of Meteorology.

ern edge of Melville Water. During the summer months, when precipitation runoff from the Swan catchments is low, the Swan–Canning river system is saline throughout the area downstream of the Kent Street Weir (Thompson et al., 2003), as shown in Fig. 1. In winter, the flushing effect of the catchment runoff reduces salinity in the Swan–Canning system to brackish levels. A representative annual cycle of local river salinity is shown in Fig. 2, calculated from measurements of river water conductivity according to the 1978 Practical Salinity Scale (PSS, 1978; Fofonoff and Millard, 1983). Under the PSS, a salinity of 35.00 describes seawater at 15 C; we will assume this salinity is equivalent to a dissolved inorganic salts concentration of 35 g L1 to approximate the concentration of seawater (temperature variations are ignored hereafter). Thus in this work we will use PSS to measure salinity (a dimensionless quantity), and dimensions of g L1 to measure concentration, noting the assumed numerical identity between the

Canning Bridge is built across a narrowing of the Canning River near the entrance to Melville Water. Land either side of the river in the vicinity of the bridge is zoned as traffic reserve, light industrial/commercial or public open space. The western bank forms part of a broader peninsula of land, approximately 3 km in length and 2 km in width and bounded by the Canning River to the east, Melville Water to the north and Alfred Cove to the west. The river acts as a discharge sink for the aquifer at the site. Surface water bodies are absent in the peninsula, and much of the area is zoned residential. Peak topographic elevations in the peninsula range from approximately 20 m Australian height datum (AHD) in the north to above 30 m AHD in the south. Throughout the peninsula the superficial aquifer consists largely of Tamala Limestone and associated sands of the Bassendean Formation, and is bounded below at an elevation of approximately 20 m AHD by subcrops of the Kings Park Formation and the Mullaloo Sandstone Member (Davidson, 1995, p. 45), which may be regarded as providing barriers to vertical flow. The particular geologic units in the superficial aquifer are transmissive, with conductivities greater than 1 m/d; in some locations nearer the coast solution features and preferential flow channels have been observed in the underlying Tamala Limestone unit yielding significantly greater conductivities. Groundwater at the site is impacted by dissolved-phase petroleum hydrocarbons, including the BTEX (benzene, toluene, ethylbenzene, and xylene isomers) compounds and naphthalene. These impacts originate from a past unintentional release of light nonaqueous phase liquid (LNAPL) from underground storage tanks approximately 80 m from the river bank. Monitoring activities commenced at the site in 1995. A dissolved-phase BTEX groundwater plume emanates from the source (Fig. 1) and follows an easterly flow path, penetrating under public open space at the western shore of the Canning River just south of the bridge (see Section ‘‘Flow models in plan projection’’). In this vicinity, the superficial aquifer consists of fine- to medium-grained sands, underlain by a clay unit at an approximate elevation of 10.5 m AHD. Porosity of the soil was measured to be 0.2–0.3. The hydraulic gradient from the source to the river foreshore varies seasonally between 5 and 10 mm m1 (0.5–1%).

Tidal observations A limited record (three weeks) of river stage heights in Canning River (Westbrook et al., 2005) exhibited tidal signatures at diurnal and semidiurnal frequencies broadly similar to those observed at Fremantle harbour (Hamilton et al., 2001; Trefry and Bekele, 2004) and elsewhere in the Swan River (e.g. Smith, 1999). Signal correlation analysis showed that fluctuations in river levels at Canning Bridge and at the Perth Gauging Station were well correlated, with the Canning Bridge signal apparently lagging the Perth signal by 30 min, and with the mean elevation

334 at Canning Bridge being several centimetres higher than the Perth Gauge reading over the measurement period. The standard deviation of the river level fluctuations at Canning Bridge is 97% of the value recorded at the Perth Gauge. Hence, with suitable transformation, the Perth Gauge data could be used as a surrogate river stage height signal for the Canning site. Water level oscillations were observed in all site monitoring wells. In topographically simple and constant density systems, the amplitudes and phases of these oscillations may be compared with the corresponding river oscillations to deduce estimates of aquifer hydrological properties. Simple one-dimensional linear propagation is most relevant to situations involving semi-infinite and confined aquifers bordered by fresh water bodies with long straight coasts. These propagation models use the homogeneous aquifer diffusivity, D  T/S = KB/S, as a key parameter. Here, T is the aquifer transmissivity, K is the hydraulic conductivity, B is the (constant) saturated thickness of the aquifer, and S is the storativity. Analytical solutions for these simple propagation models were reported by Jacob (1950) and have been used widely ever since. In the present case, the Canning River provides neither a freshwater (nor constant density) boundary to the aquifer, nor a perfectly linear reach and therefore the simple linear tidal-propagation models are not entirely valid. However, as shown recently by Trefry and Bekele (2004), corrections due to density effects are modest, as are nonlinear effects associated with unconfined aquifers. Shallow bank slopes

M.G. Trefry et al. (i.e. near horizontal) can provide significant perturbations to linear tidal propagation theories, especially within a few metres of the shoreline (Ataie-Ashtiani et al., 1999; Teo et al., 2003). If we confine attention to distances greater than a few metres from the shoreline and less than a 100 m or more (negating perturbations from curving river reach), we may expect that a simple tidal propagation analysis at the Canning River site may provide us with at least an order of magnitude estimate of local aquifer diffusivity. Fig. 3 shows the amplitude spectrum measured over 77 days at bore MW33, located 8 m inland from the river bank (see Fig. 5). The spectrum, corrected for finite-sample tapering, is dominantly forced at the diurnal O1 and K1 frequencies, with minor semidiurnal constituents. This spectrum can be compared with the surrogate river elevation spectrum to generate amplitude attenuation factors, yielding average diffusivity estimates of 370 m2 d1 and 350 m2 d1 for the O1 and K1 tidal constituents, respectively. For wells MW42 (48 m) and MW18 (50 m), which have time series of shorter duration (<10 days), the O1 and K1 peaks were not resolved adequately to permit a spectral analysis. Using the measured aquifer saturated thickness of approximately B = 10.5 m, the MW33 spectral estimates yield K/S values of 35 m d1 and 33 m d1 for the O1 and K1 modes, respectively.

Flow models in plan projection A simple idealised groundwater flow model for the peninsula was developed in plan view in order to calibrate large-scale transmissivities and to assess likely average flow directions near the plume. The groundwater modelling package Feflow (Diersch, 2004) was used as a solver for both the present flow simulations and the following flow and density simulations. Feflow has a well tested range of solution methods for density coupled systems (Oswald and Kinzelbach, 2004; Diersch and Kolditz, 2002). The peninsula shoreline was digitised from maps and entered into the model as a fixed head (first type) boundary condition with elevation 0 m AHD, which is slightly above the Fremantle mean sea level elevation of 0.04 m AHD, and below a measured mean elevation at Canning Bridge of 0.11 m AHD. The southern boundary was chosen to lie along an 5 m AHD contour of head (Water and Rivers Commission, 1997, p. 62). The grid contained 2982 nodes and 5676 triangular elements. Spatially uniform recharge was estimated to be distributed evenly throughout the year, amounting to 25% of the annual rainfall (equivalent to 0.6 mm d1), with a uniform aquifer basement elevation of 20 m AHD. Elevation of the phreatic surface was calibrated as a function of hydraulic conductivity.

Effective peninsula hydraulic conductivity

Figure 3 Amplitude spectrum of detrended water levels measured at the Perth Gauging Station and at well MW33 during the period 7 July to 22 September 2002.

Simulations were focused on reproducing the 1 m AHD head contour in the immediate vicinity of Canning Bridge (Water and Rivers Commission, 1997, p. 62). Fig. 4 shows the elevation of the phreatic surface in a northern subregion of the peninsula grid calculated using a uniform hydraulic conductivity of K = 5 m d1 throughout the model domain. The model shows tolerable agreement with the head contours indicating that the chosen conductivity is a reasonable

Hypoaigic influences on groundwater flux to a seasonally saline river

335

Figure 4 Steady heads and streamlines for a two-dimensional phreatic model of the peninsula. Predictions of annual mean head are shown by solid contours (expressed in metres AHD), and are compared with estimated (dashed) maximum contours published in Water and Rivers Commission (1997). Also shown are mean annual streamlines relevant to the Canning Bridge area (arrowed lines). The heavy black line shows the path of Canning Highway and the shaded grey region shows the local site area.

approximation for the assumed recharge. The position of the 1 m AHD contour is insensitive to changes in K of ±30%, and to changes in the river boundary stage height of ±0.1 m. The estimated value of K = 5 m d1 is somewhat lower than a value of 40 m d1 used by Smith (1999) in a study of flow processes in an alluvial (river bank) aquifer 10 km to the north-east on the Swan River, but is consistent with values for aquifer material consisting of a mixture of fine-medium sand with carbonate rock (Kruseman and de Ridder, 1994, p. 21). Matches between the simulated and estimated 1 m AHD contours are closest in the vicinity of Canning Bridge. Fig. 4 also shows the mean streamlines calculated from the simulation. The streamlines show spatial curvature that is controlled by the topographic salient in the river bank that forms the narrows at the bridge location.

Site groundwater flow model The previous results for the peninsula model were used to define boundary conditions for a localised site grid (in plan projection), allowing processes relevant to the hydrocarbon plume to be studied more efficiently and at higher resolution. In the absence of local pumping, the 1 m AHD contour provides a useful datum for a fixed head or fixed flux boundary condition along the western edge of a site grid. The northern and the southern boundaries can be chosen to lie along mean streamlines, which define zero flux conditions. The river bank conveniently defines the eastern boundary with a fixed head condition. The two-dimensional model for the site will be used to calibrate local hydrological

Figure 5 Part (a) shows the location of the local site model, together with boundaries, mean annual streamlines, and the estimated location of the hydrocarbon plume (pale zone). Part (b) shows the calibrated steady head solution for the site, and the locations of the calibration wells. Heads are expressed in units of m AHD.

parameters in anticipation of more detailed three-dimensional models of local flow and transport processes. Fig. 5a shows the site grid chosen using the 1 m AHD head contour and relevant streamlines as groundwater boundaries, together with the local topographic setting. The river provides a fixed head boundary at elevation 0.11 m AHD, corresponding to the local mean elevation. The grid contains 1834 nodes and 3450 triangular elements. Also shown in Fig. 5a, by pale fringes, is a simulation of the plume extent generated using these boundary conditions. This simulated plume is congruent with the measured plume at the scale of the figure. The use of a fixed head upgradient boundary for groundwater flow allows large influxes or effluxes of fluid to the model in the case that the aquifer water balance is stressed, e.g. through local pumping, or

336 river flooding. A common alternative would be to specify a fixed flux through the upgradient boundary, i.e. a secondtype condition; however this requires a greater knowledge of local hydraulic conductivity. In both cases, the assumed vertical recharge rate will affect the steady water table levels throughout the model domain. In this work, the fixed head condition was chosen, with the proviso that any predictions of the modelling would depend on the absence of local pumping and flooding.

Calibration of the hydraulic conductivity As indicated earlier, the site grid was used to calibrate an effective uniform aquifer hydraulic conductivity by comparison with local site measurements of standing water levels in boreholes. Fig. 5b shows the locations of 6 bores available to be used in the calibration process. Standing water levels were measured in these bores in October 2000 and January 2001, showing consistent reductions of approximately 0.14 m over the 3 month period. For calibration purposes, the measurements were averaged over the dates to yield a mean elevation for each bore (MW7 0.54 m AHD, MW14 0.53 m AHD, MW32 0.28 m AHD, MW33 0.13 m AHD, MW35 0.24 m AHD, MW42 0.49 m AHD). The hydraulic conductivity parameter K for the steady site model was calibrated (using the PEST module for parameter estimation in Feflow) against these mean values. The calibrated heads were as follows: MW7 0.50 m AHD, MW14 0.56 m AHD, MW32 0.31 m AHD, MW33 0.20 m AHD, MW35 0.25 m AHD, MW42 0.43 m AHD, yielding a root-mean-square residual error of 0.04 m. The best estimate for the conductivity was K = 2.6 m d1, which is slightly lower than the value estimated from the peninsula model, but still in line with the range of values usually indicated for medium-fine sand. The 95% confidence interval is large (0.3–30 m d1), implying that the global-property model used here cannot describe perfectly the detailed water table slopes at the site. Returning to the estimated tidal diffusivities (350– 370 m2 d1; Section ‘‘Tidal observations’’), inserting the above K values leads to storativity values below the measured porosity of the aquifer matrix of 0.2–0.3. For example, with K = 2 m d1 a storativity of S = 0.057 is indicated, while K = 10 m d1 yields S = 0.28; larger values of K are contraindicated with S values above the local porosity. Hence, within a reasonable error margin, the near-shore tidal propagations are not inconsistent with these head calibrations of K and indeed the conductivities expected for medium-fine grained sands. We must be mindful that this groundwater spectral data comes from a well that is close to the river bank, and hence may be affected by departures from the assumptions of the ideal one-dimensional tidal propagation theory. Even so, the various estimates are in reasonable agreement and the indicated interval of K is within reasonable bounds. Uncertainties in the global K value will alter the scale of simulated flow velocities and hence solute distributions in transport simulations. For example, larger K values will increase flow rates and dispersion of the saline wedge. However, the main objective of this study is to understand the nearshore hypoaigic processes and their potential influences on the fate of dissolved hydrocarbon species. Hence, if

M.G. Trefry et al. we confine attention to the physical processes involved and remember that the global K value is simply an approximation of necessity, we can learn much from the consequent simulations of the flow and transport in the hypoaigic zone. In particular, we are interested in seasonal processes and variations. In this sense, accurate absolute values for flow velocities are of less concern than understanding how flow velocities may change throughout the annual cycle.

River-aquifer models in vertical section Although the detailed cross sectional bathymetry of the Canning River is not known, it is well understood that, where bank slopes are low and in the absence of confining beds, most discharging groundwater flux from shallow aquifers occurs near the river/lake water line, and that discharge fluxes beneath the centre of surface water bodies are very low (Winter, 1999; Townley and Trefry, 2000). Hence, for the present purposes, it will be sufficient to determine the river bank slope; the full details of the river bed elevation far from shore are not as important. The bank slope values used in the vertical section models were 6.6 at the intertidal face, extending horizontally for 4 m toward the river below the mean water line, followed by a much longer and gentler slope of 1.6 (Fig. 6). These slopes are representative of measured local river bank topography; the wedge dynamics are not thought to be sensitive to the deep river bed slope. Seepage dynamics will depend on the slope of the intertidal face (see Mao et al., 2006).

Freshwater flow As a first step in understanding the hydrodynamic processes at the study site, we consider steady two-dimensional (vertical section) simulations. The sections are aligned with a notional streamline from the contaminant source to the river. Key issues are the steady dynamics of the groundwater system subject to a saline boundary condition. The model grid contains 12,826 nodes and 24,896 triangular elements, and is refined for greater spatial resolution near the river bank and shore line. The uniform hydraulic conductivity K was set to 5 m d1. The constant river elevation was assumed to be 0.11 m AHD, and the aquifer basement was set at 11.0 m AHD. A simple van Genuchten saturation model was used to describe the capillary fringe and vadose zone water saturations, with parameter values of n = 1/ (1  m) = 1.964, scaling coefficient As = 4.1 m1, residual water saturation Sr = 0.0025 with zero air entry pressure. Contributions from recharge in the 2D section near-field were ignored. For this simulation the effect of the river was represented as a set of fixed head boundary conditions along the river bed. In order to assist with interpretation of fluxes in these two-dimensional models, the inset in Fig. 6a defines three sections (A, B, C) for which fluxes were extracted from the simulation results. Focusing and seepage Fig. 6b shows a steady phreatic simulation from a twodimensional (vertical section) model constructed for the site, with measured ground and river bed elevations near

Hypoaigic influences on groundwater flux to a seasonally saline river

337

a

b

c

Figure 6 Two-dimensional vertical section models on a transect through the study area, showing the effect on the groundwater flow regime of including density effects from the saline river. Arrowed lines indicate flow paths calculated from the numerical solutions using a particle tracking procedure.

the shore. The vertical focusing effect of the near-shore river bed on groundwater flow paths (calculated by particle tracking) is clear in Fig. 6b. This effect is most pronounced for shallow bed slopes, and is absent where surface water bodies penetrate vertically to the base of the aquifer, as is sometimes the case for mining voids. Importantly, groundwater fluxes are greatest where the flow paths converge at the river bank, and least underneath the river further offshore. Based on this simple freshwater model, we may expect that any dissolved conservative species originating upgradient and advecting along groundwater streamlines would preferentially discharge to the river along a quite narrow zone on the river bed near the river water line. The water table gradient in this freshwater model is 2.8 mm m1 (0.28%), which is a factor of two lower than the measured gradients at the Canning River site. The horizontal groundwater flux in Section A of the model is q = 0.168 m2 d1. Of this total groundwater flux, approximately 60% discharges at high fluid velocity to the river in a narrow seepage zone (Section B, 5 m width) at the shoreline. The maximum velocity at the shoreline is 0.2 m d1, where the fluid velocity is calculated by v = K/h oh/ox, with the porosity h = 0.3. The remainder of the groundwater travels extremely slowly along flow paths at the bottom of the aquifer before moving vertically into the river bed at greater distances from the shore.

Saline flow Inclusion of density effects from the saline river waters leads to a more realistic vertical section model, especially in defining the groundwater seepage face (Zhang et al., 2001). Density driven flow and transport was simulated using the default Boussinesq model for density coupling in Feflow,

with a linearised concentration–density relationship for a single conservative saline species, and viscous flow corrections neglected. Due to the relatively small number of elements, a direct matrix solver was used in preference to iterative techniques. The river boundary nodes were set to fixed concentration conditions equivalent to the estimated maximum river salinity (35), yielding a river water density of 1.026 g cm3. Temperature was assumed to be constant at 25 C. Since Feflow assumes that boundary conditions are expressed in terms of equivalent fresh water heads, each of the fixed head nodes along the river bed boundary were adjusted to account for the weight of the dense water column above, assuming a steady river elevation of 0.11 m AHD. Longitudinal and transverse dispersivities for the saline transport were initially set at aL = 5.0 m and aT = 0.5 m, respectively. Explicit steady solution techniques were unable to provide accurate solutions. Instead, the simulations were run in transient flow and transport mode with an automatic predictor–corrector time stepping algorithm for at least 15 years of simulated time, until the saline distribution had stabilised, see Fig. 6c. The shading in Fig. 6c indicates the calculated fluid salinity. Convection The simulation yields a wedge penetration of 19 m (measured at the salinity 20 location on the aquifer base), in reasonable agreement with field measurements indicating a wedge extending more than 15 m into the aquifer. A conventional vertical section sharp-interface model (Glover, 1959) yields the depth z of the interface of the dense wedge at a distance x into the aquifer from a dense boundary, as given below:  2 2Axrh Arh z2 ¼ ð1Þ þ   a a

338

M.G. Trefry et al.

Here $h is the effective slope of the water table generating fresh water flow to the boundary, q=qmin ¼ 1 þ a represents the normalized density of the saline fluid, A is the unit thickness of the section (in the y-direction), and the aquifer depth is assumed to be infinite. Applying the simulated (Section ‘‘Focusing and seepage’’) effective slope of $h = q/AK = 0.168/5 = 0.0336 and density of 1.0264 yields an intercept at 11.1 m AHD of x = 48 m. The disparity with the simulation prediction of 20 m is usually reconciled in terms of the diluting effect of dispersive mixing in miscible flow models (Volker and Rushton, 1982). The presence of this wedge deflects the fresh groundwater flow upwards toward the river shore. The groundwater dilutes the front of the saline wedge via mechanical dispersion, creating a zone of variable density, with low density and upward flow on the landward side of the fresh-saline interface, and medium density on the river side. Further offshore, high density saline water flows downwards. Thus the density effects generate a recirculating fluid convection cell in the aquifer under the shoreline. In this steady simulation, fluid velocities are greatest at the margins of the convection cell, and least (near stagnant) at the centre of the cell. Mass transfer through the centre of the convection cell is largely diffusional. This general hydrodynamic picture is consistent with the results of Smith and Turner (2001), which were supported by field measurements at the nearby Swan River. Density-driven convection cells arise from the gravitational instability of saline (dense) waters overlying fresh (less dense) waters, and are commonly classified in terms of the equivalent cellular Rayleigh number, Ra, given by (Wooding, 1989; Diersch, 2002): Ra ¼

KaB hDmol

ð2Þ

where K is the hydraulic conductivity of the aquifer (assumed isotropic), B is the aquifer thickness, h is the effective porosity and Dmol is the effective coupled free diffusion coefficient of the density-determining ions in water. This dimensionless group is appropriate to studies of free (unforced) convective overturn in idealised rectangular porous domains with the dense fluid source at the top boundary and zero lateral pressure gradient. For the present study involving mixing of saline and fresh fluids, the density coefficient a is given by (qmax  qmin)/qmin, where qmin is the density of the reference fluid, i.e. fresh water. Typically, in this model, free convection cells occur for Ra P 4p2, with simple diffusional processes holding at lower values of Ra. For higher values, above Ra  390, the convective processes become unstable. This result neglects the influence of lateral gradients (e.g. regional groundwater flow) on the convection cell. Simmons and Narayan (1997) provide a useful correction for lateral flow, yielding a modified Rayleigh number, Ra 0 , given by Ra0 ¼

KaB ðDmol þ aL V amb Þ hDmol ðDmol þ aT V amb Þ

ð3Þ

where aL (aT) is the longitudinal (transverse) dispersivity and Vamb is the ambient regional groundwater flow speed in the horizontal plane. Applying this formulation to saline disposal basin scenarios, Simmons and Narayan (1997) found

that stable convection was obtained for Ra 0 6 1250; for higher values unstable fingering occurs. It is instructive to consider how the Canning River site rates in terms of these convection numbers. In Section ‘‘Flow models in plan projection’’, the hydraulic conductivity of the superficial aquifer was estimated to be close to 5 m d1 and B = 10.5 m. The following parameter values are assumed: porosity h = 0.3, molecular diffusion coefficient Dmol = 1.3 · 104 m2 d1 (Smith and Turner, 2001),  a ¼ 0:026 and a regional hydraulic gradient of 7 mm m1. With these values, the cellular Rayleigh number Ra  35,000 and the modified Rayleigh number Ra 0  350,000, are well above the respective critical values of 390 and 1250 for the onset of unstable fingering. The order of magnitude difference between the two Rayleigh numbers results from the ratio of dispersivities aL/aT because of the large regional groundwater flux to the river bank (Simmons and Narayan, 1997). This dispersion ratio has not been measured at the site and was assumed to be equal to 10, in line with common practice (Gelhar et al., 1992). Convective instability The high values of Rayleigh number indicate that the nearbank aquifer is likely to be prone to unstable fingering and convective circulation. The simulations presented in this paper do display evidence of fingering processes, especially in the later three-dimensional simulations. However, discussion of the fingering dynamics is kept to a minimum. This de-emphasis is deliberate – so far we have been most interested in understanding the long-term behaviour of the saline distributions. Instabilities in convection predominantly arise at the initial intrusion of dense fluids into the aquifer; the presence of established wedges and convection cells would inhibit fingering zones since the vertical density contrasts would be reduced. Nevertheless, there is potential for small-scale fingering to occur at a semidiurnal time scale due to tidal flooding of the seepage face. This, in turn, may have the potential to perturb the groundwater flows to the seepage face at this rapid time scale. Large seasonal variations in river salinity may also induce fingering from overtopping of less dense pore fluids by saline surface water. However, as Simmons and Narayan (1997) remark, even small variations in aquifer conductivity, such as might occur through layering or lenticular features in the porous medium, may either promote or reduce the instabilities. Numerical instabilities are also potential sources of fingering effects in density coupled models. The present simulations use a relatively coarse grid structure and, in the absence of more detailed information on the fine structure of the porous medium, we confine attention to bulk and longer-term behaviour. Where fingering is apparent in numerical results it will be noted in the text without further emphasis. The reader is referred to Otero et al. (2004) for a discussion of high-Rayleigh number convection. Dispersive effects We may compare the long-term flux estimates from this density driven simulation with those of the earlier freshwater model. The dense wedge causes all the aquifer ground-

Hypoaigic influences on groundwater flux to a seasonally saline river

339

Figure 7 Steady convection cells under the river bank calculated in two-dimensional section with different dispersivity values (a), (b) and (c). Part (a) corresponds to the density-driven solution of Fig. 6c. Solid lines indicate flow paths calculated using a particle tracking procedure. Salinity distributions are depicted by shaded contours. The location of the salinity 20 contour at the base of the aquifer is indicated in each case by an arrow. The dashed line shows the sharp interface solution of Glover (1959). Part (d) shows salinity variations calculated along a common vertical profile.

water flux to exit into the river at a narrow seepage face. However, the flux through Section B is significantly higher than for the fresh water simulation (see Fig. 6).

A series of simulations with smaller dispersivities show that the toe of the saline wedge penetrates further into the aquifer with less dispersion (26 m and 33 m penetrations

340 for (aL, aT) = (2.5, 0.25) and (1.25, 0.125), respectively). Reassuringly, the calculated wedge shape also tends toward the Glover sharp-interface prediction (dashed lines in Fig. 7) as the dispersivities are reduced. However, a key feature of the sloping bed simulations is a relatively fresh (i.e. brackish, salinity 10) cell underneath the river bank. This cell becomes more pronounced as the dispersivities decrease. The implication is that a sampling well drilled through the near-shore river bed would sample at first highly saline water, then brackish water, and then highly saline water again as the sampling depth increased. The disparity between the saline and brackish zones depends on the dispersivities – lower dispersivities yield brackish cells that tend toward fresh water concentrations, as shown in Fig. 7d for a vertical profile extending from the shoreline. Brackish cells below the river bed were also seen in the field investigations (Westbrook et al., 2005). Measurements indicate a zone of near-fresh water several metres offshore and several metres below the river bed, see Fig. 8 (cf. Fig. 11 of Mao et al. (2006), where the re-entrant zone is generated by tidal saline flooding). Salinities as low as 1 were estimated from fluid conductivities in this zone, perhaps indicating that local dispersion is small. Similarly, there was evidence of saline fluid several metres thick at the bottom of the aquifer 10 m inland from the Canning River shore. This is consistent with the simulated wedge thicknesses, especially for dispersivities at the lower end of the simulated range, before the wedges pinch out further into the aquifer. There is no independent data from the site which can be used to estimate the dispersivities for the surficial aquifer, although other evidence from medium sand aquifers seems to support the use of saline dispersivities at much lower values. Lebbe (1999) showed how bore resistivity measurements of hypoaigic saline distributions could best be fitted with millimetre-scale dispersivities and significant vertical anisotropy in K. Zhang et al. (2001) measured sub-millimetre saline dispersivities for a fabricated glass bead aquifer. For the present simulations, the grid discretizations effectively controlled the choice of dispersivities. Numerical stability requires that dispersivities be larger than characteristic dimensions of the finite elements so that grid Peclet numbers remain small.

Figure 8 Measurements (circles) of BTEX concentration (with estimated contour lines) and electrical conductivity (with estimated shaded contours) taken at the field site in September 2001. Note the re-entrant zone of fresh water between the saline river bed and the wedge below. Figure reproduced from Westbrook et al. (2005). For the Canning site at 25 C, one salinity unit is approximately equal to 1615 lS cm1.

M.G. Trefry et al. Seepage flux and convective recirculation Seepage fluxes were predicted for the density-driven simulations. In all these simulations, the regional groundwater flux to the river, as measured across a vertical section upstream from the river bank, was 0.146 m2 d1, as compared to 0.168 m2 d1 for the earlier fresh water simulation. This is attributed to the displacement effect of the saline wedge effectively providing an increased head boundary condition (at the seepage face) for the discharge of regional groundwater flow. In order to account for this change in absolute flux, relative fluxes are listed here for Section B (expressed relative to the flux through Section A, cf. Fig. 7). The relative fluxes were 60% for the freshwater simulation and, for the density driven simulations, 97%, 104% and 108% for (aL, aT) = (5.0, 0.5), (2.5, 0.25) and (1.25, 0.125), respectively. Here, it is seen that whilst 60% of the regional flux discharges through the 5 m long Section B in the fresh water model, almost all of the flux discharges through that section in the saline models. In fact, there is an apparent (but slight) oversupply of fluid to the discharge section as the dispersivities decrease. The difference is contributed by the convection cell as a recirculation of river water. As discussed by Smith (2004), this recirculation effect is driven largely by dispersive mixing between the fresh and saline fluids across the radius of the convection cell. In particular, Smith (2004) defines a percentage saline circulation ratio PSC as Kz afW ð4Þ q where q = Kxoh/ox is the regional groundwater flux moving horizontally toward the surface water body, Kz is the vertical component of the hydraulic conductivity tensor, and the empirical dispersion functional fW is given by  1  eaW ; 0 6 W 6 40 ð5Þ fW ¼ 40 < W < 20 b=W c ;

PSC ¼

and W = B/aT. For aL/aT = 10, the empirical parameters a, b, c have the values 0.11, 3 and 0.3, respectively. In the absence of evidence to the contrary, we can assume isotropy of the hydraulic conductivity, that is Kx = Kz = 5 m d1; this yields a value of PSC  2% for the Canning River site (the derived PSC value declines for Kx > Kz). This is interpreted to mean that 98% of the fluid exiting the aquifer through the river bed originates from upgradient (inland) in the aquifer; the rest is recirculating (and variably dilute) river water returning to the river via the convection cell. This value of saline recirculation is much lower than the 20% value estimated for the Swan River bank near Ron Courtney Island (Smith, 2004), however this is explained by the proportionally lower groundwater flux at the latter site. The key controls are topographic – the Ron Courtney Island site is situated at a part of the river where head gradients are much lower. At the Canning River site, the contamination plume (Fig. 5a) is driven by a significant inland groundwater mound and meets the river essentially at a point of inflection of bank curvature just to the south of the Canning Bridge salient. Groundwater fluxes are much higher there and tend to dominate the seepage budget. In any event, the PSC value of 2% is in approximate agreement with the modelled estimates of saline contributions to seepage supplied above.

Hypoaigic influences on groundwater flux to a seasonally saline river We can gauge the length of seepage face utilised by the regional groundwater flux by post-processing the simulation solutions. Here, a particle track commencing in fresh water on the bottom of the aquifer near the end of the saline wedge was traced to the river bed. The terminal point of this track was then used to calculate a distance from the simulated shore waterline. The seepage face lengths thus calculated are 3.7 m, 2.2 m and 2.0 m for (aL, aT) = (5.0, 0.5), (2.5, 0.25) and (1.25, 0.125), respectively. These results show that this ‘‘fresh seepage’’ distance declines with dispersivity, mainly because the transition between fresh and saline water becomes less blurred. Note that the total length of seepage face available to the regional groundwater is longer than the quoted distances because of the possibility of seepage from above the river watermark. However, the present model is not designed to predict this effect accurately (see, e.g., Thorenz et al. (2002) for a discussion of the complexities of near water table convective flow), hence the use of the waterline as a measurement datum. Overall, the density-driven simulations show that virtually all the regional groundwater flux discharges to the river along a narrow zone near (above and below) the waterline on the shoreline (as discussed by Volker et al. (2002)). The length of this discharge zone is affected by the choice of dispersivities in the simulation. As remarked earlier, field observations of salinity distribution tend to favour dispersivities at the lower end of the range considered, implying that the discharge zone may be less than 3 m in length.

Characteristic time scales in the saline system It was seen earlier how diurnal fluctuations in river stage height contributed to correlated small fluctuations in heads measured at locations in the aquifer. These fluctuations were used with the fresh water models to estimate aquifer hydraulic conductivity. It is natural to consider what effect such tidal fluctuations and longer-term influences have on the gross flux dynamics within the saline aquifer models. Long-term influences Firstly, we consider processes at a length scale of tens of metres. These include regional fluxes, saline wedges and convection cells. Simulations show that the characteristic relaxation times for these attributes are of the order of years, meaning that the simulations would be unable to resolve variations in these attributes on diurnal timescales. In other words, the saline wedges and convection cells would appear essentially unchanged as tides vary from day to day. The wedges and convection cells are much more sensitive to seasonal and inter-annual variations. For example, consider the steady saline wedge penetration distance of 33 m (Section ‘‘Dispersive effects’’ and Fig. 7) calculated for a steady river salinity of 35. Whilst providing a useful conceptual scale for the system, this wedge calculation is based on unchanging river water salinity. In reality, the river salinity varies throughout the year, providing a river bed boundary condition of varying salinity. In turn the wedge will constantly be perturbed, meaning that the 33 m penetration distance will be an upper bound

341

for a wedge generated by stationary (constant mean) salinity variations. We can estimate a lower bound for the wedge penetration distance by taking the high-salinity numerical solution and exposing it to a river boundary condition corresponding to the minimum annual river salinity of 5. This will erode the wedge slowly, with the loss of mass from deep in the aquifer to the river bed being increasingly diffusion-limited. Running such a model for a simulated period of 10 years yields neither a true steady state (i.e. with maximum concentration in the model domain of 5 g L1), nor a short term (e.g. 1 month) perturbation of the 35 g L1 wedge. The 10-year perturbation shifts the 20 g L1 location of the toe of the wedge almost 20 m towards the river bank. Therefore we should expect the wedge penetration distance in a seasonal model with varying river salinity to fluctuate in a range bounded (loosely) by 13 and 33 m from the mean river shoreline. A more precise estimate of the seasonal effects upon the saline wedge can be gained by performing an explicit simulation of a seasonal river salinity cycle. A corresponding simulation was performed using Feflow. This simulation allocated a unique seasonal salinity variation to each of the 124 nodes along the river bed; zero recharge variation was assumed. This salinity variation was taken from the salinity cycle of Fig. 2 and incorporated the depth below mean river elevation at each node. The simulation commenced with the steady wedge for salinity 35 and was integrated for 30 annual cycles until the total saline mass in the aquifer had reached a quasi-stationary state, changing by less than 0.07% per annual cycle. The aquifer saline mass started at 13,445 kg at year zero and converged to 11,544 kg at the end of the simulation, i.e. the quasi-stationary state of the wedge contains 14% less mass than the maximum concentration state. This simulation predicted the inter-annual variation in wedge toe location (salinity 20 contour) to have amplitude no more than 0.14 m around a mean penetration distance of 27.8 m from the river shoreline. At the same time, there were significant variations in saline interface shape and position higher in the aquifer near the river shoreline. In general the simulated saline distribution higher in the aquifer responded more quickly to variations in river salinity than did the wedge toe. The wedge toe seemed to expand in the winter (6 months after summer’s saline river conditions), and contract during summer. These modest variations in wedge shape were generated purely by river salinity cycles. A more representative model would also include seasonal variations in regional groundwater discharge. Such a simulation was also run. The variations in discharge were generated by altering the upstream fixed head boundary condition by ±30% (0.7–1.3 m AHD) over an annual cycle with extrema in March (minimum) and October (maximum). This simulation was again run to a quasi-steady state after 10 annual cycles, with the total wedge mass converging towards 11,375 kg, i.e. 1.5% less than the river-only variation simulation. In this simulation, the salinity 20 location at the wedge toe experienced fluctuations of amplitude 0.6 m. The wedge exhibited a mean penetration distance of 27.4 m, 0.4 m less than for the river-only variation simulation. Fig. 9 shows the variation of salinity contours as a function of season for the seasonal model. The con-

342

M.G. Trefry et al. mean transition point, and may also cause the groundwater seepage face to migrate up and down the shore. Mixing will be of two main types. Fresh-saline mixing will occur as the seepage face moves. A more subtle mixing involving ‘‘young’’ saline river water and ‘‘old’’ saline groundwater will also occur within the convection zone. Assuming that the steady flow regime is perturbed by a change in elevation of Dh = 0.1 m for Dt = 6 h, and that the K value of 5 m d1 applies, an exchange process may occur in the first 1 m of sediment over a zone of the order of Dx  K$hDt  10 cm in length. This is similar to the length scale discussed by Webster (2003) in a study of wave-enhanced diffusion in marine sediments, and by Huettel et al. (2003) in connection with sediment biogeochemical reactions. These concepts are supported by the numerical study of AtaieAshtiani et al. (1999), where it was found that tidal fluctuations had greatest effect near the seepage face and minimal effect at depth in the aquifer, especially where the aquifer depth is much greater than the tidal amplitude. These arguments set a scale at which further studies of high frequency transients may be carried out, both in the field or by simulation. It is beyond the scope of the present simulation study to address such issues in detail.

Figure 9 Groundwater salinity contours calculated from seasonal variations in river water salinity and regional groundwater discharge. The deep saline wedge is effectively immobile, although seasonal variations in the saline distribution near the river bed are strong.

tours just below the river bed vary strongly throughout the year, as the river water salinity moves between summer highs and winter lows. The domain most affected by these seasonal variations is the top 2 m of the saturated zone beneath the river bed. Similarly, since the regional fluxes depend on the location and concentration of the wedge, we may expect that the gross regional fluxes and seepage locations (i.e. at the shoreline) will not be varying significantly from day to day with normal tidal variations. Short-term influences Now, consider processes at sub-metre scales near the river bed. Far from the shoreline, the river bed experiences pressure from a saline column of water above it. Tidal variations in this column height vary from 10% to 1% or less, according to the local depth. This change in column height will induce small changes in infiltration rate through the bed into the slowly moving saline convection cell. The sense of infiltration will not change. As we move closer to the shoreline there is a change in sign of the infiltration rate, from a recharge (to the aquifer) sense to a discharge sense (to the river). The change of sense is located at the centre of the intersection of the convection cell with the bed. The location of this transition point is sensitive to the regional groundwater flux, river elevation and density, dispersivity and bed slope. As the river elevation fluctuates, the transition point will shift on the river bed, toward the shore as the river deepens, and back again as the river shallows. This process will induce mixing in the bed sediments near the

Fluid ages As an independent check on these arguments, we calculate the fluid age or residence time using an elegant algorithm introduced by Goode (1996). The fluid age calculation may provide extra information on the origins of fluids mixing in the hypoaigic zone which, in turn, may potentially be important in assessing the natural hydrocarbon degradation potential of the riverine seepage face. For example, older fluids may have significantly higher reductive capacities than younger and more recently oxygenated recharge fluids. Saline riverine fluids will also mix in the hypoaigic zone, potentially making this zone one of complex chemical and biological activity even without the contribution of a dissolved contaminant load. Feflow implements the fluid age algorithm for density coupled simulations by introducing a third dependent variable (in addition to the standard head and concentration equations) for the residence time in the guise of a temperature equation. The solid phase heat conduction and heat capacity is set to zero, the fluid heat capacity is set to unity and the fluid heat conduction is set to zero. The remaining non-zero parameters are the porosity and thermal dispersivities and a source term set globally across all elements equal to the porosity. This source term is the fluid aging mechanism in the thermal transport equation. Boundary conditions are first type with value zero (age in days) along all boundaries where fluid enters the domain and zero-flux conditions on the other boundary nodes. The Feflow implementation cannot support fluid age calculations for unsaturated problems, so we make a fully saturated assumption for our two-dimensional slice, and apply a steady recharge rate to the phreatic surface at 25% of the mean annual total, calculated daily. The river salinity is set to 35. Fig. 10 shows two renderings each of the solutions for two different dispersivity settings; renderings (a,c) with linear contour shading and renderings (b,d) with logarithmic contour shading. The shaded distributions indicate the mean age of fluid at each point, incorporating processes of advec-

Hypoaigic influences on groundwater flux to a seasonally saline river tion, mechanical dispersion, molecular diffusion and physical dilution. Renderings (a) and (b) pertain to the default dispersivity case, aL = 5.0 and aT = 0.5 m. Rendering (a) shows that fluid age increases with depth below the land surface up to approximately 10 years toward the base of the aquifer. Young recharge water is clearly visible with light grey tones, grading to darker tones at the base of the aquifer. The vertical gradient of fluid age is determined by the dispersive mixing of fluids as they migrate toward the river bank. Similarly, young river water is shaded in light grey, but there is a zone of much older (darker shades) water underneath the river at the base of the aquifer. Fluid ages here reach 10,000 days (almost 30 years), significantly older than the oldest water in the terrestrial zone. As we have already seen, this aged fluid zone corresponds to the dense wedge which is stationary at the aquifer base. The reason why the ages within this zone are not larger is due to the combined effects of mechanical dispersion and molecular diffusion. Even though the saline concentration profile is stationary there is still sufficient dispersion and diffusion to refresh the fluids within the dense zone every 30 years or so. Rendering (b) uses a logarithmic age scale for the same solution as (a), allowing short-term processes to be seen more easily. In this case, it is clear that convective recirculation in the top few metres underneath the seepage face occurs at a much shorter time scale than the dispersive and diffusive processes in the saline wedge. From Fig. 10b, we can estimate the fluid age in the convective recirculation near the sediment–water interface as being of the order of a month or more. As we look closer to the seepage face this recirculation time lessens, however the previous arguments are confirmed – full convective recirculation times vary between months and years depending on which travel path is chosen. If we move toward a finer dispersivity model with aL = 1.25 and aT = 0.125 m, as presented in renderings (c) and (d), we see a more realistic story. Rendering (c) shows a sharper and more persistent delineation between aged fluids deeper in the aquifer and younger fluids near the water table. This is especially pronounced further upgradient from the river bank: there is likely to be significant age disparities between deep aquifer fluids and the recharged fluids. For example, the deep aquifer fluids reach ages of almost 90 years and as the groundwater flux exits the seepage face the mean fluid age is still of the order of 20 years. In this rendering, the age of the saline wedge is approximately equivalent to the age shown in rendering (a) at 30 years. Rendering (d) shows fluid in the convection cell to have age of several months to a year. Overall, the fluid age calculations show that the average age of the fluid at the seepage face is several decades or more, consisting of much older groundwater fluids diluted with younger recharge waters. The convection cell incorporates relatively young saline fluid (less than one year of age), whilst the saline wedge itself is replenished every few decades by dispersive and diffusive mechanisms. Finally, we note that diurnal tidal fluctuations, together with fluctuations in solar radiation, wind stress and vegetative uptake also have the potential to influence the water distribution in the unsaturated zone. The water table defines the upper extremity of the seepage face. Again, due

343

to the fine spatial and temporal resolutions required to measure and model these features, diurnal transients in the unsaturated zone are also beyond the scope of this study, as are simulations at lower values of the dispersivities.

Three-dimensional river-aquifer models The preceding vertical section models provided insight into flow processes normal to the river bank. However, as was seen from the two-dimensional simulations in plan projection, the spatial curvature of the river bank was sufficient to induce convergence of regional groundwater flows in the vicinity of the study site. Three-dimensional models of flow and transport are therefore required in order to simulate the combined influences of topography and density-driven flow. In this section, several three-dimensional models are discussed. First, steady simulations of the saline wedge are considered for three different river salinities (with identical mean recharge rates to the land surface). Second, the river salinity is allowed to vary on a seasonal basis (in the presence of constant recharge), and finally both the river salinity and recharge rate are varied on a seasonal basis. Before commencing this discussion, we note that dynamics occurring at diurnal timescales or shorter are neglected for the reasons outlined in the previous section. The bank slopes defined in the three-dimensional model vary spatially throughout the grid. This is because the elevations of the ground surface were generated by spatial interpolation from surface contour data and from spot river depth measurements. However, bank slopes near the intersection of the notional contaminant plume and the shore line were 5.6 and 1.0 for the intertidal and deeper bed zones, respectively. These values are broadly representative of the slopes throughout the three-dimensional model and compare well to the slopes specified in the earlier vertical section models.

Grid development A three-dimensional grid was developed from the previous two-dimensional plan and section models discussed. The three-dimensional saturated flow model contains 17 slices defining 16 layers. The number of slices was chosen by calculating the effect on steady head distributions and solute concentrations of increasing the net number of slices. By the time 17 slices were included in the model the net change of head and concentration was negligible. Each slice contains 4846 nodes, giving a total of 82,382 nodes and 152,208 triangular prism elements. The grid spatial triangulation is irregular, with more resolution in the vicinity of the contaminant plume (see Fig. 5). The layers are approximately evenly spaced in the vertical, although the Feflow ‘‘best adaptation to stratigraphic data’’ (BASD) scheme moves the layers dynamically to minimise solution error during time stepping. Hydraulic conductivity distribution The river was idealized as a three-dimensional zone of arbitrarily large hydraulic conductivity (K = 2160 m d1)

344

M.G. Trefry et al.

Figure 10 Fluid ages calculated along a transect through the study area for a steady density coupled simulation with recharge and two dispersivity models. The full simulated transect is not shown, 157 m has been truncated from the left side to focus on convective processes. Parts (a,b) show fluid age distributions for aL = 5.0 m and aT = 0.5 m, and parts (c,d) for aL = 1.25 m and aT = 0.125 m. Parts (b,d) are logarithmic transformations of parts (a,c), respectively.

and unit storativity (S = 1.0), i.e. representing a fully saturated and highly conductive medium. The base of the aquifer was set at 15 m AHD, with the river water level set to be a constant 0 m AHD. Elevations for the first (top) and second slices were set to mimic the general ground surface and river bed profiles used in the two-dimensional vertical slice simulations. The volume between the first and second slices and east of the river bank was set to the river K and S values, forming a wedge shape as shown in Fig. 11. The vertical and eastern boundary nodes of this zone were assigned first-type head conditions, to set the river stage height, and first type mass conditions to set the river salinity. The remainder of the aquifer was assigned the mean conductivity value of 12.1 m d1 and a porosity h = 0.3. Saline dispersivities were set to aL = 5.0 and aT = 0.5 m. The simulations were run in Feflow’s phreatic mode, i.e. with a free and moveable top surface. The remainder of

the problem parameters are discussed in the following sections.

Steady recharge and river salinity In this case, a uniform recharge rate equivalent to an annual total of 200 mm to the water table was specified over the land surface. Zero recharge was applied to the river surface. Three simulations were run for river saline concentrations at the maximum (35 g L1), median (15 g L1) and minimum (5 g L1) annual concentrations (denoted by Cmax). These simulations were run to quasi-steady state using the transient flow–transient transport method in Feflow; several centuries of simulation time were elapsed to ensure quasi-steady conditions. These simulations were used to define theoretical maximum and minimum bounds for the saline convection dynamics near the river bank.

Hypoaigic influences on groundwater flux to a seasonally saline river

345

Figure 11 River idealization in the three-dimensional simulations (cut-away view). The river is represented by a zone of high K value (darker, K = 2160 m d1) with the aquifer represented by a zone of lower K value (lighter, K = 12.1 m d1). River boundary conditions are set along the top face of the three-dimensional grid and vertically down the eastern boundary as far as the river zone extends.

Topographic influences on seepage Possible relationships between landform shape and groundwater flow were investigated. Smith and Turner (2001) showed that periodic spatial meanders in river banks are sufficient to define spatially periodic variations in groundwater seepage flux, even in situations of saline-driven convection. In the present situation there are no meanders; however the high curvature of the river bank in the horizontal plane just south of Canning Bridge was seen to generate a focussing effect on the groundwater flows. Calculations of maximum seepage velocity were made along the river bank from the steady three-dimensional solutions and are summarized in Fig. 12 as functions of along-bank arc length using the following conventions and definitions specific to the Canning River site: (i) The river bank is parameterised as the locus of points (x(y), y), where y is the northing coordinate and x(y) is the easting coordinate. There are N such points along the discretised bank. (ii) The land lies below the x(y) function value and the river above. P (iii) Arc length is defined as si ¼ ij¼2 jðxðy j Þ; y j Þ ðxðy j1 Þ; y j1 Þj for i P 2, . . ., N. (iv) The bank curvature and hypercurvature quantities are d2 x d3 x defined as dy 2 and dy 3 , respectively, evaluated at each of the arc lengths si. The data points for each simulation were gathered from the numerical output for the top (phreatic) slice at locations 1 m landward from the shoreline. The raw seepage data contain peak features that appear to be generated by sudden changes in direction of the discretised shoreline. If these features are reduced by filtering, more regular trends emerge. In Fig. 12, such peaks have been removed from the data sets and the remainder smoothed by an exponential kernel. Broadly speaking, for each of the three simulations, the seepage velocities are reasonably constant and below 0.2 m d1 for arc lengths less than 250 m, i.e. along the

Figure 12 Variation of the magnitude of groundwater velocity at the river bank as a function of bank arc length, measured from south to north (top plate). The plume intersects the bank between the two vertical dashed lines. Results are shown for three steady river salinities. Also shown are the bank curvature and hypercurvature functions. In the lower plate, the groundwater velocity magnitude distribution is shown for the high salinity simulation, in three-dimensional view. The zone of maximum velocity magnitude is circled on the top face, corresponding to maximum seepage velocity.

southern reach of the shoreline. The velocities then trend upward to 0.25–0.3 m d1 by 330 m arc length, and then rapidly decline to values near zero further northwards along the shore. The topographic curvature of the river bank is also plotted as a function of bank arc length in Fig. 12. In physical terms, the maximum seepage velocity occurs at the beginning of the Canning Bridge salient where bank curvature changes rapidly from near zero to 0.02 m1. Then, as the salient develops the seepage velocity declines to small values in the vicinity of the bridge. This is fully consistent with the flow field diagrams in Figs. 4 and 5. Coincidentally, the hydrocarbon plume reaches the river bank at the same location as the peak seepage velocities are

346

M.G. Trefry et al.

predicted to be (see Fig. 5 and vertical bars in Fig. 12). As a side note, there is statistical evidence of a correlation between bank curvature and local seepage velocity; calculations on the data presented in Fig. 12 yield a linear correlation coefficient of 0.6 with the presence of a positive correlation strongly supported at the 95% confidence level by t-test. The correlation between seepage and bank hypercurvature, i.e. the third spatial derivative, is lower at 0.4, although there is clearly a strong positive hypercurvature peak at the same arc length as the seepage peak. A second t-test indicated that the hypercurvature–seepage velocity relationship was also significant. There is good statistical evidence for a purely topographic shoreline control on seepage velocity; however there are influences from many scales on seepage velocity, including regional boundary conditions affecting flow, and local and regional topographic curvatures.

Seepage and river salinity Also notable in Fig. 12 is the definite trend in seepage velocities for different river salinities. As salinities decline, seepage velocities also decline. This is rationalized in terms of the saline wedge behaviour under differing river salinities. At steady state, wedge penetration distances against regional flow are strongly affected by river saline load: high salinities result in long penetration distances which, in turn, act to focus regional groundwater flux into shorter seepage faces along the shoreline, hence boosting local seepage velocities. The reverse is true, with less saline river conditions yielding shorter wedges and long seepage faces, and hence smaller seepage velocities. Fig. 13 shows the variation of wedge penetration distance along the river bank arc, as calculated for a small set of control planes by taking the difference of the intersection locations of the 2/3 Cmax contour with the river bed and aquifer base. For high river salinity (equivalent concentration Cmax = 35 g L1) the wedge penetrates well beyond the shoreline (positive penetration) while for low salinity (equivalent concentration Cmax = 5 g L1) the regional discharge to the river forces the dilute wedge back underneath the river bed (negative penetration). At the median salinity (equivalent concentration Cmax = 15 g L1) the wedge penetrates almost 6 m under the river bank in the vicinity of the contaminant plume. Statistical comparisons show a positive linear correlation (r = 0.50) between bank curvature and wedge penetration distance, with an 83% level of confidence. However there was no evidence for correlation with seepage velocity. This is counter to expectations, but the number of control sections defined in Fig. 13 is small, so it is possible that a significant correlation with seepage velocity does exist and would be uncovered by a more exhaustive set of control sections. As mentioned above, these results do not represent actual dynamics at the site, since they explicitly ignore seasonal and other effects. It is expected that the seasonal dynamics will be represented by a much smaller range of seepage velocities and wedge penetration distances than outlined above, since the wedges respond at timescales much longer than the normal seasonal periods. Nevertheless, the steady results discussed here have focussed atten-

Figure 13 Variation of saline wedge penetration distance as a function of bank arc length, measured from south to north. Results are shown for three steady river salinities. The top plate shows the location of plane sections (with associated bank arc lengths) used for measuring wedge penetrations. The plume intersects the bank between the two vertical dashed lines.

tion on some prime spatial influences on seepage velocity in the hypoaigic zone.

Seasonal recharge and river salinity Allowing the recharge and river salinity to vary on seasonal time scales provides a better representation of the site dynamics. In the absence of field measurements of recharge and local river salinity, data of Fig. 2 were used to establish cyclic recharge and salinity variations which were applied to the model grid. The model was then run for many annual cycles to yield a quasi-stationary state (i.e. a periodic state of the model in which predictions at a given time of the year are effectively constant from year to year despite the seasonal variations within each year) after which time monthly snapshots of the system were saved. The snapshots give a means of estimating dynamical effects in the hypoaigic zone throughout the annual cycle. In fact, three such models were used:

Hypoaigic influences on groundwater flux to a seasonally saline river

347

• Model A: both the recharge rate and the river salinity were held constant at mean annual values throughout the year. The annual mean values were assumed to be 24 and 0.6 mm d1 for the salinity and recharge rate, respectively. • Model B: the salinity varied according to Fig. 2 and the recharge was held constant at the mean annual value. • Model C: both the salinity and recharge varied according to Fig. 2. Before comparing the simulation results from these three models, we first consider the convergence of the model results with respect to time.

Stationarity Stationarity is the consistency of model predictions from year to year in the presence of repeated intra-annual cycles. The approach to stationarity of three-dimensional convective models is complex, essentially because of the widely different characteristic time scales of the seasonal, buoyancy and diffusive processes. Hence some local measures may converge to a stationary state at a faster rate than other measures. For this reason, a range of indicators was used to assess stationarity. First, the total saline mass load in the model was calculated at the ends of annual cycles. Second, local salinities were calculated at the ends of annual cycles at a total of 165 monitoring points distributed regularly at five elevations in the profile near the intersection of the contaminant plume and the river bank. Finally, local heads were calculated at the same 165 points. These measures were then collected at annual intervals over many annual cycles, providing information on the convergence of the measures to stationary limits. Typically in these convective simulations, heads attained quasi-stationary states faster than concentrations, but then their rates of convergence diminished to be congruent with the concentration stabilization measures. In this section we show how the stabilization curves can be modelled empirically to provide estimates of suitable stabilisation times for production simulations. It was found that inadequate stabilization times led to quite different predictions of intra-annual cycles for fluxes and wedge characteristics, even if successive annual measures were quite similar. The approach to stationarity is easiest to assess when boundary conditions are steady, i.e. for Model A. Starting with a steady (high salinity) solution, the mean annual recharge and salinity values were applied to the model (with constant upgradient head boundary set at 1.0 m AHD and the river elevation at 0.0 m AHD) and the relaxation towards a steady state was monitored. Fig. 14a shows the relaxation of the total mass within the grid together with the local head (b) and saline concentration (c) measured at a location just above the aquifer base 10 m inland from the river bank along the flow line toward the contaminated site. This location is within the penetrating wedge and so will sample changes in wedge concentration during the approach to steady state. Results for this location are given for all three models. Let us confine attention to the results for Model A (filled squares in Fig. 14). The total mass reduces rapidly at first as excess saline mass exits the river through the imposed first

Figure 14 Stabilization of integrated (total mass, (a)) and point (head (b) and concentration (c)) measures for the threedimensional models A, B and C. The concentrations at 10 m AHD have been divided by 10 for graphical convenience.

type condition. The mass stabilization is well represented by a short-lived exponential form plus a longer, approximately linear decline. This empirical function is drawn as a dashed line in the figure: the exponential part is 99.99% relaxed after 42 years, as listed in Table 1. Thereafter the rate of mass loss stabilizes to a slow but steady value over a century timescale or longer. This stage represents a diffusive/dispersive process where fresh groundwater slowly dilutes the saline wedge, stripping mass and exiting through the seepage face. Also, a legacy saline zone directly under the river is slowly decaying via diffusion. This shows nicely

348

M.G. Trefry et al.

Table 1 Relaxation times estimated from the stabilization simulations for the three models of Section ‘‘Seasonal recharge and river salinity’’ Elevation

Model A

Model B

Model C

0 m AHD

Head: 22 years Conc: 15 years

Head: 24 years Conc: 10 years

Head: 37 years Conc: 12 years

5.5 m AHD

Head: 21 years Conc: 16 years

Head: 20 years Conc: 10 years

Head: 44 years Conc: 10 years

10 m AHD

Head: 24 years Conc: 30 years

Head: 21 years Conc: 35 years

Head: 45 years Conc: 55 years

Total Mass

42 years

15 years

20 years

The times are calculated from a simple negative exponential fit, assuming 99.99% decay.

how episodic saline events in surface waters can lead to long-lived saline distributions in underlying aquifers where the rate of advective stripping is very low. Similarly, Fig. 14b and Table 1 show that the local head stabilizes very rapidly (99.99% relaxation after 22 years, 21 years and 24 years at high, intermediate and low elevations in the aquifer, respectively), implying that long-term dynamics under saline rivers are driven essentially by saline processes rather than landscape hydrology. The local concentration stabilization, too, supports this conclusion with 99.99% relaxation times of 15 years, 16 years and 30 years at high, intermediate and low elevations in the aquifer, respectively. In summary, near the water table the solution relaxes within a decade or so, and then enters a phase of very slow equilibration. However, at depth in the aquifer within the saline wedge, the initial relaxation is significantly slower and the long-term equilibration has a slightly higher magnitude of slope. At a larger scale, we can expect that rates of model stabilization will depend on local advective conditions throughout the grid. Where groundwater fluxes are low, and hence dynamics are purely diffusion-limited, any saline distributions will take exceedingly long time to reach steady state. In the present situation, this is likely to be the case under the peninsula close to the Canning Bridge. For practical reasons, we choose to focus on saline dynamics in the contaminated zone. Guided by the above results for the stepchange Model A, the transient models (Models B and C) were run for 60 years commencing from the same high salinity initial condition. Such a limited stabilization period may lead to potential over-estimation of stationary salinities directly under the river (since the dynamics there will be diffusionlimited). However the advection-dominated zone near the seepage face and wedge will approach a quasi-stationary state much more rapidly, and it is this zone that we are most interested in. As an aside, we note that it is likely that many (saline) riverine aquifers are in non-stationary conditions since climatic variations can occur at timescales of the order of decades, often well within the longer stabilization times indicated by Fig. 14. Fig. 14 also shows the convergence to stationary conditions of indicators from Models B (filled triangles) and C (open circles), with stabilization times listed in Table 1. Note that since Models A, B and C each correspond to different physical problems, it is not necessary for their predic-

tions to approach one another in the limit of increasing stabilization time. The data presented in Fig. 14 for Models B and C were collected at annual boundaries, so that they occurred at identical times in the annual cycle. Truly stationary conditions would result in constant values for the various measures as a function of stabilization time. It is clear from Fig. 14 that the solutions for both Models B and C are reasonably stable and stationary even through they have only been integrated over 60 annual cycles. Shorter model stabilization times may lead to inaccuracy in predicted seasonal trends, even at the high-flux contaminated site. In all three Models the predicted stabilization times at the lowest elevation (within the wedge) are the longest, and certainly longer than the thirty year period recommended for fresh groundwater systems characterisation (see Von Asmuth and Knotters, 2004). Seasonal dynamics in the contaminated zone Models A, B and C were run for 60 years of simulated time. For Models B and C, seasonal variations were captured by storing intra-annual solutions at 30 day intervals for the last year of simulation. Large spatial variations were seen in local seepage rates and wedge shapes along the river bank, and large temporal variations were seen throughout the annual cycle. In particular, as suggested by the earlier steady results, an extensive saline wedge persists under the bridge throughout the year, whereas further south the wedge waxes and wanes according to seasonal influences. Fig. 15 provides shaded contours of the saline concentration distribution in June and December (simulation time). The cut-away sections in the figure clearly show how the saline wedge penetrates far under the Canning Bridge river bank where groundwater fluxes are least. There is significantly less penetration under the river bank to the south of the bridge. It is also apparent how the variations in river salinity happen much faster than the deeper saline distributions can respond. For example, Fig. 15a shows the lowsalinity river condition at June 30 of the annual cycle as a light grey shade, yet just beneath the river there is a zone of darker shading, indicating a layer of high salinity. The situation is reversed for the December 31 case shown in Fig. 15b. In both plots, it is clear that the wedge penetrates much further into the aquifer underneath the Canning Bridge salient than is seen normal to the bank further to the south. The subtle interplays between the variations of

Hypoaigic influences on groundwater flux to a seasonally saline river

Figure 15 Simulated seasonal concentration distributions at June 30 (a) and December 31 (b). The view looking from the north-west to the south-east; the west-east section line passes just south of the point of intersection of the plume with the river bank.

river salinity and aquifer recharge, together with the lagged and slow responses of the aquifer saline wedges at locations along the riverbank, are key to understanding the potential for degradation of advecting contaminants near surface water bodies. Whilst the variations all along the bank are of interest, focus is now directed to the seasonal fluctuations near the

Figure 16 Simulated seasonal variation of wedge penetration length (top curves) and maximum seepage rates (bottom curves) for the three-dimensional models A (solid lines), B (circles, dotted) and C (squares, dashed).

349

intersection of the subject contaminant plume with the river. The central flow line of the predicted plume (Fig. 5) was approximated by a linear section of length 215 m; the centre of the section was close to the river shoreline. Groundwater flow rates and concentration distributions were calculated at monthly intervals along the vertical slice defined by the 215 m section. Fig. 16 shows the seasonal variation of wedge penetration along the section (defined by the intersection of the salinity 20 contour with the aquifer base) and also the peak groundwater speed for the three Models. The groundwater speed is maximal at the seepage face at the river bank, so the values in Fig. 16 represent the predicted rate (speed) of seepage immediately beneath the sediment–water interface. Fig. 16 shows that the wedge penetration distance ranges between 7.3 and 8.3 m (6.4% fluctuation) for Model B, and between 6.4 and 8.5 m (14% fluctuation) for Model C. This is an important result; Model B only incorporates saline variability whereas Model C incorporates both saline and recharge variability. Since Model C demonstrates higher variation in wedge penetration than does Model B, it is clear that variations in seasonal recharge rates are a dominant control on wedge persistence. Interestingly, the peak wedge penetrations are approximately coincident with the time of greatest recharge, i.e. July–August. This indicates the presence of significant lags in the relaxation of the concentration distribution in response to peak river salinities which occur in February–March. Of the two seasonal Models, only Model C attains the wedge penetration distance of the mean-condition Model A. The results for the seepage rates in Fig. 16 also deserve some study. Here the rates are highest during the summer period and lowest in late winter and early spring, which seems to be in direct opposition to the recharge signal. The reason for this is a subtle interplay between the river salinity and the accessible seepage face. We must draw a distinction between the groundwater seepage flux, i.e. the net amount of seepage measured in terms of volume per unit time, and groundwater seepage rate, i.e. the speed at which groundwater moves to the seepage face. Regardless of recharge amount, when the river salinity is increased, the accessible length of river bank in the models through which the groundwater flux can exit is reduced. This is because high-density river water can easily displace fresh-to-brackish groundwater by simple buoyancy processes, and because the shoreline location is fixed (not movable) in the model. Thus a given flux of groundwater must exit the aquifer over a short interval of bank slope, leading to high seepage rates. When river salinities decline in the winter/spring period, saline densities also decline and a longer interval of the bank is accessible to groundwater flux, which in turn yields a reduction in seepage rate. The sensitivity of the modelled seepage rate to the river salinity is demonstrated by the sharp increase in seepage rate at day 365. At this time the salinity boundary condition re-establishes a dense overtop to the aquifer as the wedge retreats, impelling the incident groundwater flux to a narrow seepage face. With this in mind, we can revisit the wedge penetration effect noted above. Since the maximum recharge rate coincides with the minimum river salinity, the peak wedge penetration correlates with the time of maximum accessible seepage face. We

350 can rationalize the situation as follows: when the river salinity is least, there is also the least resistance to groundwater seepage at the river bank. Groundwater flux is then attracted to the permissive seepage face, thereby reducing the local pressure on the deep wedge, which responds by expanding. As the river salinity increases again, local pressures on the deep wedge increase and the wedge contracts. We now consider detailed maps of the wedge geometry through the annual cycle, as presented in Fig. 17 for Model C. The annual cycle commences at day 0 with a dense layer of saline river water (salinity >25 g L1) overlying a re-entrant zone of moderate density fluid (15 < salinity < 25) which itself overlies the saline wedge at the bottom of the aquifer. Up until day 120, the overtopping saline layer strengthens and the less-saline re-entrant zone diminishes.

M.G. Trefry et al. Further to the right, away from the river bank, there is evidence of unstable fingering of the saline distribution as indicated by the variation of the salinity 30 contour. After day 120 the river salinity drops rapidly to brackish levels, which eliminates the dense overlying layer and tends to dilute the top surface of the brackish wedge. By day 300, the river salinity is increasing fast, and a new dense layer is forming just beneath the river bed. Crucially, this new dense layer forms faster than the less dense layer below it can re-equilibrate. The result is a sequence of trapped zones of reduced salinity, perhaps isolated from each other by dense saline fingers falling from above. In the day 300 (late October) plot, there are two such zones visible as elongate closed contours of salinity value 15. This complex stratification in the salinity distribution is reminiscent of the September 2001 field measurements

Figure 17 Simulated seasonal variation of wedge distribution for Model C. The wedge undergoes significant and complex changes throughout the seasonal cycle. The river lies to the right of the plots, and the land to the left, as indicated by the overbars. Contours are numbered according to salinity.

Hypoaigic influences on groundwater flux to a seasonally saline river plotted in Fig. 8. Although the qualitative resemblance is strong, the field data do show very low salinities in the reentrant zone, in contrast to Fig. 17 where the re-entrant zone has salinities everywhere higher than 10. This is likely an artefact of the numerical simulation. For practical reasons, finite element dimensions were constrained to be larger than 1 m, and often much larger. This, in turn, constrains the feasible size of dispersivities to be larger than 1 m, resulting in potentially exaggerated smearing of the saline distributions in the aquifer. Greater resolution in the wedge contours would be obtained by re-gridding the model and using lower dispersivities, but this is computationally expensive and beyond the scope of this study. In any event, these three-dimensional results have shown that seasonal effects are sufficient to account for saline distributions of the kind measured at the site, especially where dispersion is low. This supports the earlier time scale analysis which identified seasonal influences with spatial scales larger than 1 m and tidal influences with smaller spatial scales in the aquifer porous medium. At the other end of the concentration scale, the salinity five contour is relatively immobile during the annual cycle in Fig. 17. This contour represents the most dilute fringe of the saline distribution and is particularly sensitive to the assumed dispersivities. Nevertheless, close inspection reveals that the upper part of the contour moves towards the river in the winter flush and in the opposite direction as the effects of seasonal recharge on regional groundwater flux diminish and river salinity rises. The toe location of the salinity five contour shows a similar qualitative trend, although the trend is lagged with respect to the upper contour location movement.

Discussion A variety of steady and transient three-dimensional models were used to study topographic, climatic and estuarine controls on seepage rates at the contaminated site. In the earlier steady simulations, it was seen how the contaminant plume reaches the river bank at a local high in seepage rate. This local extremum in seepage rate is controlled by the long-shore curvature of the river bank, i.e. it is a topographic effect. Seepage rates on the promontory near the bridge are much lower, and further south on the almost linear reach of the bank they take on intermediate values. Some statistical correlations between long-shore distance and seepage rates were apparent in the simulation results, in line with earlier work by Smith and Turner (2001). These kinds of correlations are difficult to establish using twodimensional (plan and section) models alone, because they ignore the kinds of couplings between vertical and horizontal influences that are expressed so clearly in three dimensions with the strong horizontal variations in wedge penetration near Canning Bridge. Groundwater fluxes and seepage rates are inextricably linked with saline distributions and wedges arising from river saline variations. The seasonal simulations showed how the dense river water controlled seepage rates by narrowing the seepage face, which in turn led to increased pressure on the deep saline wedge. This gave rise to the counter-intuitive observations that

351

1. times of low river salinity were likely to correlate with expansions of the saline wedge at the aquifer base by approximately 15%; 2. times of high recharge were likely to correlate with reductions in seepage rate by up to 50%. The central factor in both these observations is the issue of delay: the saline wedge responds to changes in the hydraulic environment in a complex manner. No single equilibration/stationarity time applies to all parts of the wedge, so we face a situation where some indicators of wedge state (e.g. toe location) may vary out of phase with the dominant driving variables (e.g. river salinity and regional groundwater flux), i.e. the observations above are unlikely to be causal. In this respect, it becomes easier to rationalize the observations above – the delays in response of some wedge properties to changing river salinities mean that we must understand the cyclic dynamics of the wedge in order to draw meaningful correlations with important driving variables. Because of delayed responses, phenomena that appear implausible for steady systems may be natural occurrences in cyclic systems. Underlying the two observations made above is the assumed anti-correlation between the seasonal river salinity variation and the seasonal recharge rate (see Fig. 2). In the absence of field data to constrain these seasonal variables and their correlation, the variations in Fig. 2 were based on the mean local monthly rainfall data. Thus the observations presented are simple idealizations, yet they do describe a rich variety of effects that may well operate in similar hypoaigic circumstances. Seasonal recharge is seen to be a dominant control on wedge penetration, and seasonal river salinity is a dominant control on groundwater seepage rates. In general, seasonal influences appear to be sufficient to induce re-entrant salinity distributions of the kind displayed in Fig. 8, especially if local dispersion is small. We also remark that the estimated groundwater seepage rates of between 0.28 and 0.5 m d1 at the contaminated river bank are the best estimates available. As a word of caution, we repeat that these results have been generated using a fixed boundary idealization for the river water–river bank location. This idealization constrains the seepage face to be at a fixed location in space, regardless of upgradient fluxes (which depend on season). The model would be improved by incorporating a moving seepage boundary in future work.

Conclusions Even though the Canning River site has a relatively simple geological stratigraphy and the aquifer material is reasonably uniform and well-sorted, a variety of effects compound to generate a complex hypoaigic environment. Firstly, the Canning River has a strong seasonal variation in salinity, which induces significant and periodic convective flows in the adjoining aquifer. Secondly, the river bank describes a complex curve in horizontal projection, the curvature of which induces important variations in local groundwater flux along the bank. Coincidentally, the contaminant plume advects in a south-easterly direction towards the river bank, intercepting the river bank at a point of maximum ground-

352 water flux and maximum bank curvature. This has important ramifications for natural attenuation of the advecting dissolved species. For example, dissolved species will encounter a range of fluid types and flow speeds before exiting the aquifer at the river seepage face. Fluid types encountered include the mature groundwater advecting from further inland, the relatively oxygenated recharge water, and potentially much older fluids, originating from saline river waters, dispersing into the dissolved contaminants stream from the saline wedge. As this chain of mixing types progresses, fluid velocities increase to maxima at the seepage face where fresh river waters are met. Penetration of the saline wedge under the river bank is controlled by the relative densities of the river water and groundwater, and by the water table gradient, which is a function of the regional flux and K, discharging to the seepage face. Measurements showed that the wedge extended more than 15 m inland from the river bank in September 2001. Two-dimensional simulations in vertical section were able to mimic this penetration distance by using relatively low values of dispersivity. On the other hand, three-dimensional simulations produced smaller penetration distances (10 m). The three-dimensional simulations were restricted to larger values of dispersivity, and also incorporated effects of groundwater flow focussing due to topographic bank curvature. It is likely that the three-dimensional simulations could be adjusted to yield closer cosmetic agreement to the observed wedge penetrations, e.g. by reducing the hydraulic conductivity and dispersivity values. This has not been attempted due to the scarcity of independent field measurements to support the changes, and since the important features of the seasonal hypoaigic processes are already apparent in the three-dimensional solutions. Another important effect is provided by the annual variation of recharge. Recharge modulates the groundwater flux, resulting in strong periodic fluctuations in the saline distributions arising from the interaction of the river water with the aquifer. As the saline distributions fluctuate throughout the year, they induce significant variations in groundwater seepage rates through the river bank; the variations have amplitudes up to 50% with late winter and spring having the lowest seepage rates and summer and early autumn the highest. Note that this is the reverse of the seepage flux trends, when greatest fluxes are in late winter and spring and lowest fluxes in summer and autumn. So we have a situation where greatest seepage rates accompany minimum groundwater fluxes, and lowest seepage rates accompany maximum groundwater fluxes. The dynamic range between these variations is governed largely by the location of the contaminant plume within the region of maximum long-shore bank curvature.

Acknowledgements The authors gratefully acknowledge the work of P.L. Bjerg, S.J. Fisher, N. Innes, C.D. Johnston, J.L. Rayner, B.S. Robertson, S. Westbrook and R.J. Woodbury in the collection, preparation and analysis of field data referenced in this manuscript, and B. Robson, D. Rassam and A.J. Smith for

M.G. Trefry et al. useful suggestions on an earlier version of this manuscript. The suggestions of an anonymous reviewer were particularly helpful in improving this work. The WA Dept for Planning and Infrastructure is thanked for supplying the Perth Gauging Station tidal data, as is the Bureau of Meteorology for supplying Perth rainfall data used in this work. This work was partially funded by BP and Shell.

References Ataie-Ashtiani, B., Volker, R.E., Lockington, D.A., 1999. Tidal effects of seawater intrusion in unconfined aquifers. J. Hydrol. 216, 17–31. Becker, M.W., Georgian, T., Ambrose, H., Siniscalchi, J., Fredrick, K., 2004. Estimating flow and flux of groundwater discharge using water temperature and velocity. J. Hydrol. 296, 221–233. C ¸ elik, M., 2002. Water quality assessment and the investigation of the relationship between the River Delice and the aquifer systems in the vicinity of Yerko ¨y (Yozgat, Turkey). Environ. Geol. 42, 690–700. Clavero, V., Izquierdo, J.J., Fernandez, J.A., Niell, F.X., 1999. Influence of bacterial density on the exchange of phosphate between sediment and overlying water. Hydrobiologia 392, 55– 63. Dai, Z., Samper, J., 2006. Inverse modelling of water flow and multicomponent reactive transport in coastal aquifers. J. Hydrol. 327, 447–461. Davidson, W.A., 1995. Hydrogeology and groundwater resources of the Perth region, Western Australia. Western Australia Geological Survey, Bulletin 142. Diersch, H.-J., 2002. Feflow Reference Manual. WASY, Berlin. Diersch, H.-J., 2004. Feflow 5.1 – Finite Element Subsurface Flow and Transport Simulation System. WASY, Berlin. Diersch, H.-J.G., Kolditz, O., 2002. Variable-density flow and transport in porous media: approaches and challenges. Adv. Water Resour. 25, 199–944. Duncan, T., Shaw, T.J., 2003. The mobility of rare earth elements and redox sensitive elements in the groundwater/seawater mixing zone of a shallow coastal aquifer. Aquat. Geochem. 9, 233–255. Fetter, C.W., 1994. Applied Hydrogeology, third ed. Macmillan, New York. Fofonoff, N.P., Millard Jr., R.C., 1983. Algorithms for computation of fundamental properties of seawater. UNESCO Technical Papers in Marine Science, 44, UNESCO. Gelhar, L.W., Welty, C., Rehfeldt, K.R., 1992. A critical review of data on field-scale dispersion in aquifers. Water Resour. Res. 28, 1955–1974. Golosov, S.D., Ignatieva, N.V., 1999. Hydrothermodynamic features of mass exchange across the sediment–water interface in shallow lakes. Hydrobiologia 409, 153–157. Glover, R.E., 1959. The pattern of fresh-water flow in a coastal aquifer. J. Geophys. Res. 64, 457–459. Goode, D.J., 1996. Direct simulation of groundwater age. Water Resour. Res. 32, 289–296. Hamilton, D.P., Chan, T., Robb, M.S., Pattiaratchi, C.B., Herzfeld, M., 2001. The hydrology of the upper Swan River Estuary with focus on an artificial destratification trial. Hydrol. Process. 15, 2465–2480. Helena, B., Pardo, R., Vega, M., Barrado, E., Fernandez, J.M., Fernandez, L., 2000. Temporal evolution of groundwater composition in an alluvial aquifer (Pisuerga River, Spain) by principal component analysis. Wat. Res. 134, 807–816. Huettel, M., Røy, H., Precht, E., Ehrenhauss, S., 2003. Hydrodynamical impact on biogeochemical processes in aquatic sediments. Hydrobiologia 494, 231–236.

Hypoaigic influences on groundwater flux to a seasonally saline river Hughes, C.E., Binning, P., Willgoose, G.R., 1998. Characterisation of the hydrology of an estuarine wetland. J. Hydrol. 211, 34–49. Jacob, C.E., 1950. Flow of ground water. In: Rouse, H. (Ed.), Engineering Hydraulics. Wiley, New York, pp. 321–386. Kim, J.-H., Lee, J., Cheong, T.-J., Kim, R.-H., Koh, D.-C., Ryu, J.S., Chang, H.-W., 2005. Use of time series analysis for the identification of tidal effect on groundwater in the coastal area of Kimje, Korea. J. Hydrol. 300, 188–198. Kruseman, G.P., de Ridder, N.A., 1994. Analysis and Evaluation of Pumping Test Data, 2nd ed. Publication 47, International Institute for Land Reclamation and Improvement, Wageningen, The Netherlands. Lambs, L., 2004. Interactions between groundwater and surface water at river banks and the confluence of rivers. J. Hydrol. 288, 312–326. Lebbe, L., 1999. Parameter identification in fresh-saltwater flow based on borehole resistivities and freshwater head data. Adv. Water. Resour. 22, 791–806. Linderfelt, W.R., Turner, J.V., 2001. Interaction between shallow groundwater, saline surface water and nutrient discharge in a seasonal estuary: the Swan–Canning system. Hydrol. Process. 15, 2631–2653. Mao, X., Enot, P., Barry, D.A., Li, L., Binley, A., Jeng, D.-S., 2006. Tidal influence on behaviour of a coastal aquifer adjacent to a low-relief estuary. J. Hydrol. 327, 110–127. Moore, W.S., 1999. The subterranean estuary: a reaction zone of ground water and sea water. Marine Chem. 65, 111–125. Naegeli, M.W., Uehlinger, U., 1997. Contribution of the hyporheic zone to ecosystem metabolism in a prealpine gravel-bed river. J.N. Am. Benthol. Soc. 16, 794–804. Neilson-Welch, L., Smith, L., 2001. Saline water intrusion adjacent to the Fraser River, Richmond, British Columbia. Can. Geotech. J. 38, 67–82. Oswald, S.E., Kinzelbach, W., 2004. Three-dimensional physical benchmark experiments to test variable-density flow models. J. Hydrol. 290, 22–42. Otero, J., Dontcheva, L.A., Johnston, H., Worthing, R.A., Kurganov, A., Petrova, G., Doering, C.R., 2004. High-Rayleigh-number convection in a fluid-saturated porous layer. J. Fluid Mech. 500, 263–281. Pitman, A.J., Narisma, G.T., Pielke, R.A., Holbrook, N.J., 2004. Impact of land cover change on the climate of southwest Western Australia. J. Geophys. Res.–Atmos. 109. Article No. D18109. PSS, 1978. Practical salinity scale. IEEE J. Oceanic Eng. OE-5, 14. Simmons, C.T., Narayan, K.A., 1997. Mixed convection processes below a saline disposal basin. J. Hydrol. 194, 263–285. Simpson, S.L., Maher, E.J., Jolley, D.F., 2004. Processes controlling metal transport and retention as metal-contaminated groundwaters efflux through estuarine sediments. Chemosphere 56, 821–831. Smith, A.J., 1999. Application of a tidal method for estimating aquifer diffusivity. Technical Report 13/99, CSIRO Land and Water, Australia. Smith, A.J., 2004. Mixed convection and density-dependent seawater circulation in coastal aquifers. Water Resour. Res. 40. Article W08309. Smith, A.J., Turner, J.V., 2001. Surface water–groundwater interaction and nutrient discharge. Hydrol. Process. 15, 2595–2616. Stanford, J.A., Ward, J.V., 1988. The hyporheic habitat of river ecosystems. Nature 335, 64–66.

353

Taniguchi, M., Turner, J.V., Smith, A., 2003. Evaluation of groundwater discharge rates from subsurface temperature in Cockburn Sound, Western Australia. Biogeochemistry 66, 111–124. Teo, H.T., Jeng, D.S., Seymour, B.R., Barry, D.A., Li, L., 2003. A new analytical solution for water table fluctuations in coastal aquifers with sloping beaches. Adv. Water Resour. 26, 1239– 1247. Thompson, P.A., Waite, A.M., McMahon, K., 2003. Dynamics of a cyanobacterial bloom in a hypereutrophic, stratified weir pool. Mar. Freshwater Res. 54, 27–37. Thorenz, C., Kosakowski, G., Kolditz, O., Berkowitz, B., 2002. An experimental and numerical investigation of saltwater movement in coupled saturated-partially saturated systems. Water Resour. Res. 38. Article 1069. Townley, L.R., Trefry, M.G., 2000. Surface water–groundwater interaction near shallow circular lakes: flow geometry in three dimensions. Water Resour. Res. 36, 935–949. Trefry, M.G., Bekele, E., 2004. Structural characterization of an island aquifer via tidal methods. Water Resour. Res. 40. Article W01505. Ullman, W.J., Chang, B., Miller, D.C., Madsen, J.A., 2003. Groundwater mixing, nutrient diagenesis, and discharges across a sandy beachface, Cape Henlopen, Delaware (USA). Estuar. Coast. Shelf S. 57, 539–552. Volker, R.E., Rushton, K.R., 1982. An assessment of the importance of some parameters for seawater intrusion in aquifers and comparison of dispersive and sharp-interface modelling approaches. J. Hydrol. 56, 239–250. Volker, R.E., Zhang, Q., Lockington, D.A., 2002. Numerical modelling of contaminant transport in coastal aquifers. Math. Comput. Simulat. 59, 35–44. Von Asmuth, J.R., Knotters, M., 2004. Characterising groundwater dynamics based on a system identification approach. J. Hydrol. 296, 118–134. Water and Rivers Commission, 1997. Perth Groundwater Atlas, ISBN 0 7309 7345 X. Webster, I.T., 2003. Wave enhancement of diffusivities within surficial sediments. Environ. Fluid Mech. 3, 269–288. Westbrook, S.J., Rayner, J.L., Davis, G.B., Clement, T.P., Bjerg, P.L., Fisher, S.J., 2005. Interaction between shallow groundwater, saline surface water and contaminant discharge at a seasonally- and tidally-forced estuarine boundary. J. Hydrol. 302, 255–269. Winde, F., van der Walt, I.J., 2004. The significance of groundwater-stream interactions and fluctuating stream chemistry on waterborne uranium contamination of streams – a case study from a gold mining site in South Africa. J. Hydrol. 287, 178–196. Winter, T.C., 1999. Relation of streams, lakes, and wetlands to groundwater flow systems. Hydrogeol. J. 7, 28–45. Wooding, R.A., 1989. Convective regime of saline groundwater below a ‘dry’ lakebed, Technical Report 27, CSIRO Centre for Environmental Mechanics, Australia. Young, R.G., Huryn, A.D., 1996. Interannual variation in discharge controls ecosystem metabolism along a grassland river continuum. Can. J. Fish. Aquat. Sci. 53, 2199–2211. Zhang, Q., Volker, R.E., Lockington, D.A., 2001. Influence of seaward boundary condition on contaminant transport in unconfined coastal aquifers. J. Contam. Hydrol. 49, 201–215. Zhang, Q., Volker, R.E., Lockington, D.A., 2002. Experimental investigation of contaminant transport in coastal groundwater. Adv. Environ. Res. 6, 229–237.