On variable density surface water–groundwater interaction: A theoretical analysis of mixed convection in a stably-stratified fresh surface water – saline groundwater discharge zone

Journal of Hydrology (2006) 329, 390– 402

available at www.sciencedirect.com

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

On variable density surface water–groundwater interaction: A theoretical analysis of mixed convection in a stably-stratified fresh surface water – saline groundwater discharge zone Gudrun Massmann a,*, Craig Simmons b, Andrew Love a, James Ward b, Julianne James-Smith a a

South Australian Government, Department of Water, Land and Biodiversity Conservation, Knowledge and Information Division, Research and Innovation Group, GPO Box 2834, Adelaide, SA, Australia b Flinders University, School of Chemistry, Physics and Earth Sciences, GPO Box 2100, Adelaide, SA, Australia Received 9 June 2005; received in revised form 15 February 2006; accepted 21 February 2006

KEYWORDS

Summary Understanding the discharge behaviour of saline groundwater into fresh surface water can be critical for the effective management of water resources. While variable density flow has been studied intensely in a number of settings, the role it plays on the discharge behaviour of saline groundwater into freshwater streams is often neglected when calculating salt loads into a stream. The aim of this study was to determine what role variable-density flow behaviour plays in surface water/groundwater interaction in a stably-stratified fresh surface water/saline groundwater interface. The mixed convection ratio M, a measure of the ratio of density driven flow to advective driven flow, was defined for a matrix of one-dimensional numerical simulations that employed both varying hydraulic and density gradients. Vertical salt breakthrough into the surface water only occurred in the advection dominated cases (M < 1) and the salt flux into the surface water increased with increasing groundwater concentration until M reached a value of 1. Beyond this, when the flow was driven by the density difference between the two fluids (M > 1) vertical discharge of salt into surface water did not occur and the saltwater/freshwater interface migrated downwards with increasing density differences between the two fluids. This study therefore shows that there is a critical concentration difference that maximises salt loads to a surface water body and that a density-invariant approach to estimate the salt flux into the surface water (as the product of flow velocity determined

Variable density flow; Mixed convection ratio; Surface water/ groundwater interaction; River salinisation; Equivalent freshwater head

* Corresponding author. E-mail addresses: [email protected] (G. Massmann), [email protected] (C. Simmons), Love.Andrew@saugov. sa.gov.au (A. Love), [email protected] (J. Ward), [email protected] (J. James-Smith).



0022-1694/$ - see front matter c 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2006.02.024

On variable density surface water–groundwater interaction

391

through a potentiometric analysis and groundwater concentration) may be inadequate, especially where large density differences exist between the fresher surface water body and the underlying saline groundwater. The study is a purely theoretical approach and conclusions were drawn from simplified 1D simulations. Hence, further laboratory and modelling work is needed to confirm and test the plausibility of these findings for more realistic 2D and 3D cases. c 2006 Elsevier B.V. All rights reserved.

 Introduction

Variable density driven flow processes are important in a number of natural systems, where the density varies as a function of suspended solid content, temperature and pressure of a fluid. Fresh surface water bodies often overly more saline groundwater in many arid natural settings, for example along stretches of the River Murray in Australia, (Allison et al., 1990; Jolly et al., 1998) – what is termed a stably stratified hydrodynamic system. Hydrogeological field settings of studies that have previously emphasised the importance of variable density flow include, amongst others, upconing below wells (Diersch et al., 1984; Reilly and Goodmann, 1987), seawater intrusion in coastal regions (Volker and Rushton, 1982; Huyakorn et al., 1987; Cheng and Chen, 2001) and waste disposal in deep salt formations (Oldenburg and Pruess, 1995; Kolditz et al., 1998) or unstable flow phenomena where saltier and hence more dense fluid overlies fresher less dense groundwater as is often typical of, for example, saline disposal basins and salt lakes (Simmons and Narayan, 1997, 1998; Wooding et al., 1997a,b; Simmons et al., 2002) and dense contaminant plume migration (e.g., Liu and Dane, 1996; Oostrom et al., 1992a,b; Schincariol et al., 1997). Extensive reviews on the subject of variable density flow have, for example, been given by Holzbecher (1998), Simmons et al. (2001), Diersch and Kolditz (2002) and a recent article by Simmons (2005) highlights the current and future research challenges and needs in this research area. Understanding the complex interaction between groundwater and surface water is essential for the effective management of water resources (Sophocleus, 2002). Whilst numerous authors have examined groundwater discharge behaviour into rivers or lakes (see for example, reviews given in Winter, 1999; Sophocleus, 2002) the influence of density differences between a fresh surface water body and saline groundwater, which often exist, are typically not discussed. Most research has focussed on hydraulic (e.g., studies summarised in Winter, 1999; Schubert, 2002) or biogeochemical aspects (e.g., Jacobs et al., 1988; Bourg and Bertin, 1993; Dousson et al., 1997) of surface water/ groundwater interaction. Generally, for hydraulically connected surface water/groundwater systems, the flow is assumed to be controlled by the same mechanisms as those that occur in inter-aquifer leakage (Rushton and Tomlinson, 1979). The interactive flow is assumed to be a direct function of the head difference between the surface water body (river or lake) and aquifer and of the hydraulic conductivity of the semi-permeable river sediments only. In an extensive and detailed review on interactions between groundwater and surface water, Sophocleus (2002) points out that this simple approach of treating the flow between a river or lake

