Skip to main content Accessibility help


  • Access
  • Cited by 2



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Hydrochemical coupling of a glacial borehole–aquifer system
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Hydrochemical coupling of a glacial borehole–aquifer system
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Hydrochemical coupling of a glacial borehole–aquifer system
        Available formats
Export citation


Measurements of the electrical conductivity of subglacial water provide a useful complement to measurements of pressure and turbidity. In the summer season, fluctuations of conductivity can be attributed to changes in water transport, water provenance and subglacial residence time. These explanations are unlikely to apply during the winter season because surface melt sources are not active and the subglacial water system is predominantly unconnected. Thus, fluctuations in water conductivity during the winter months seem paradoxical. To introduce a quantitative basis for comprehending such phenomena, we develop an interpretative model of the hydrochemical interaction between a water-filled borehole and a subglacial aquifer. The electrical conductivity of water near the borehole–aquifer contact is affected not only by diffusion but also by advective transport of solute between the two reservoirs in response to pressure forcing of the system. Using records of ice strain, water pressure and electrical conductivity from unconnected boreholes in Trapridge Glacier, we demonstrate that changes in borehole geometry caused by ice-strain events provide a plausible mechanism for at least some of the observed fluctuations of electrical conductivity. Conductivity records provide information regarding advective coupling of the borehole–aquifer system that is not available from pressure records alone.


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 (Stone, 1993; Stone and others, 1993; Stone 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 Clarke (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.

Fig. 1. Postulated response to impulse pressurization. (a) Compressional strain event. Sketch (i) depicts a conductivity sensor (C) placed in an unconnected borehole. Shading is used to indicate the upward-decreasing concentration of solute in the hole. Sketch (ii) illustrates how compression of the borehole would drive fluid transport from the borehole to a subglacial aquifer and modify the concentration profile ofsolute in the borehole. Sketch (iii) illustrates how compression of the subglacial aquifer might drive the transport of solute-concentrated water from the aquifer into the borehole. (b) Uniaxial strain rate corresponding to unrecovered compression of the borehole. (c) Conductivity time series as measured by sensor C for two different cases. Flushing of mineralized water from higher in the borehole as in (a-ii) should lead to decreased conductivity detected by the sensor (solid line). Fluid flow from the aquifer as in (a-iii) should lead to increased conductivity at the sensor (dashed line). Local conductivity maxima indicate times when the strengths of advective and diffusive transport are balanced.

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 (Stone and Clarke, 1993; Waddington 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.

Fig. 2. Model geometry. The borehole is cylindrical with cross-sectional area A(t) and water-column height L(t). The aquifer is tabular and cylindrically symmetric with thickness δ aq. The radius R0 of the interface is assumed constant with time.

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




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 Bear and Verruijt (1987) and Stone (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 (Bear 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 = −Dc (e.g. De 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. Malvern, 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 (Stone 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 (Domenico 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 Clarke (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 (Patankar, 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.


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 aqc eq) |c aqc 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 (Stumm, 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 (Dobos, 1975, p. 88) and (calculated from Glasstone, 1942, p. 61; Dobos, 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 (Freeze and Cherry, 1979, p.106). Calcite solubility increases by about 25% as water temperature drops from 25° to 0°C (Langmuir, 1997, p. 206), but increased pressure increases calcite solubility only slightly (Krauskopf and Bird, 1995, p. 63). Thus, in subglacial waters, we expect equilibrium concentrations of calcite of 0.1–0.7 kg m−3. Following White (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 Clarke (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 Patankar (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 (Patankar, 1980, p. 87) which is more rigorous than upstream schemes (e.g. Bear and Verruijt, 1987, p. 322). Integration in time is performed using the Fortran DDRIV2 integration routine developed by Kahaner and Sutherland (Kahaner and others, 1989), a routine which employs the integration method of Gear (1971). This combination of finite-difference spatial derivatives and numerical integration in time is referred to as the method of lines (e.g. Schiesser, 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 Oldenborger (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 (Stone, 1993; Stone and others, 1993; Stone and Clarke, 1996) and can be adequately represented by lumped-element models (Clarke, 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 (Clarke 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 Stone 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 (Sjö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.

Table 1. Constant simulation parameters

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.

Fig. 3. Subglacial electrical-conductivity and water-pressure measurements for Trapridge Glacier during late summer 1994 and winter 1994/95. (a) Electrical conductivity recordfor sensor 94C08. (b) Water-pressure record for sensor 94P04 (in units ofpressure head). Day 225 corresponds to 23 August 1994.

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.

Fig. 4. Simulated response of electrical conductivity, interface pressure and water velocity in response to a simplified seasonal strain event. (a) Electrical conductivity at 0.0, 0.1 and 0.5 m height within the borehole (0.10 m corresponds to the sensor height for the data presented in Figure 3). Initially, advection is the dominant transport process as unmineralized water from near the glacier surface replaces mineralized water near the glacier bed. After approximately 120 days, diffusion overtakes advection as the dominant transport process. (b) Interface pressure (in units of pressure head). (c) Interface water veloci y. (d) Assumed strain-rate forcing , corresponding to unrecovered uniaxial compression of the borehole.

Table 2. Model parameters for case-study 1

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.

Fig. 5. Subglacial electrical conductivity and waterpressure measurements for Trapridge Glacier during winter 1993/94. (a) Electrical conductivity record for sensor 92C15. (b) Water-pressure record for sensor 93P07 (in units of pressure head). Day 300 corresponds to 29 October 1993.

Fig. 6. Simulated response of electrical conductivity, interface pressure and interface water velocity for a borehole in contact with a low-permeability aquifer. (a) Electrical conductivity at 0.0 and 0.1 m height within the borehole. Zero time corresponds to day number 225 (Fig. 7). (b) Interface pressure (in units of pressure head). (c) Interface water velocity. (d) Assumed strain-rateforcing , corresponding to unrecovered uniaxial expansion of the borehole.

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.

Table 3. Model parameters for case-study 2

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 (Harrison 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.

Fig. 7 Subglacial electrical-conductivity and pressure measurements for Trapridge Glacier during late summer and early winter 1994. (a) Electrical conductivity record for sensor 92C09. (b) Water-pressure record for sensor 93DP03 (in units of pressure head). (c) Vertical strain record for sensor 92VS04. Day 220 corresponds to 8 August 1994.

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 (Murray and Clarke, 1995; Waddington 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 (Sjö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).

Fig. 8. Simulated response of electrical conductivity, interface pressure and water velocity in response to a simplified impulse strain event. (a) Electrical conductivity at 0.0, 0.15 and 1.0 m height within the borehole. Zero time corresponds to day number 225 (Fig. 7). (b) Interface pressure (in units ofpressure head). (c) Interface water velocity. (d) Assumed strain-rate forcing , corresponding to unrecovered uniaxial compression of the borehole.

Table 4. Model parameters for case-study 3

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.

Fig. 9. A complicated event that cannot be simply explained using our interpretation model. (a) Englacial vertical strain ( sensor 97VS03) and subglacial water-pressure head (sensor 97P32) observed in summer 1998. (b) Subglacial water conductivity measurements for a trio array of conductivity sensors (from the bottom upward these are, respectively, 97C15, 97C16 and 97C17).

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 (Oldenborger, 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.


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.


Bear, J. and Verruijt, A.. 1987. Modeling groundwater flow and pollution. Dordrecht, D. Company.
Clarke, G. K. C. 1996. Lumped-element model for subglacial transport of solute and suspended sediment. Ann. Glaciol., 22, 152159.
Clarke, G. K. C. and Blake, E. W.. 1991. Geometric and thermal evolution of a surge-type glacier in its quiescent state: Trapridge Glacier, Yukon Territory, Canada, 1969–89. J. Glaciol., 37(125), 158169.
De Wiest, R. J. M. 1969. Flow through porous media. NewYork, Academic Press.
Dobos, D. 1975. Electrochemical data. A handbook for electrochemists in industry and universities. Amsterdam, etc., Elsevier Scientific.
Domenico, P. A. and Schwartz, F. W.. 1998. Physical and chemical hydrogeology. Second edition. Chichester, John Wiley and Sons.
Freeze, R. A. and Cherry, J. A.. 1979. Groundwater. Englewood Cliffs, NJ, Prentice-Hall.
Gear, C. W. 1971. Simultaneous solution of differential-algebraic equations. IEEE Trans. Circuit Theory, 9, 8995.
Glasstone, S. 1942. An introduction to electrochemistry. New York, D. Van Nostrand Co.
Harrison, W. D., Echelmeyer, K. A. and Engelhardt, H.. 1993. Short-period observations of speed, strain and seismicity on Ice Stream B, Antarctica. J. Glaciol., 39(133), 463470.
Kahaner, D., Moler, C. and Nash, S.. 1989. Numerical methods and software. Englewood Cliffs, NJ, Prentice-Hall.
Krauskopf, K. B. and Bird, D. K.. 1995. Introduction to geochemistry. Third edition. NewYork, McGraw-Hill.
Langmuir, D. 1997. Aqueous environmental geochemistry. Upper Saddle River, NJ, Prentice-Hall.
Malvern, L. E. 1969. Introduction to the mechanics of a continuous medium. Engle-wood Cliffs, NJ, Prentice-Hall.
Murray, T. and Clarke, G. K. C.. 1995. Black-box modeling of the subglacial water system. J. Geophys. Res., 100(B7), 10,23110,245.
Oldenborger, G. A. 1998. Electrical conductivity in glacial boreholes. (B.Sc. thesis, University of British Columbia, Earth and Ocean Sciences (Geophysics).)
Patankar, S. V. 1980. Numerical heat transfer and fluid flow. New York, Hemisphere Publishing.
Schiesser, W. E. 1991. The numerical method of lines. San Diego, CA, etc., Academic Press.
Sjöberg, E. L. and Rickard, D. T.. 1984. Temperature dependence of calcite dissolution kinetics between 0 and 62°C at pH 2.7 to 8.4 in aqueous solution. Geochim. Cosmochim. Acta, 48(3), 485493.
Stone, D. B. 1993. Characterization of the basal hydraulic system of a surge-type glacier: Trapridge Glacier, 1989–92. (Ph.D. thesis, University of British Columbia.)
Stone, D. B. and Clarke, G. K. C.. 1993. Estimation of subglacial hydraulic properties from induced changes in basal water pressure: a theoretical framework for borehole-response tests. J. Glaciol., 39(132), 327340.
Stone, D. B. and Clarke, G. K. C.. 1996. In situ measurements of basal water quality and pressure as an indicator of the character of subglacial drainage systems. Hydrol. Processes, 10(4), 615628.
Stone, D. B., Clarke, G. K. C. and Blake, E. W.. 1993. Subglacial measurement of turbidity and electrical conductivity. J. Glaciol., 39(132), 415420.
Stumm, W. 1990. Aquatic chemical kinetics. Reaction rates of processes in natural waters. New York, etc., JohnWiley and Sons.
Waddington, B. S. and Clarke, G. K. C.. 1995. Hydraulic properties of sub-glacial sediment determined from the mechanical response of water-filled boreholes. J. Glaciol., 41(137), 112124.
White, W. B. 1988. Geomorphology and hydrogeology of karst terrains. New York, etc., Oxford University Press.

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 24. 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.