Introduction
Spatial and temporal measurements of the electrical conductivity of subglacial and englacial water provide useful indicators of hydrological conditions. Interpretation of these data requires a framework for analyzing the coupling between a borehole, containing one or more conductivity sensors, and the subglacial water–sediment system with which it chemically and hydraulically interacts. The purposes of this paper are to explore this interaction and, guided by a theoretical model, to interpret field observations from Trapridge Glacier, Yukon Territory, Canada, where year-round measurements of subglacial hydraulic properties have been collected since summer 1988 (Reference StoneStone, 1993; Reference Stone, Clarke and BlakeStone and others, 1993; Reference Stone and ClarkeStone and Clarke, 1996).
In holes that are well connected to the subglacial hydraulic system, temporal variations in water conductivity observed near the ice–bed contact can largely be explained in terms of the subglacial residence time of the water. A lumped-element model to describe this situation was elaborated in Reference ClarkeClarke (1996). In holes that are unconnected or weakly connected, the explanation of temporal variations is less clear. For example, during the winter season, we occasionally observe rapid changes of water conductivity in unconnected holes. Because these fluctuations cannot be attributed to the operation of an active hydraulic system, other explanations must be sought. We propose that some of these puzzling fluctuations result from ice- or bed-deformation events that drive water exchange between the borehole and subglacial bed. During these events, advection, rather than diffusion, can become the dominant mechanism for chemical transport. Figure 1 illustrates possible hydrological and chemical responses when a sealed, water-filled borehole is subjected to sudden deformation (negative longitudinal strain rate corresponds to compression of the hole). As the hole is pressurized, downward water flow flushes comparatively unmineralized water from the borehole to the sub-glacial bed. Upon termination of a strain event, water flow ceases and ion diffusion from the bed to the borehole eventually restores chemical concentrations to their former state.
We begin by developing a physical model that describes solute transport in a borehole–aquifer system. The utility and limitations of the model are assessed by using it as a tool for interpreting electrical conductivity variations beneath Trapridge Glacier.
Physical Model
We assume that the glacier rests upon water-saturated unlithifled permeable sediment that can be regarded as a somewhat permeable aquifer, a situation that is thought to apply to Trapridge Glacier (Reference Stone and ClarkeStone and Clarke, 1993; Reference Waddington and ClarkeWaddington and Clarke, 1995). The aquifer is taken to be tabular and underlain by an impermeable boundary. A quasi-cylindrical borehole, coaxial with the z axis and having cross-sectional area A and water column height L, is in hydrological and chemical contact with the subglacial aquifer; for geometrical simplicity, the borehole is assumed to completely penetrate the subglacial aquifer (Fig. 2). Balance relations for this coupled borehole–aquifer system are developed below.
Water balance
Borehole hydraulics
The cross-sectional geometry of the borehole is described by an ellipse having a and b as its semi-axes with a aligned in the x direction (down-flow) and b aligned in the y direction (cross-flow). The cross-sectional area and water volume of the hole are therefore given by A = πab and V H = πabL, where a, b and L can vary with time. We take a(0) = R 0, b(0) = R 0 and L(0) = L0 as the initial dimensions of the water column and assume that distortions of the cross-sectional area A = πab are primarily attributed to englacial strain; thus
where and are longitudinal and transverse strain rates. The cause of variations in L depends upon whether the hole is open or sealed. For an open hole
with
where Q in and Q out are, respectively, the water discharge into the top of the hole and the discharge out of the bottom. For a sealed hole assumed to be completely water-filled, the causes of variation in L are vertical strain and freeze-down. To distinguish these, we define dL/dt = (dL/dt)ε + (dL/dt)f and write
where R f is the freezing rate. We take R f as a given forcing rather than attempting to calculate it from thermodynamics. If one chose to calculate R f, the basis of the calculation would be the expression
where K I is the thermal conductivity of ice, ΔH f is the enthalpy of fusion for ice and ∂T I(L, t)/∂z is the ice temperature gradient at the top of the sealed water column.
Combining the foregoing results, we have
for the sealed hole. In both cases, V H(0) = V 0 = πa 0 b 0 L 0. From Equation (6), it is apparent that, for an incompressible material, the contributions from vanish. For glacial ice, complete cancellation does not necessarily occur and, even when R f ≡ 0, ice strain can cause volume changes.
The mass of water contained in the borehole is simply m H = ρV H, where ρ is the density of water, and it follows that the time rate of change of water mass in the borehole is
where Q in depends upon whether the hole is open or closed. For an open hole, Q in is determined by meltwater inflow from the glacier surface; for a sealed hole, freeze-down introduces an apparent input dischage
where ρ I is the density of ice. Regardless of whether the hole is open or sealed, the output discharge is given by
where R 0 is the initial borehole radius at the hole–aquifer contact (neglecting any change with time), δ aq is the thickness of the tabular subglacial aquifer and qr (R0, t) is the radial component of the volume flux of water at the borehole–aquifer contact.
Hydrostatic pressure of the water column relative to the ambient pressure of the aquifer is the drive for water flow in an open hole. If the borehole–aquifer system is initially in equilibrium, p(r, 0) = p 0 = ρgL 0 and the equation for the pressure perturbation in the borehole is simply
For the sealed hole, the pressure perturbation is solely due to compression of the water phase and we require an equation of state for water in the form
where p is pressure, ρ R is the density at some reference value of pressure p R and β is the compressibility. Noting that
and , Equation (7) can be written as
leading to
where Q in = A(ρ − ρ I)R f/ρ.
Aquifer hydraulics
In the static case, there is no water flow and the hydraulic potential φ = p + ρgz is constant. For φ = φ 0, the pressure is given by p 0(z) = P − ρgz, where P is the pressure on the z = 0 plane defined by the top of the thin aquifer. For water to flow between the borehole and the aquifer, it is necessary that φ ≠ φ0 and p( r , t) ≠ p 0(z). The non-hydrostatic pressure contribution p′( r , t) = p( r , t) − p 0 drives water flow. For this case, the fluid potential is φ(r, t) = φ0 + p′(r, t).
The diffusion equation for transient Darcian flow is well known and will not be rederived here (see Reference Bear and VerruijtBear and Verruijt (1987) and Reference StoneStone (1993) for details). The dilatation of the solid aquifer structure can be decomposed into two components. The first is the conventional recoverable elastic component α ∂p/∂t viewed as a response to pressurization. The second component is the non-recoverable component ∂Δ s/∂t where Δs is viewed as a forced dilatation of the solid structure. The governing equation is then
where α is the aquifer compressibility, n is the aquifer porosity and q is the volume-flux density of water. Darcy’s law, with a scalar hydraulic conductivity K, gives
Assuming constant K in Equation (16) and applying the separation φ =φ 0 + p′, Equation (15) gives
where S s = ρg(α + nβ) is the specific storativity (Reference Bear and VerruijtBear and Verruijt, 1987, p. 63). If cylindrical symmetry and isotropy are assumed, Equation (17) can be reduced in polar coordinates to the single spatial variable r such that p′ = p′ (r, t) and
Hydraulic coupling
Chemical advection calculations require knowledge of the fluid velocity field inside the borehole and within the aquifer. Within the aquifer, we follow convention and calculate the average linear velocity as
where the porosity n is assumed constant and
The calculation of borehole water velocity v H(z, t) requires closer attention. For the open borehole, there is no z dependence of velocity and the magnitude is fixed by the rate of discharge of water from the borehole to the aquifer; thus we can calculate the cross-sectionally averaged fluid velocity as
In contrast, for a sealed borehole in the absence of freeze-down, the water velocity varies linearly over the length of the water column such that v H(L, t) = 0 and v H(0, t) = −Q out(t)/A(t); thus
If freeze-down occurs but no water flows from the bottom of the borehole, the resulting velocity function is linear with z but v H(0, t) = 0 and v H(L, t) = −(ρ − ρ I)R f/ρ; thus
If both processes operate simultaneously, the two effects superpose and
Borehole chemistry
The total amount N (be it mass or moles) of a given chemical species within some volume V can be expressed in terms of the concentration c (mass or moles per unit volume) as
Within the borehole, a global balance equation of the form
applies, where Φ is the volumetric rate of production of the chemical species and J is the diffusive flux (mass or moles) of the species through the bounding surface S of the volume V. In the present model, production of the chemical species is by dissolution of the basal sediments; within the borehole, Φ is assumed to vanish. This simplification is equivalent to an assumption of no suspended sediment.
Fick’s law, J = −D∇c (e.g. Reference De WiestDe Wiest, 1969, p.115), relates diffusive mass flux to a concentration gradient where c = c( r , t) is solute concentration and D is the chemical diffusivity of the medium which is assumed to be constant. From Fick’s law and Gauss’s theorem with vanishing Φ, Equation (26) reduces to
Applying Reynolds’ transport theorem (e.g. Reference MalvernMalvern, 1969, p. 211) to Equation (25) yields the material time rate of change of solute,
which, when combined with Equation (27), leads to the local form balance equation
In the borehole, we assume that c = c(z, t) and v = v H (z, t) k so that Equation (29) can be recast in a single spatial dimension as
Aquifer chemistry
The aquifer component of the system (Fig. 2) is assumed to be cylindrically symmetric (coaxial with the z axis), homogeneous and isotropic with thickness δ aq ∼ 0.10 m (Reference Stone and ClarkeStone and Clarke, 1993). Because δ eq is small, vertical variations within the porous horizon are neglected. As for the borehole, a general balance equation applies, but attention must be given to production of the chemical species and to liquid vs solid surfaces and volumes. Following the previous approach for the borehole and assuming constant porosity, the local form balance equation
describes solute concentration within the aquifer. Provided n is constant, one can substitute the porosity-adjusted diffusivity D(n) = D 0 n (Reference Domenico and SchwartzDomenico and Schwartz, 1998, p. 218) directly into Equation (31); D 0 is the diffusivity for n = 1. Recasting Equation (31) in cylindrical polar coordinates and neglecting variations with respect to z and θ gives
where qr is the radial component of water flux calculated from the pressure gradient. Within the aquifer, the volumetric rate of chemical production is non-zero and, following Reference ClarkeClarke (1996), we assume
where the c eq is the equilibrium concentration, v is the kinetic order of the reaction and sgn x = x/|x| is the algebraic sign function. The production parameter n depends upon multiple factors including the reaction rate parameter, particle size and particle geometry.
Chemical coupling
At the interface between the borehole and the aquifer, the flow field deviates from both of the idealized one-dimensional situations. To proceed, we regard the axial cylinder directly below the borehole as comprising a well-mixed control volume (Reference PatankarPatankar, 1980) into which all flow from the borehole is perpendicular to its upper surface and out of which all flow into the aquifer is perpendicular to the vertical cylindrical boundary. An expression for the time rate of change of concentration in the control volume couples the two system components. The time rate of change of N cv(t), the amount of solute within the control volume, is governed by the chemical production, the diffusive flux plus the advective flux through the upper borehole surface and the diffusive flux plus the advective flux through the vertical control volume wall. Thus
where diffusive fluxes are evaluated using Fick’s law and advective fluxes are evaluated as the product of fluid velocity and concentration at the control volume surfaces.
Synthesis
Water transport in the borehole–aquifer system is governed by the ordinary and partial differential equations describing volume and pressure evolution:
where Equations (35e) and (35f) are subject to the coupling equality . The chemical balance equations for the hole and aquifer are respectively
where, from Equation (33), Φ = −κ sgn (c aq − c eq) |c aq − c eq| v and subscripts H and aq denote concentrations within the borehole and aquifer, respectively. Finally, the two systems are subject to the coupling Equation (34).
In summary, we consider two primary ways to drive fluid motion and solute transport within the borehole-aquifer system: compressing or decompressing the hole via ice strain should lead to fluid flux into or out of the aquifer; pressurizing the hole via bed compression should lead to fluid flux into the borehole.
Scaling analysis and parameter sensitivity
The physical constants appearing in Equations (35) and (36) give a misleading impression of the number of freely adjustable parameters in the model. Dimensional analysis (see Appendix) indicates that the following five dimensionless parameters govern the behaviour of the physical system:
where we name the borehole seepage parameter, the borehole chemical-advection parameter, the aquifer chemical-diffusion parameter, the aquifer chemical-advection parameter and the aquifer chemical-production parameter. Judicious assignment of parameter values based on previous theoretical and field investigations of Trapridge Glacier reveals that the governing equations, and thus the model responses, are dominated by the advection parameters and which in turn are determined by the constants ρ, g, α, β, n, K and D 0. Of these, ρ, g and β are well-known physical constants that are not free parameters of the model, and n varies over a limited range (say 0.2 < n < 0.4). In principle, D 0 should also be well determined, but we find it necessary to treat this as an adjustable parameter.
Conductivity from concentration
For Trapridge Glacier, CaCO3 is the chemically dominant bed mineral. Hence, the major dissolved ions are those associated with calcite dissolution. If the amount of aqueous carbon dioxide is low, rate constants for the various dissolution pathways of CaCO3 are pH-dependent such that the dominant species present for acidic waters would be the anion HCO3 − and the cation Ca2+ at equal molar concentrations; for alkaline waters, CO3 2− will partially replace HCO3 − in solution (Reference StummStumm, 1990, p. 433). The former situation would be most suitable for meteoric waters and possibly borehole water at some distance above the bed; the latter might apply to longresidence-time water at and within the bed.
Conductivity σ of an ionic solution can be written as
where Λ i is the equivalent molar conductivity of the ith ion at infinite dilution and Ci is the molar concentration of the ith ion. Because molar concentrations of cations and anions must balance,
where C is the molar concentration of either species. At 0°C, ΛCa 2+ = 61.6 × 10−4 S m2 mol−1 (Reference DobosDobos, 1975, p. 88) and (calculated from Reference GlasstoneGlasstone, 1942, p. 61; Reference DobosDobos, 1975, p. 87).
Thus far, the development of the physical model is insensitive to whether concentrations are expressed in terms of mass or moles. We regard mass concentration as the more tangible quantity, and henceforth will follow this preference. Since one mole of calcium carbonate yields one mole of each of HCO3 − and Ca2+, the molar concentration C can be expressed in terms of the mass concentration of CaCO3 and its molar mass of approximately 0.10 kg mol−1. Similarly, the ionic molar conductivities can be combined and converted to a mass conductivity for calcite. Thus, the conductivity of this simplified calcareous solution can be expressed as
where Λ m = 0.0816 S m2 kg−1 is the equivalent mass conductivity for CaCO3 and c(r, t) is the mass concentration of the calcium carbonate solute as a function of space and time. For alkaline subglacial waters, the contribution of CO3 2− will increase the equivalent conductivity up to a maximum of Λ m = 0.1336 S m2 kg−1.
With respect to the chemical-production Equation (33), the equilibrium concentration of calcite for ground-waters at 25°C and 1 bar (105 Pa) pressure can be considered to fall between 0.1 and 0.5 kg m−3 for partial pressures of carbon dioxide between 10−3 and 10−1 bar (Reference Freeze and CherryFreeze and Cherry, 1979, p.106). Calcite solubility increases by about 25% as water temperature drops from 25° to 0°C (Reference LangmuirLangmuir, 1997, p. 206), but increased pressure increases calcite solubility only slightly (Reference Krauskopf and BirdKrauskopf and Bird, 1995, p. 63). Thus, in subglacial waters, we expect equilibrium concentrations of calcite of 0.1–0.7 kg m−3. Following Reference WhiteWhite (1988, p.144), we take v = 2 as an appropriate value for CaCO3 dissolution in many natural situations and we assign k based on the estimates of Reference ClarkeClarke (1996).
Numerical analysis
Equations (35) and (36) together with Equation (34) comprise an initial value problem in time and a non-homogeneous boundary value problem in space. The partial differential equations are discretized in space based on the control volume integration method of Reference PatankarPatankar (1980) applied to a logarithmically stretched, node-centred grid in r and a regularly spaced, though time-varying, grid in z. To maintain accuracy we employ a specialized scheme to model advective transport (Reference PatankarPatankar, 1980, p. 87) which is more rigorous than upstream schemes (e.g. Reference Bear and VerruijtBear and Verruijt, 1987, p. 322). Integration in time is performed using the Fortran DDRIV2 integration routine developed by Kahaner and Sutherland (Reference Kahaner, Moler and NashKahaner and others, 1989), a routine which employs the integration method of Reference GearGear (1971). This combination of finite-difference spatial derivatives and numerical integration in time is referred to as the method of lines (e.g. Reference SchiesserSchiesser, 1991).
The concentration signature of borehole-aquifer interactions is most pronounced near the interface between the system components, and we expect remote boundary conditions to have a negligible influence. Within the aquifer, at the maximum radial distance of the computational grid, we assume zero pressure perturbation and vanishingly small diffusive transport; advective transport will naturally approach zero owing to the cylindrical geometry. At the top of the borehole water column, the boundary condition is tailored to suit the active processes.
As an initial test of the theoretical development, we used the numerical model to interpret laboratory benchtop borehole–aquifer experiments described in Reference OldenborgerOldenborger (1998). Although some degree of qualitative agreement was achieved, the effort was undermined by unsatisfactory experimental technique and requires repetition before detailed discussion is warranted.
Application to Field Data
Long-term measurements of subglacial water pressure, electrical conductivity and ice strain for Trapridge Glacier provide a substantial data archive from which the present study draws. Few of the available records demonstrate a self-evident coupled response of borehole water pressure and electrical conductivity, but, for those that do, it is useful to distinguish between summer phenomena and winter phenomena. By summer phenomena, we refer to variations that are associated with the flow of surface meltwater through a connected water system. The causes of these variations are reasonably well understood (Reference StoneStone, 1993; Reference Stone, Clarke and BlakeStone and others, 1993; Reference Stone and ClarkeStone and Clarke, 1996) and can be adequately represented by lumped-element models (Reference ClarkeClarke, 1996). Winter phenomena, less obviously associated with meltwater flow through an established subglacial drainage network and more challenging to explain, are the principal concern of this paper. Trapridge Glacier has a subpolar thermal regime with temperature increasing from roughly −5°C near the surface to the melting point near the bed (Reference Clarke and BlakeClarke and Blake, 1991). The implication for our study is that holes drilled to the bed will rapidly seal themselves from surface water input; thus Qin = 0 within hours of completing a drillhole. Four case histories have been selected for special attention; of these, three are readily comprehended using the interpretation model.
Model initialization
Table 1 summarizes the values of physical constants and model parameters that are common to all of the following model simulations. We have previously remarked that aquifer compressibility is a free parameter, but, consistent with Reference Stone and ClarkeStone and Clarke (1993), we have fixed this value at α = 1.0 × 10−8 Pa−1. Initial concentration profiles for the system are determined from steady-state simulations. An exponential concentration gradient in the borehole will result from upward chemical diffusion from the aquifer in competition with a steady downward water flux associated with decreasing borehole volume and a permeable bed. Equation (35) indicates that for constant pressure, volumetric flow rate will exactly balance volumetric changes regardless of K. The time required to reach steady state and the magnitude of pressure increase will depend on K; the steady-state concentration profile, however, will be a function only of the volumetric strain rate and the chemical diffusivity. Using the measured diffusion coefficient for calcite of D0 = 1.66 × 10−10 m2 s−1 (Reference Sjöberg and RickardSjöberg and Rickard, 1984), volumetric strain rate can be adjusted to yield a steady-state concentration profile in the borehole that matches the observed data immediately preceding an impulse event.
For simplicity, we invoke a uniaxial strain rate and exclude downhole freezing. For a sealed hole with no bed strain, any processes that combine to yield the same resultant hole dilatation will lead to essentially equivalent results. In practice, the strain rate is applied to a system for which the borehole is at some ambient concentration c a and the aquifer is at some equilibrium concentration c eq (Table 1). The ambient concentration is that of continental rainwater, and the equilibrium concentration is set in accordance with observed electrical conductivities in excess of 90 mS m−1. This necessary equilibrium concentration is greater than the expected value of 0.7 kg m−3 but is not of major concern since the assumed conductivity-concentration relationship is a simple proportionality. High observed conductivities may be indicative of a transition from acidic to alkaline waters for which lower concentrations would result in higher conductivities.
Another point that deserves comment is that we take the initial state of the borehole to be the unstrained state. Thus, the initial borehole pressure is dictated by the initial height of the water column. Here, and in subsequent examples for sealed holes, we adjust the initial water-column height L0 to match the initial pressure in the hole. Thus, the choice of L0 is unrelated to the actual length of the water column in the frozen-over borehole. For a uniaxial strain rate applied to a cylindrical borehole, the fractional volume change, and hence the induced pressure variation, is independent of L. Therefore, this simplifying assumption has no effect on the simulation results provided that the true hole length exceeds the thickness of the chemically altered water layer.
Month- to seasonal-scale events
Case 1. Seasonal-scale pressure event
As our first example, we consider the conductivity record for sensor 94C08 and the pressure record for sensor 94P04, both presented in Figure 3. The sensors are located in the same borehole which was hydraulically connected at the time of drilling and penetrated 64 m of ice. The conductivity sensor was installed so that its electrodes were 0.10 m above the bed. The hole was drilled in July 1994 and, like all holes drilled inTrapridge Glacier, sealed by freezing within 1 day. Both records commence in late August 1994, near the end of the summer melt season, and terminate in late March 1995, before the onset of summer melt. The main feature of the conductivity record is a pronounced conductivity drop from ∼6 mS m−1 to ∼2 mS m−1 beginning at day 275 and lasting roughly 10 days. The main feature of the pressure record is a slow increase in pressure head from ∼60 m to ∼82 m followed by a slow decrease to ∼70 m. The pressure maximum coincides with the conductivity drop. Throughout the observation interval, borehole water pressure exceeded the ice flotation pressure at this site. Examination of the vertical ice-strain record (not shown) for a nearby site reveals a monotonic increase in vertical strain, which is consistent with monotonic longitudinal compression. There is no obvious indication of the occurrence of an impulse ice-strain event; thus we cannot point to a glacier-scale strain event as the cause of the broad peak in the pressure record. Nevertheless, the reality of the water-pressure peak cannot be dismissed and we postulate that it is an expression of near-borehole ice deformation processes.
Table 2 and Figure 4 present the fitted model parameters and results of our attempt to simulate the essential features of Figure 3. In our judgment, the model adequately reproduces the quantitative and qualitative features of the field data. For simplicity, we invoke a uniaxial strain rate and exclude downhole freezing. The postulated ice deformation event is represented by a Gaussian strain-rate function (Fig. 4d) that corresponds to approximately 1 mm of unrecovered down-flow compression of the borehole. Width of the strain-rate function is adjusted to yield a simulated pressure signal (Fig. 4b) that has similar onset duration to that of the observed pressure impulse (Fig. 3b). Hydraulic conductivity, in conjunction with strain-rate amplitude, dictates the magnitude and decay of the simulated pressure perturbation and also the magnitude and duration of the conductivity response.
The most conspicuous discrepancy between the observations and the simulation is that the field example exhibits a rapid drop in electrical conductivity coincident with the pressure maximum but occurring well after the onset of the pressure increase. The simulated system, however, exhibits a nearly immediate conductivity response to the pressure forcing. If the conductivity record is ignored, the pressure response can be simulated very well. The simulations (Fig. 4) represent a balance between the accuracy of the simulated pressure and conductivity responses (with special attention paid to matching the onset of the pressure signal). The example illustrates the coupling between variables: any number of scenarios may produce a given pressure response, but a reasonable scenario must also honour the conductivity record which acts as a proxy for advective solute transport.
Scaling analysis of the governing equations (Appendix) confirms that K, D0 and α are the primary controls on the model system, and an exploration of a range of model behaviours yields further insight into model parameters. By isolating hydraulic conductivity, it can be seen that reduction of K increases the delay time of the conductivity response with respect to the pressure signal. However, the same reduction in K results in a conductivity response of increased duration which is not supported by the field data. Similarly, increasing K results in a conductivity response of shorter duration, but the delay time with respect to the pressure signal is decreased. Sensitivity analysis shows that values of D 0 ≤ 10−9 m2 s−1 yield acceptable results. For D 0 > 10−9 m2 s−1, the electrical conductivity of the borehole water at z = 0.10 m returns too rapidly to pre-event values.
Case 2. Effect of bed permeability
Correlation between borehole water pressure and electrical conductivity records is expected only if pressure variations produce substantial water-flow rates. Lack of correlation is expected for a low-permeability substrate. The observations (Fig. 5) and the simulation results (Fig. 6) for this case-study illustrate that the degree of coupling of the unconnected borehole–aquifer system, and thus the magnitude of the conductivity variation, is highly sensitive to bed permeability.
Figure 5b shows a 10 m fluctuation in the pressure head indicated by sensor 93P07 occurring over an interval of 25 days in winter 1993/94. The hole was drilled in July 1993 and was hydraulically unconnected at the time of drilling. The conductivity record for sensor 92C15 (Fig. 5a), installed with its electrodes 0.11 m above the bed and in the same hole as the pressure sensor, reveals that only minor fluctuations in conductivity accompany this significant pressure event.
Model parameters are listed inTable 3. Again, the column height is dictated by the pressure record, and the chemical diffusivity is taken as the measured value for calcite (sensitivity to variations in D is found to be similar to that of the case 1 simulations). Conditions differing from case 1 are the initial concentration profile, the hydraulic conductivity and an assumed strain event (Fig. 6d) corresponding to down-flow expansion of the borehole. Compared to case 1, there is an order-of-magnitude reduction in the assumed value for hydraulic conductivity of the aquifer. This reduced hydraulic conductivity value is the limiting value for which the pressure impulse does not result in significant fluctuation of electrical conductivity within the hole. The conductivity record acts to constrain simulations such that the pressure record must be matched subject to very little exchange of water between the hole and aquifer. For example, a more permeable aquifer and a larger strain event could yield the same pressure signal but would be accompanied by significant flow which would be reflected by the conductivity record. In this fashion, the simulations can establish an upper bound on K for the glacial substrate.
Impulse events
Case 3. Pressure impulse
Conductivity sensor 92C09 was installed with its electrodes 0.155 m above the bed in the same hole as pressure sensor 93DP03. The hole was drilled in July 1993 and was hydraulically connected at the time of drilling. A vertical-wire ice-strain meter 92VS04 (Reference Harrison, Echelmeyer and EngelhardtHarrison and others, 1993) was installed at a nominal depth of 15 m in a nearby shallow hole drilled in July 1992. Observations spanning late summer and early winter of 1994 are presented in Figure 7.
The most prominent feature of the records for the three sensors is the event occurring on day 238 (26 August 1994). The event is marked by a rapid drop in electrical conductivity (Fig. 7a) and two distinct positive pressure pulses (Fig. 7b). The double-peak structure of the pressure pulse is also apparent in the fine structure of the conductivity variation. Close examination of the time series for conductivity and pressure reveals that the conductivity decrease occurred over an interval of ∼4 hours, whereas the pressure increase occurred in <20 min (the winter-mode sampling rate of our data loggers). The observed step-like increase in vertical ice strain (Fig. 7c) indicates that vertical elongation of glacier ice coincided with an increase in borehole water pressure. This combination of circumstances could be achieved if the volume change in the borehole was dominated by horizontal compression.
Results of our attempt to simulate the observed features of Figure 7 are presented in Table 4 and Figure 8. The magnitude of the pressure pulse in the model corresponds to a 20 m increase of hydraulic head and has similar magnitude to the observed pressure pulse. The magnitude of the observed conductivity fluctuation is also well matched by the model. In order to accommodate rapid flushing, successful modelling requires a relatively permeable aquifer with a hydraulic conductivity on the order of 10−8 m s−1, which is within the range expected for Trapridge substrate (Reference Murray and ClarkeMurray and Clarke, 1995; Reference Waddington and ClarkeWaddington and Clarke, 1995). Similarly, to realize a quick return to equilibrium, a solute diffusivity on the order of 10−6 m2 s−1 is required, a value far greater than the measured diffusion coefficient for calcite of 1.66 × 10−10 m2 s−1 (Reference Sjöberg and RickardSjöberg and Rickard, 1984). Resulting discharge rates are high (Fig. 8c), and disturbance of bed sediment is conceivable. Therefore, plausible explanations for the high diffusion coefficient might involve suspended sediment within the borehole. If sediment is present in suspension, chemical dissolution will be active within the water column rather than localized in the sediment layer. Interpretation of the observed conductivity profile solely in terms of upward diffusion from the bed would therefore lead to an erroneously high estimate of chemical diffusivity. While dimensional analysis suggests a relative insensitivity to chemical production, increasing the production parameter by three orders of magnitude does provide a slightly increasing asymptote for diffusive equilibrium at late times (not shown).
Case 4. Complex phenomena
Figure 9 shows records from pressure sensor 97P32 and a vertical array of three conductivity sensors 97C17 (top), 97C16 (middle), 97C15 (bottom) placed in the same hole drilled in July 1997. The hole was hydraulically unconnected at the time of drilling. The conductivity sensors were placed with their electrodes at 0.10, 0.25 and 0.40 m above the bed. The record for vertical-wire strain meter 97VS03, installed at a nominal depth of 15 m in a nearby hole, is also shown.
Striking features of the observations presented in Figure 9 are the simultaneous increase in electrical conductivity on 30 June 1998 recorded by the three conductivity sensors, the accompanying increase in vertical ice strain and the curious staircase-like decrease in water pressure. Closer examination reveals that there is no apparent water-pressure fluctuation occurring at the time of the dramatic increases in water conductivity and that the tendency for increase in vertical strain is arrested, even slightly reversed, at the onset time for the electrical conductivity events. Rapid increase in the conductivity of borehole water suggests inflow of mineralized water from the bed consistent with volumetric expansion of the sealed borehole. It appears that the conductivity signal is largely driven by variations in ice strain that are not immediately expressed as variations in water pressure; the absence of a strong association between water-pressure variations and conductivity fluctuations is consistent with a bed that is comparatively permeable.
Note that before the onset of the conductivity event, the electrical conductivity profile in the borehole is consistent with the expectation that solute concentration decreases with distance from the bed. Sensor 97C15, 0.10 m above the bed, indicates a substantially greater water conductivity than sensor 97C17 located 0.40 m above the bed. The electrical conductivity indicated by sensor 97C16, 0.25 m above the bed and equidistant from sensors 97C15 and 97C17, is near that for the uppermost sensor, suggesting the presence of an exponentially fading chemical boundary layer having a characteristic length scale of ∼0.2 m. Following the event, all conductivity sensors show an increase in conductivity. However, the middle sensor shows the greatest increase and, in fact, increases to a value exceeding that indicated by the lowermost sensor. This observation is also noted qualitatively in an experimental environment involving a benchtop borehole-aquifer system (Reference OldenborgerOldenborger, 1998). For both situations, it seems that the chemical stratification has been disturbed by factors that are not included in the interpretation model. We have not attempted to simulate the observations for case 4.
Concluding Remarks
Translating subglacial sensor measurements into convincing explanations that involve the operation and interaction of sub-glacial physical processes remains a challenging task. In this contribution, we have provided a starting theoretical framework for the interpretation of variations in electrical conductivity in conjunction with pressure records from unconnected boreholes. Several features of the subglacial water system have been highlighted: (i) A chemical boundary layer forms at the contact between unconnected water-filled boreholes and the glacier bed. (ii) Ice deformation processes can disturb this boundary layer and, by doing so, generate a measurable change in the electrical conductivity indicated by sensors placed in boreholes. (iii) Nearly coincident variations in borehole water pressure and electrical conductivity can occur if there is water exchange between the borehole and the sub-glacial environment, but in the limiting cases there is either a pressure signal with no conductivity fluctuation (impermeable bed) or a conductivity signal with no pressure signal (highly permeable bed). In this sense, records of electrical conductivity provide information on the advective communication of the borehole–aquifer system, information not provided by pressure and strain records alone.
Although the gross features of hydrochemical interactions between water-filled boreholes and a subglacial aquifer can be reproduced by a physically based simulation model, nature is more complicated than our model, and some observations remain unexplained. For example, the observations of case 4 defy explanation by a simple advection-diffusion model and suggest the importance of more complex diffusion and dissolution processes, possibly activated by turbulent water flux and entrainment of suspended sediment. The extent of these difficulties is only revealed by installing a vertical array of conductivity cells near the glacier bed and it will require substantial effort to clarify the details of these mechanisms.
Acknowledgements
This research was funded by the Natural Sciences and Engineering Research Council of Canada. Many years before publishing the construction details of his ingenious vertical-wire strain meters, W. D. Harrison kindly provided us with his design and assembly instructions. Since 1988 we have installed these devices in Trapridge Glacier, and we take this opportunity to acknowledge their inventor. We are grateful to Parks Canada staff for permitting fieldwork to be conducted in Kluane National Park. Base camp support was provided by the Arctic Institute of North America. We thank A. Fountain, M. Sturm, M. Tranter and J. Walder for helpful criticisms.
Appendix: Scaling Analysis
For scaling analysis of the coupled borehole–aquifer system, we shall consider the simplified physical model as directly applied to the Trapridge case-studies. For these cases, the hole is sealed and there is zero freeze-down such that Q in = 0. As well, the forced dilatation term Δs, originally included for generality, is not employed. Dimensionless variables, indicated by an asterisk, are defined by the following scalings:
The scalings of Equations (A1a–A1q) lead to the following non-dimensional equations:
subject to the initial conditions a*(0) = 1, b* (0) = 1, L*(0) = 1, , specified and with
and borehole–aquifer coupling conditions
In Equations (A2g) and (A2h) the dimensionless parameters are isolated by braces and have the following interpretations: the borehole seepage parameter
characterizes the relative magnitude of water exchange between the borehole and the aquifer; the borehole chemical-advection parameter
characterizes the importance of advective transport of solute in the borehole; the aquifer chemical-diffusion parameter
characterizes the relative effectiveness of diffusive transport of solute in the aquifer; the aquifer chemical-advection parameter
characterizes the importance of water velocity within the aquifer; the aquifer chemical-production parameter
characterizes the importance of chemical production within the aquifer.
The dimensionless parameters defined above were evaluated for each of the three case-studies and these values are included in Tables 2–4. For all of our case-studies the values of the production parameter and the diffusion parameters and were found to be small, whereas those of the advection parameters and were large. In terms of the dimensionless governing equations, the implication is that , and have a negligible influence on the model response. Thus from Equations (A6) and (A8) it is apparent that the most influential physical properties are ρ, g, α, β, n, K and D 0. Of these, ρ, g and β are well-determined physical constants that are not subject to tuning, and the plausible range for n is small. Thus the physical model has essentially three free parameters. We have not expected the chemical diffusivity D 0 to be a free parameter of the model but were forced to accept this outcome. As explained in the discussion of case-study 3, the presence of suspended sediment could produce a high effective value of diffusivity.