and aquifer is often too simplistic and that more appropriate models exist. However, in that recent review, there is no discussion on the possibility that density dependent flow behaviour may be an important consideration in surface/ groundwater interaction in streams and lakes, further suggesting that this subject has not previously been considered in any general conceptual or theoretical manner. Whilst a simplest approach might be to convert measured heads into equivalent freshwater heads, this approach is problematic particularly where there is vertical density stratification and vertical flows are of interest. Indeed, Lusczynski (1961) recognized that equivalent freshwater heads cannot be used to determine the vertical hydraulic gradient in an aquifer with water of non-uniform density. Thus, the use of equivalent freshwater head analyses together with a standard density-invariant version of Darcy’s Law is expected to lead to errors in estimation and prediction of surface water-groundwater interaction processes and the subsequent estimates of salt discharge to the surface water body. A system in which flow is driven by both forced (hydraulically driven) convection (also commonly called advection) as well as free (density driven) convection is a mixed convective flow regime. Mixed convective flow can be quite complex as there are many potential configurations of both the density driven and advective flow components. Mixed convection, where both advection driven by hydraulic gradients and convection driven by density gradients control resultant fluid motion and solute distributions have been studied extensively in fluid mechanics (see for example, Nield and Bejan, 1999) but also in groundwater hydrology (Schincariol and Schwartz, 1990; Oostrom et al., 1992a,b). In the laboratory, downward displacement experiments of fresh water by denser miscible solutions have been carried out (so-called unstable displacements) and can give rise to instabilities (Jiao, 2001; Jiao and Ho ¨tzl, 2004; Wood et al., 2004, amongst many others). Likewise, so-called stable displacement experiments where fresh water is displaced by brine from below have been performed by a number of authors (e.g., Kempers and Haas, 1994; Watson and Barry, 2001; Watson et al., 2002 and Jiao and Ho ¨tzl, 2004). Watson et al. (2002) found that Darcy’s law is valid for high concentration displacements. Despite a large body of both theoretical and laboratory based work relating to mixed convection phenomena, the systematic application of mixed convection theoretical developments to the study of saline groundwater discharge into freshwater streams and lakes is both largely unexplored and poorly understood. The very limited number of saline groundwater discharge studies that include density effects in a stream environment (e.g., Jolly et al., 1998) do so in a site-specific way. It therefore remains a challenge to draw generalised conclusions

392 about the hydrogeological controls and ensuing hydrodynamics processes that govern saline groundwater discharge into freshwater streams. In this paper, we study surface water/groundwater interaction processes that occur in a mixed convective system characterised by a stably-stratified fresh surface water–saline groundwater discharge zone. Within a theoretical framework, we use a groundwater flow and solute transport model to determine how saline groundwater discharge to a surface water body varies as a function of the mixed convective ratio M.

Conceptual and theoretical framework Conceptual model A conceptual model, which is representative of the situation at any gaining freshwater stream or lake underlain by saline groundwater (Fig. 1) can be established. The hydraulic head of the groundwater (h1, equivalent freshwater head, EFH, corrected according to measured piezometric head and salinity) is higher than the surface water head (h0) at any gaining surface water body. Hence, the flow direction resulting from the head gradient (EFH gradient) only would be upwards towards the surface water feature if it is a discharge zone. The salinity of the groundwater is greater than the surface water salinity. Therefore, the water with the highest concentration of total dissolved solids (c1) and likewise the largest density underlies the less saline and less dense water with a total dissolved solids (TDS) concentration c0.

G. Massmann et al. We begin by considering the simplest possible geometrical configuration, namely a purely one-dimensional column that represents a vertical flow line directly beneath the discharge zone. This leads to a numerical model implementation in the form of a vertical column with a height of z = 10 m and a width of x = 1 m, representing the final groundwater discharge point (Fig. 1). Since groundwater flow paths directly beneath the surface water body are sub-vertical to vertical, horizontal flow is disregarded in this initial analysis. The top and bottom of the column represent the boundaries of the surface water body (top) and the regional groundwater (bottom). The surface water hydraulic head and salinity are often well known and less variable in space and time than the groundwater head and salinity, which may be highly variable. The groundwater head (h1) and TDS content (c1) were therefore varied systematically in a matrix of simulations (Fig. 2), whilst the equivalent freshwater head of the surface water body (h0) and TDS concentration (c0) were held constant. The aim here was to systematically study the effects of the hydraulic gradient and the concentration (i.e., density) gradient on the groundwater discharge behaviour. This model assumes homogeneous and isotropic porous media conditions, isothermal and viscosity-invariant conditions. The surface water system is not explicitly modelled but is rather represented by a constant head boundary condition.

FEFLOW numerical model and governing equations Numerical modelling was performed using FEFLOW (Diersch, 2002a), a numerical finite element simulation system for two-dimensional (2D), horizontal, vertical and axisymmet-

Figure 1 Development of simplified generalised modelling approach from the schematised ‘‘real’’ situation at any gaining stream or lake underlain by saline groundwater. Numerical model implementation is shown to the left; the groundwater equivalent freshwater head (h1) and concentration (c1) are larger than those of the river (h0, c0), the shape of the depth profiles is assumed unknown and non-linear with both represented by some function of space and time.

On variable density surface water–groundwater interaction

393

Figure 2 Matrix of simulations with mixed convection ratio M given above each cross, each representing a simulation performed under the corresponding hydraulic and density gradient. The threshold of M = 1 marks the boundary between advection and density driven flow dominated simulations.

ric, or three-dimensional (3D) flow and transport processes in porous media. FEFLOW simulates transient or steadystate, fluid density-coupled or linear, flow and mass, flow and heat or completely coupled thermohaline transport. It has successfully been tested for a number of benchmark examples in variable density flow, such as the Elder problem, the Henry problem, the salt dome problem and the salt pool problem (Diersch, 2002b). The fundamental mathematical basis of FEFLOW (Diersch, 2002b) is formed by the physical principles of: mass conservation of fluids: o o ðea qa v ai Þ ¼ ea qa Q aq ðea qa Þ þ ot ox i

ð1Þ

mass conservation of solids: o o o a ðea Cak Þ þ ðea v ai Cak Þ þ ðj Þ ¼ ea Rk ot ox i ox i ik momentum conservation of fluid and solid continua:   kaij opa a a vi þ  q gj ¼ 0 ea la oxj

ð2Þ

ð3Þ

and energy conservation: o o o a ðea qa v ai E a Þ þ ðj Þ ¼ ea qa Q aT ðea qa E a Þ þ ot ox i ox i iT

ð4Þ

FEFLOW uses the following density dependent formulation of the Darcy law (Diersch, 2002b):   oh qf  qf0 qfi ¼ K ij fl ð5Þ þ e j oxj qf0 In practise, the density differences are accounted for by the fluid density difference ratio a entered in the flow material menu, which is defined as (Diersch, 2002b):

a¼

ðqf ðCs Þ  qf0 Þ qf0

where: t x ea

ð6Þ

= time [T] = Eulerian spatial coordinate vector [–] = volume P fraction of a-phase, where 0 6 ea 6 1 and a ea ¼ 1 [1] qa = density of a-phase [ML3] a vi = velocity vector of a-phase [LT1] Qaq and QaT = mass and heat supply of a-phase [ML3T1] and [L2T3] a Ck = concentration of component k in phase a [ML3] a jik = Fickian mass flux vector of a-phase [ML2T1] jaiT = Fourierian heat flux vector of a-phase [MT3] Rk = macroscopic reaction rate for homogeneous and heterogeneous chemical reactions of component k: Rk ¼ rak þ Rhet [ML3T1] k k aij = permeability tensor of a-phase [L2] la = dynamic viscosity of a-phase [ML1T1] a p = thermodynamic pressure of a-phase [ML1T2] gi = gravity vector [LT2] a E = internal (thermal) energy of a-phase [L2T2] f qi = Darcy velocity of fluid phase [LT1] Kij = tensor of hydraulic conductivity of fluid phase [LT1] fl = viscosity relation function [–] qf = density of fluid phase [ML3] f q0 = reference density of fluid phase related to the reference head h0, the reference concentration C0 and the reference temperature T0 [ML3] qf(Cs) = density of fluid phase with maximum concentration CS [ML3] ej = gravitational unit vector [1] h = hydraulic head referenced to q0 [L]

394 For further information on the FEFLOW numerical code the reader is referred to Diersch (2002b).

Mixed convection ratio In our study, both free convection (purely driven by density differences) and forced convection or advection (purely driven by hydraulic head gradients) are important. In such a mixed-convective system, the mixed convection ratio can serve as a first simple measure of the dominance of either of the two driving forces. The mixed convection ratio M is a dimensionless parameter and is the ratio of the density driven flow speed to the advective flow speed and for a homogeneous and isotropic system can be expressed as:  f      q  qf0 ohf Dq Dhf M¼ ¼ ð7Þ q0 ox j Dl qf0 where Dq is the density difference between the heavy and lighter fluids [ML3], q0 is the density of the lighter fluid [ML3], Dhf is the hydraulic EFH difference [L], Dl is the length over which hydraulic head (EFH) difference is measured [L]. Note here that this is clearly found from the variable density version of Darcy’s Law (Eq. (5)) as the ratio of the density driving force (right-hand side of bracketed equation) to the advective driving force (left-hand side of bracketed equation). Importantly, the denominator is therefore clearly the gradient of EFH, and not any other ‘‘point’’ or ‘‘measured’’ head. This point is not entirely clear in earlier literature on mixed convection ratios (e.g., Nield and Bejan, 1999; Bear, 1972) and it is ambiguous as to how to correctly define the mixed convection ratio in a mixed convective system, especially in defining the appropriate hydraulic head (measured or corrected to equivalent freshwater head) to be used. Care must therefore be exercised in calculating this ratio, paying attention to ensure EFH are used in the analysis. The mixed convection ratio can then be used to describe the relative important of density driven and advectively driven flow. Critically, where M < 1, advection dominates and where M > 1 density driven flow dominates.

Numerical modelling Model discretisation, boundary conditions and parameters The model domain was discretised to form a completely uniform mesh containing 2703 mesh nodes and 2528 square mesh elements, 16 in the x direction (Dx = 0.0625 m) and 158 in z-direction (Dz = 0.0625 m). Whilst the model is quasi 1D, the simulations were performed as truly two-dimensional problems, using the FEFLOW vertical (cross-sectional) problem projection mode. All simulations were performed as transient flow, transient transport simulations and were run in transient mode until steady-state conditions were reached. For the temporal approximation, the forward Adams–Bashforth/backwards trapezoid rule predictor corrector scheme with adaptive time stepping (error tolerance 1005, initial time step length 0.001 days, maximum time-

G. Massmann et al. step length 0.1 days) was used (Diersch, 2002b) which helped to eliminate numerical oscillations that may have otherwise occurred. A sensitivity analysis to spatial and temporal discretisation revealed that the numerical solution had converged for the values chosen and that obtained solutions were independent of the choice of spatial and temporal discretisation. The top and bottom boundaries were implemented as Dirichlet or 1st kind boundaries for flow and transport with a constant head (h0, h1) and concentration (c0, c1). The left and right-hand side of the model were implemented as impermeable (no flow) boundaries. The initial concentration was set to 0 mg L1 throughout the entire model domain. The resultant pressure distribution of a steady-state simulation performed with the chosen boundary conditions and including density effects was used as the initial pressure condition for transient transport simulations. The model domain was considered to be homogeneous, with uniform hydraulic conductivity, porosity and dispersivity. The densities were taken from Weast et al. (1984–1985). A list of complete model parameters is given in Table 1. In order to perform model runs with a number of hydraulic gradients, i = (h1  h0)/z, as well as concentration gradients, Dc = c1  c0, the head (EFH) and the TDS concentration of the groundwater were varied systematically, while all other parameters were kept constant. Molecular diffusion was ignored, since its effects are small compared with those of dispersion in typical field scale settings.

Model simulations A matrix of simulations was developed to cover a realistic range of both hydraulic and density gradients and hence mixed convection ratios. Hydraulic gradients of 0.1, 0.01, 0.001 and 0.0001 were adjusted by varying the constant head boundary at the bottom (h1). The groundwater concentration was set to 297,200, 100,200, 29,500, 10,100, 3000 and 1000 mg L1, resulting in densities of 1188.7, 1066.2, 1018.9, 1005.3, 1000.4 and 998.9 kg m1 (Weast et al., 1984–1985) and density ratios (a) of 0.19081, 0.06809, 0.02007, 0.00708, 000217 and 0.00067, respectively. These density ratios in combination with the hydraulic gradients resulted in 24 simulations, each marked with a cross in Fig. 2. Above each cross, the mixed convection ratio of the corresponding combination of i and a is given. For the simulations above the dotted line (M < 1), forced convection should be dominant, while for those below the line (M > 1), free convection should dominate. In this way, the full implications of density dependent analysis are explored to ultimately allow the distinction between advection and density dominated cases. The entire set of simulations was performed for dispersivities (b) of 0.05, 0.5 and 1 each, covering a realistic range of b under the current settings, without violating the grid Peclet Number. The dispersivities were varied in order to evaluate the effect of dispersivity on the shape of depth profiles and breakthrough curves. All simulations were finally repeated with a density difference ratio (a) of 0 while maintaining the EFH values, to fully compare density dependent results with constant density (density invariant) cases.

On variable density surface water–groundwater interaction Table 1

395

Model input parameter

Parameter

Symbol

Value

Unit

Model: Height/aquifer depth Width Cell width

z x Dz, Dx

10 1 0.0625

m m m

Solid properties Hydraulic conductivitya Effective porositya Long. dispersivityb Transv. dispersivityb Molecular diffusivity

Kf e bL bV Df

1a E  03 0.2 0.05, 0.5, 1 0.05, 0.5, 1 0

m s1 – m m m2 s1

Fluid properties Surface water headc Groundwater headd Hydraulic gradient Surface water TDS, reference conc. Groundwater TDS Concentration difference Surface water density Groundwater density Density ratio

h0 h1 i c0 c1 Dc q0 q1 a

9.9 Varied, 9.901–10.9 1E  04 to 1E  01 0 varied, 1000–29,7000 1000–29,7000 998.23 998.9–1188.7 6.71E  04 to 1.91E  01

m AHD m AHD – mg L1 mg L1 mg L1 kg m3 kg m3 –

a b c d

Based upon realistic values for a coarse sand (Freeze and Cherry, 1979). Based upon realistic values for a coarse sand and a domain length of 10 m (Gelhar, 1984). Based upon water level of River Murray which was initially the focus of the study. Equivalent freshwater head.

Model testing In order to check model performance, simulations were initially performed without the upper Dirichlet boundary and without consideration of density effects (a = 0). The results were compared to the analytical solution by Ogata (1970). Excellent agreement between breakthrough curves from the analytical solution and the numerical modelling were obtained (not shown), as expected from a similar benchmark problem given in Diersch (2002b).

Modelling analysis In order to facilitate a quantitative analysis, the results of numerical simulations were compared in terms of the: – Occurrence of mass breakthrough throughout the upper model boundary. – Shape of the depth profiles of mass (TDS) at steady-state conditions: In order to find a quantitative, comparable measure of the depth profile, the solute present (SP) and the center of mass (c.o.m.) were calculated from the steady-state depth profiles as suggested by Prasad and Simmons (2003). The solute present is proportional to the amount of solute present in the aquifer at any time. It is the area under the concentration versus depth graph, equivalent to the amount of solute present in the aquifer per unit cross-sectional area (per unit coefficient of density change). It is the cumulative sum of the products computed by multiplying the average concentration with the length of the corresponding depth interval.

From the solute present, the c.o.m. can be calculated. The c.o.m. of a depth profile is the vertical center of the total solute present, i.e., the depth (below river or lake base) above and below which 50% of the solute present is equally distributed. – Time dependent breakthrough curves at an observation point at the centre of the model domain (x = 0.5 m; z = 5 m): The breakthrough curves are described by the times t0.1, t0.5 and t0.9 where 10%, 50% and 90% of the input concentration have broken through, respectively. The time t0.5 is the time where water travelling at the average pore water velocity would arrive (Freeze and Cherry, 1979). The 10% and 90% breakthrough are considered in order to account for the slope of the breakthrough curves which is an indicator of breakthrough rate. – Mass discharge into the surface water body (i.e., across the upper model boundary) at steady-state conditions taken from the FEFLOW budget analyser tool. All depth profiles and breakthrough curves were normalised for ease of comparison by dividing the results obtained by the variable input TDS concentration of the groundwater c1 and are therefore given as c/c1.

Results Mass breakthrough The first observation that can be made from each simulation is whether or not a breakthrough of mass occurs through the

396

G. Massmann et al.

Figure 3 Breakthrough of mass throughout model domain observed at steady-state conditions in the simulations plotted on top of simulation raster (Fig. 2).

upper model boundary (Fig. 3). This is determined by the dominant direction of flow. These observations fundamentally dictate the spatiotemporal patterns of salt fluxes and discharge patterns and hence bear direct influence over the results obtained for other quantitative indicators such as c.o.m. or solute present. In simulations with mixed convection ratios M < 1, the flow direction is vertically upwards towards the surface water and mass breakthrough through the upper model boundary occurs after some time before steady-state conditions are finally reached. In contrast, when M > 1 and flow is expected to be density dominated, no breakthrough of mass occurs. In the cases with substantially larger M (M larger than 50, depending slightly on the dispersivity), the flow direction actually reverses and is directed downwards towards the groundwater. Hence, in these cases, the saline groundwater does not discharge into the surface water, even though the advective hydraulic gradients alone would suggest discharge should otherwise occur. In cases which plot below and close to the critical threshold of M = 1 in Fig. 2, the flow is essentially stagnant and numerical problems occur as soon as steady-state conditions are reached. This is exacerbated in simulations performed with lower dispersivities.

Mass depth profiles Density invariant profiles In simulations where density effects were neglected (a = 0) whilst maintaining the applied EFH gradient, the flow direction is upwards towards the upper boundary following the hydraulic gradient. The mass shape of the normalised depth profiles at steady-state is then a function of the dispersivity b only. Hence, all steady-state density-invariant simulations performed with identical dispersivity values lie on top of each other, even though they were performed at different hydraulic gradients i. From a simple analysis of the advection–dispersion equation under steady-state conditions, it

Figure 4 Density invariant depth profiles for simulations performed with a range of dispersivities (a), showing the effect of dispersivity on the shape of the depth profile of mass.

follows that the resultant concentration solution in space is velocity independent but dispersivity dependent. The analysis is simple enough not to warrant derivation here. Depth profiles at steady-state for simulations performed with a range of dispersivities are shown in Fig. 4, in order to show their effect on the slope and shape of the profiles and clearly demonstrate that the smaller b, the steeper the concentration gradient towards the river boundary. Density variant profiles Fig. 5 shows the normalised depth profiles of mass in steadystate for the simulations performed with a density difference ratio a > 0 and a dispersivity b of 0.05 m. Similar to simulations performed with a = 0, the depth profiles of those simulations with M < 1 (advection dominated) plot on top of each other towards the upper model boundary, exactly where the density invariant cases plot. Hence, their c.o.m. is located closest to the upper model/surface water boundary (i.e., c.o.m. shallowest when given in m below surface water base). But with increasing M, the saltwater/fresh-water interface migrates vertically downwards towards the lower model boundary. Finally, in simulations with a very large M (M > 68.1), the profiles (and likewise the centre of mass) plot within the lower twentieth part of the domain (below 9.5 m depth below river). The semi-logarithmic plot of c.o.m. versus mixed convection ratio (M) illustrates that these effects are essentially identical for simulations with similar M, but

On variable density surface water–groundwater interaction

Figure 5

397

Steady-state depth profiles for all simulations with a density difference ratio of a > 0 and a dispersivity of b = 0.05 m.

different dispersivities (b). The only apparent difference is that the salt-water/freshwater interface is smoother when b is larger (not shown) and rather sharp for smaller b (Fig. 6). Consequently, the c.o.m. is smaller (closer to the upper boundary, i.e., the surface water body) when a is smaller and M < 1 and deeper (i.e., further towards the lower boundary, i.e., the groundwater) when M > 1. Fig. 6 shows that the greatest changes in the shape of the depth profile and the location of the c.o.m. occur when M is just

above 1, while the c.o.m. versus M graph approaches a plateau in the lower tenth of the domain for larger values of M.

Mass breakthrough curves A breakthrough curve of local mass can only be determined in the advection dominated cases (M < 1) because in the density dominated simulations (M > 1), breakthrough of mass does not occur (Fig. 2). Density invariant and variant

Figure 6 Centre of mass versus mixed convection ratio on a semi-logarithmic plot. The threshold of M = 1 is indicated by the dotted line.

398 cases were compared. In order to maintain clarity, only three examples for breakthrough curves with and without density are shown in Fig. 7. They represent the cases for the three hydraulic gradients were breakthrough occurs (i = 0.1, 0.01 and 0.001) for simulations with dispersivities of 1 m (M = 0.67, 0.71 and 0.68; see Fig. 2). For those cases where the hydraulic gradient is similar but the density differences and hence M are smaller but non-zero, the breakthrough curves lie between the two lines shown in Fig. 7. Hence, as M approaches (but does not exceed 1), the breakthrough occurs later than in the density invariant cases. The density variant breakthrough curves also have a lower inclination (i.e., they are less steep). It can generally be said that the effect of density is to suppress or retard the discharge and that the breakthrough at the observation point occurs later when density effects are considered. Fig. 8 demonstrates graphically the idea that at a given hydraulic gradient, the breakthrough occurs later as the concentration gradient (i.e., density difference ratio and therefore mixed convection ratio) gets larger and M approaches 1. Breakthrough time is plotted against the mixed

G. Massmann et al. convection ratio M. Instead of using t0.5, both t0.1 and t0.9 are given in order to account for the curvature of the breakthrough curves (solid and hollow symbols, respectively). The figure is grouped into plots of similar hydraulic gradient (4 graphs for each row in Fig. 2). Grey symbols and dotted lines represent density invariant simulations and black symbols and solid lines denote density variant simulations. For simplicity, only results from simulations with a dispersivity (b) of 1 are shown. In the density invariant cases, the density difference does naturally not matter and the breakthrough curves lie on top of each other as discussed previously. Hence, the slope of the dotted lines is always 0. Compared to these density invariant cases, the time difference (t0.9  t0.1) is larger when the density is considered, i.e., the black symbols (and solid lines) are further apart from each other than the grey symbols. This reflects the smaller slopes in Fig. 7. It also becomes clear by the increasing slopes of the solid lines that at a given hydraulic gradient, the breakthrough occurs later with increasing mixed convection ratios (i.e., with increasing density differences between surface water and

Figure 7 Breakthrough curves at x = 0.5 and z = 5 for three example cases below but close to threshold of M = 1 (see Fig. 6). The time-scale is logarithmic. Simulations were performed for both density variant and density invariant cases.

Figure 8 Time of breakthrough of 10 (solid) and 90 (hollow)% of the total solids with (black symbols, solid lines) and without (grey symbols, dotted lines) accounting for density. The graphs are grouped into graphs of similar hydraulic gradients i. Diamonds for i = 0.001, triangles for i = 0.01 and circles for i = 0.1.

On variable density surface water–groundwater interaction

399

groundwater). Increasing groundwater TDS concentrations would therefore slow down the breakthrough of salt until the threshold of M = 1 is passed and the breakthrough is completely suppressed. Column studies of stable miscible displacement experiments conducted by Jiao and Ho ¨tzl (2004) showed that effective dispersion coefficients decrease as the density difference increases. Similar observations were made in column studies by Watson and Barry (2001). While the depth profiles Fig. 5 also show that dispersion is smaller for the cases with larger stable density differences between surface and groundwater (M large), the time breakthrough curves of mass as discussed above appear to show the opposite behaviour and are also rather asymmetric in comparison to the density invariant cases. The reason for this inconsistency is somewhat unclear. Density differences used by Jiao and Ho ¨tzl (2004) are much smaller in comparison to those used in this paper and this may, at least in part, account for some of these differences.

water TDS concentration. However, very large TDS concentrations suddenly have the reverse effect, since the induced density dependent flow causes the water to sink (or stay where it is), which inhibits vertical discharge into the surface water. In comparison to the effect of hydraulic and density gradients, the effect of variable dispersivity on mass flux calculations is very minor. These results show that a maximum in salt discharge does indeed occur at about M = 1 but that as M exceeds one, density effects result in a reduced net upwards vertical flow (or ultimately no flow) and hence there is an immediate reduction in salt discharge. Thus, it is immediately apparent that increasing groundwater concentrations further does not increase salt discharge but rather suppresses it. This is a critical observation that would not be seen in a density invariant case where salt discharge would continue to increase linearly with increasing groundwater concentration. It is interesting to characterise the different cases that would control salt discharge behaviour. By noting that the salt flux Sf through the upper model boundary (i.e., into the surface water) is simply the fluid flux (Eq. (8)) multiplied with the groundwater TDS concentration:   oh qf  qf0 Sf ¼ qfi  c1 ¼ K ij fl  c1 þ e ð8Þ j ox j qf0

Salt discharge into the fresh surface water From a practical point of view, the amount of salt discharged into the river or lake is often an important concern. Fig. 9 shows the salt flux through the upper boundary (i.e., out of the model domain) as a function of the mixed convection ratio M separated into plots for the different dispersivities (b) and hydraulic gradients (i). This would, in the broader sense, represent the discharge of salt into the surface water. The discharge is largest for the largest hydraulic gradient i (Fig. 9, hollow symbols). Within a given i and b combination, the salt flux increases with increasing mixed convection ratio until M passes the threshold of 1. This increase is due to the fact that at a given i, the increase of M is caused by an increase of Dc (or Dq/q0) since these are directly proportional (Eq. (7), leading to a greater mass flux. When M becomes greater than 1, breakthrough of mass through the upper model boundary does not occur and the mass flux out of the model area through the upper model boundary becomes zero. Hence, there would be no discharge of salt into the surface water. This means that at a given i, the salt discharge increases with increasing ground-

There are four possible cases which can be formulated to summarise salt discharge behaviour: Dh (1) If a ¼ Dx and M = 1, then Sf = 0   j oh (2) If a = 0 and M = 0, then Sf ¼ K ij ox  c1 , which is the j traditional density invariant formulation of Darcy’s law Dh (3) If 0 < a < Dx and 0 < M < 1, then Sf ¼ qfi  C1 j Dh (4) If a > Dx and M > 1, then the term in brackets in Eq. j (8) becomes negative and there is no salt discharge.

Discussion Assumptions, limitations, usefulness In this paper, we started with the simplest 1D approach to analyse surface water groundwater interactions since these

Figure 9 Salt flux (g d1) through the upper model boundary (area of 1 m2) as a function of the mixed convection ratio (M) for different dispersivities (b) and hydraulic gradients (i). (b) Is a zoom of (a).

400 must clearly be understood before proceeding to more complex 2D and 3D analyses. Our current analysis assumes that the system being considered is strongly dominated by vertical flow phenomena, as might be the case under a large extensive lake system in a deep sediment layer where most flow is vertical to subvertical. Clearly, as the aspect ratio of the surface water feature becomes smaller, or possibly less than one (where vertical dimensions exceed horizontal dimensions) horizontal flow components become more significant relative to the vertical flow components. Thus, our study is representative of the ‘‘flow lines’’ that exist beneath large laterally extensive discharge zones where vertical flow is dominant and where horizontal flow components, say at the edges of the lake, are small in comparison. The concept of a mixed convection ratio is demonstrated theoretically in this paper as a valid means of determining where density affects are dominant (where M > 1) and where further more detailed analyses incorporating fully-coupled variable density flow analyses may be warranted. We suggest that the mixed convection ratio M is therefore a most useful indicator for determining the nature of the surface water/ groundwater interaction and whether it is advective controlled (by hydraulic head gradients), or density controlled (by density driven gradients). This paper introduces the mixed convection ratio, a dimensionless number typically used in fluid mechanics and in unstable variable density flow analyses, to the study of stable fresh surface water/saline groundwater mixing phenomena – an area that has not previously been explored. This study is the first we are aware of that explicitly demonstrates how M can be used to predict surface water groundwater interactions in a stably stratified system. Preliminary analyses in 2D suggest complex variations on these conclusions may result where large horizontal flow components occur in the mixed convective regime. It appears from a number of simple models that the nature of the interaction is highly sensitive to the aspect ratio of the domain, the aspect ratio of the discharging boundary (e.g., flat boundary, steep boundary, ‘‘L’’ shaped discharge boundary and ratio of horizontal to vertical discharge extent at that boundary). Further work will be required to fully elucidate all controlling features (geometric, hydrogeologic) of the interaction in a fully 2D (and 3D) system. However, it is expected that the multidimensional analysis will contain two major components: the vertical flow components and the horizontal flow components. The conclusions drawn about density affects in this study are likely to apply in multiple dimensions since in that case, the system may be considered as the superposition/summation of many flow lines each with horizontal and vertical components. The vertical flow components are likely to be affected in the same way that is demonstrated in this paper, namely, that as the density difference increases, the net discharge rate is reduced. However, horizontal flow components are likely to be relatively unchanged. Thus, one likely possibility although speculative, is that density differences only substantially alter vertical flow components whilst maintaining horizontal flow. Therefore, complete suppression or reversal of flow is not likely to occur in 2D and 3D systems. However, this conclusion is not yet supported by rigorous numerical modelling and clearly warrants further detailed investigation in future studies. However, even in 1D, the potential for use of

G. Massmann et al. M and the likely impact of density effects on discharge rates and mechanisms has been demonstrated. These must now be extended into subsequent multidimensional analyses coupled with a rigorous analysis of how dimensionality affects these processes and the generalised conclusions of these phenomena. At this stage, it is too early to draw generalised conclusions in light of the potentially complicating factors mentioned above.

A note on the variable-density version of Darcy’s Law and equivalent freshwater heads It is useful to provide a brief remark here on the density dependent formulation of Darcy’s Law (Eq. (5)). In physical terms, this equation shows that there are two key components: the usual head gradient (left-hand side of bracketed expression) and an additional fractional excess density term (or density gradient as seen in the right-hand side of bracketed expression). In the vertical, this density driving force acts in the direction of the gravitational unit vector as seen in the right hand part of the bracketed expression and is sufficient to drive flow even in the absence of external hydraulic head gradients. It must be noted that the hydraulic head that is used in this expression is the equivalent freshwater head (EFH). In calculating flow rates and directions in a variable density system, it is therefore not sufficient to simply correct measured heads to EFH and to use them within the familiar density-invariant form of Darcy’s Law that does not contain the additional density gradient. Rather, EFH must be employed within the fully density dependent version of Darcy’s Law described in Eq. (5). This is particularly important in variable-density flow analyses that contain any element of vertical flow or vertical stratification with respect to fluid density and is a potential source of error in variable density flow calculations. In our case studied here, vertical flow phenomena are the key, and use of EFH without care, could be particularly problematic. Indeed, Lusczynski (1961) recognized that equivalent freshwater heads cannot be used to determine the vertical hydraulic gradient in an aquifer with water of non-uniform density. The term environmental water head (EWH) was established for this purpose (Lusczynski, 1961). However, it can be demonstrated that this environmental water head is equivalent to using Equivalent Freshwater Head where the appropriate density dependent version of Darcy’s Law is used (as was actually shown in Appendix 3 of the original Lusczynski study). To that end, the variable density version of Darcy’s Law is therefore used in many current variable density groundwater flow models such as that employed in this study, FEFLOW. However, a common error that occurs in practice outside the use of these standard simulators is for equivalent freshwater head analyses to be used in the standard non-density dependent version of Darcy’s Law. This leads to errors in predicting flow rates and directions where vertical flow effects are of interest in vertically stratified layers (with respect to density) of fluid. Thus, when vertical flow phenomena are of interest in variable density flow problems, the equivalent fresh water head must be used with particular caution, and most importantly, with the correct version of Darcy’s Law that accounts for the density gradient that acts along the vertical. This point can hardly be overstated.

On variable density surface water–groundwater interaction

Summary and conclusion This study has considered the role of variable density flow on the discharge behaviour of saline groundwater into fresh surface water bodies. A numerical 1D approach was used to study how surface water/groundwater interaction processes and salt discharge behaviours were affected by the inclusion of density driven flow phenomena. The study demonstrates in a theoretical manner that: 1. The presence of density gradients can significantly modify the behaviour of the interaction between the saline groundwater and fresher surface water body; 2. It is possible for density driven flow to suppress discharge rates, where density gradients may act in the opposite direction to advection in the case where the surface water body is a discharging feature; 3. The mixed convection ratio M is an effective tool for establishing dominant hydrodynamic controls in a mixed convective flow system where density and hydraulic gradients simultaneously control flow and transport and, in this case, act in opposite directions. It may therefore be used as a simple criterion to predict the extent of the influence of variable density on salt discharge into a river or lake and can act as a useful guide as to whether more detailed analyses should incorporate density dependent processes. For M > 1, density driven flow is dominant. For M < 1, advection is dominant, a fact consistent with prior theoretical work conduced in fluid mechanics for many decades. 4. Mass breakthrough only occurs when M < 1, where advection is the driving force and the flow direction is towards the stream. When M > 1 and the flow is dominated by the density difference between the two fluids, no breakthrough of mass occurs, since the resultant flow is either so small so as to be effectively stagnant or reverses and points away from the surface water body. 5. Decreasing M moves the saltwater-freshwater interface up towards the freshwater model boundary until M becomes smaller than 1. For all advection-dominated cases (M < 1), the steady-state depth profiles are identical and their shape is dictated by the dispersivity. Larger dispersivities increase the slope and width of the saltwater-freshwater transitional interface. 6. The transient breakthrough of salt is delayed when density differences are considered. As M increases (but does not exceed 1) at a given hydraulic gradient, the breakthrough of salts is delayed in comparison to the nondensity case and the slope of the breakthrough curve is less steep. 7. A critical salinity may exist for maximum salt discharge, an effect not seen in a density-invariant analysis which assumes salt flux increases linearly with increasing salinity of groundwater and advective flow rate; the absolute amount of salt discharged is primarily a function of Dc and i (and therefore M) as well as the hydraulic conductivity (kf) of the riverbed. In simulations with constant hydraulic gradients, hydraulic conductivities and dispersivities of the sediment, the salt flux increases with increasing Dc (and hence M) until M exceeds the critical transitional threshold of 1. In simulations with constant

401 concentration differences, hydraulic conductivities and dispersivities of the sediment, the salt flux decreases with decreasing i (and hence increasing M), again until M exceeds the critical threshold of 1, where the flux then becomes zero. Thus, this study clearly demonstrates that increasing groundwater concentrations beyond some critical defined by the condition M = 1 does not increase salt discharge but rather suppresses it. This is a critical observation that would not be seen in a density invariant case where salt discharge would continue to increase linearly with increasing groundwater concentration. This ‘‘critical concentration’’ effect on salt load has, to the best of our knowledge, not been reported in previous literature. These observations are, however, only valid for vertical flow phenomena. 8. Equivalent freshwater heads alone cannot be used without a modified version of Darcy’s Law that fully accounts for the density gradient explicitly. The results of this study may find wider application in other surface water/groundwater interaction studies where variable density flow may be important in controlling the hydrodynamics of the groundwater flow and solute transport processes. Thus, whilst this study began with the motivation of examining groundwater/stream interactions, it is possible that similar analyses may be useful in applications including, but not limited to, lake–groundwater interactions, subterranean groundwater discharge, artificial recharge processes and other applications where both advection and density driven flow simultaneously control the mixed convective regime. Since the chosen approach is purely theoretical, it would be useful to conduct subsequent laboratory experiments under boundary conditions similar to those in the simulations. In addition, conclusion will need future confirmation by rigorous 2D and 3D numerical modelling as discussed above.

Acknowledgements We thank colleagues from the Department of Water, Land and Biodiversity Conservation of the Australian Government and from the Flinders University of South Australia for helpful discussions, in particular Dr. Neville Robinson from Flinders University. Many thanks to the operators of the Feflow hotline at Wasy GmbH in Berlin for their continuously friendly and patient support. Finally, we thank two anonymous reviewers for their constructive and helpful comments.

References Allison, G.B., Cook, P.G., Barnett, S.R., Walker, G.R., Jolly, I.D., Hughes, M.W., 1990. Land clearance and river salinisation in the western Murray Basin, Australia. J. Hydrol. 119, 1–20. Bear, J., 1972. In: Dynamics of Fluids in Porous Media. Elsevier Science, New York, p. 764. Bourg, C.M., Bertin, C., 1993. Biochemical processes during the infiltration of river water into an alluvial aquifer. Environ. Sci. Technol. 27, 661–666. Cheng, J.M., Chen, C.X., 2001. Three-dimensional modeling of density dependent salt water intrusion in multilayered coastal

402 aquifers in the Jahe River Basin, Shandong Province, China. Ground Water 39 (1), 137–143. Diersch, H.-J.G., 2002a. Interactive, Graphics Based Finite Element Simulation System FEFLOW for Modelling Groundwater Flow, Contaminant Mass and Heat Transport Processes. WASY Ltd., Berlin. Diersch, H.-J.G., 2002b. Reference manual for FEFLOW – finite element subsurface flow and transport simulation systemUser’s Manual/Reference Manual/White Paper Release 5.0. Wasy Ltd., berlin. Diersch, H.-J.G., Kolditz, O., 2002. Variable-density flow and transport in porous media: approaches and challenges (Review paper). Adv. Wat. Res. 25 (2002), 8–12, 899–944. Diersch, H.-J.G., Prochnow, D., Thiele, M., 1984. Finite-element analysis of dispersion-affected saltwater upconing below a pumping well. Appl. Math. Mod. 8, 305–312. Dousson, C., Poitevin, G., Ledoux, E., Detay, M., 1997. River bank filtration: modelling of the changes in water chemistry with emphasis on nitrogen species.. J. Contam. Hydrol. 25, 129–156. Freeze, R.A., Cherry, J.A., 1979. Groundwater. Prentice-Hall, Englewood Cliffs, New Jersey. Holzbecher, E., 1998. Modeling Density-driven Flow in Porous Media. Springer, Berlin. Huyakorn, P.S., Anderson, P.F., Mercer, J.W., White, J.H.O., 1987. Saltwater intrusion in aquifers: development and testing of a three-dimensional finite element model. Water Resour. Res. 23, 293–312. Jacobs, L.A., von Gunten, H.R., Keil, R., Kuslys, M., 1988. Geochemical changes along a river-groundwater infiltration flow path: Glattfelden, Switzerland. Geochim. Cosmochim. Acta 52, 2693–2706. Jiao, C., 2001. Miscible displacements in porous media with variation of fluid density and viscosity. Ph.D. Thesis, University of Karlsruhe, Karlsruhe. Jiao, C.-Y., Ho ¨tzl, H., 2004. An experimental study of miscible displacements in porous media with variation of fluid density and viscosity. Transport Porous Med. 54, 125–144. Jolly, I.D., Narayan, K.A., Armstrong, D., Walker, G.R., 1998. The impact of flooding on modelling salt transport processes to streams. Environ. Model. Softw. 13, 87–104. Kempers, L.J.T.M., Haas, H., 1994. The dispersion zone between fluids with different density and viscosity in a heterogeneous porous medium. J. Fluid Mech. 267, 299–324. Kolditz, O., Ratke, R., Diersch, H.-J.G., Zielke, W., 1998. Coupled groundwater flow and transport: 1 Verification of variable-density flow and transport models. Adv. Wat. Res. 211, 27–46. Liu, H.H., Dane, J.H., 1996. A criterion for gravitational instability in miscible dense plumes. J. Cont. Hydrol. 23 (3), 233–243. Lusczynski, NJ., 1961. Head and flow of ground water of variable density. J. Geophys. Res. 66 (12), 4247–4256. Nield, D.A., Bejan, A., 1999. In: Convection in Porous Media. Springer, New York, p. 546. Ogata, A., 1970. Theory of dispersion in granular medium. US Geol. Surv. Prof. Paper 411 (I). Oldenburg, C.M., Pruess, K., 1995. Dispersive transport dynamics in a strongly coupled groundwater-brine flow system. Wat. Resour. Res. 31, 289–302. Oostrom, M., Dane, J.H., Gu ¨ven, O., Hayworth, J.S., 1992a. Experimental investigation of dense solute plumes in an unconfined aquifer model. Wat. Resour. Res. 28 (9), 2315–2326.

G. Massmann et al. Oostrom, M., Hayworth, J.S., Dane, J.H., Gu ¨ven, O., 1992b. Behaviour of dense aqueous leachate plumes in homogeneous porous media. Wat. Resour. Res. 28 (8), 2123–2134. Prasad, A., Simmons, C.T., 2003. Unstable density-driven flow in heterogeneous porous media: A stochastic study of the Elder (1967b) short heater problem. Wat. Resour. Res. 39 (1). Reilly, T.E., Goodmann, A.S., 1987. Analysis of saltwater upconing beneath a pumping well. J. Hydrol. 89, 169–204. Rushton, K.R., Tomlinson, L.M., 1979. Possible mechanisms for leakage between aquifers and rivers. J. Hydrol. 40, 49–65. Schincariol, R.A., Schwartz, F.W., 1990. An experimental investigation of variable density flow and mixing in homogeneous and heterogeneous media. Wat. Resour. Res. 26 (10), 2317–2329. Schincariol, R.A., Schwartz, F.W., Mendoza, C.A., 1997. Instabilities in variable density flows: stability and sensitivity analyses for homogeneous and heterogeneous media. Wat. Resour. Res. 33 (1), 31–41. Schubert, J., 2002. Hydraulic aspects of riverbank filtration – field studies. J. Hydrol. 266, 145–161. Simmons, C., 2005. Variable density groundwater flow: from current changes to future possibilities. Hydrogeol. J. 13, 116–119. Simmons, C., Narayan, K., 1997. Mixed convection processes below a saline disposal basin. J. Hydrol. 194, 263–285. Simmons, C., Narayan, K.A., 1998. Modelling density dependent flow and solute transport at the Lake Tuchwop saline disposal basin complex. Victoria. J. Hydrol. 206, 219–236. Simmons, C.T., Fenstemaker, T.R., Sharp Jr., J.M., 2001. Variabledensity groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges, invited paper. J. Cont. Hydrol. 52 (1–4), 245–275. Simmons, C.T., Narayan, K.A., Woods, J.A., Herzceg, A.L., 2002. Groundwater flow and solute transport at the Mourguong salinewater disposal basin, south-eastern Australia. Hydrogeol. J. 10, 278–295. Sophocleus, M., 2002. Interactions between groundwater and surface water: the state of science. Hydrogeol. J. 10, 52–67. Volker, R.E., Rushton, K.R., 1982. An assessment of the importance of some parameters for seawater intrusion and a comparison of dispersive and sharp-interface modelling approaches. J. Hydrol. 56, 239–250. Watson, S.J., Barry, D.A., 2001. Numerical analysis of stable brine displacements for evaluation of density-dependent flow theory. Phys. Chem. Earth 26 (4), 325–331. Watson, S.J., Barry, D.A., Schotting, R.J., Hassanizadeh, S.M., 2002. On the validity of Darcy’s Law for stable high-concentration displacements in granular porous media. Transport Porous Med. 47, 149–167. Weast, R.C., Astle, M.J., Beyer, W.H., 1984–1985. Handbook Chemistry and Physics, sixtyfifth ed. CRC Press, Boca Raton, Florida. Winter, T.C., 1999. Relation of streams, lakes, and wetlands to groundwater flow systems. Hydrol. J. 7, 28–45. Wood, M., Simmons, C.T., Hutson, J.L., 2004. A breakthrough curve analysis of unstable density driven flow and transport in homogenous porous media. Wat. Resour. Res. 40 (3). Wooding, R.A., Tyler, S.W., White, I., 1997a. Convection in groundwater below an evaporating salt lake. 1. Onset of instability. Wat. Resour. Res. 33 (6), 1199–1217. Wooding, R.A., Tyler, S.W., White, I., Anderson, P.A., 1997b. Convection in groundwater below an evaporating salt lake. 2. Evolution of fingers or plumes. Wat. Resour. Res. 33 (6), 1219– 1228.