Estimation of subglacial hydraulic properties from induced changes in basal water pressure: a theoretical framework for borehole-response tests

Abstract

In this paper we develop a theoretical model describing water motion in a coupled borehole-subglacial flow system. The theory applies to basal drainage systems having multiple and extensive interconnected flow paths. Within this domain it encompasses a broad range of flow regimes, from laminar Darcian flow in a thick permeable unit to turbulent sheet flow in a very thin layer. Important terms in the model are highlighted by recasting the problem in dimensionless form. The non-dimensional formulation indicates that there are four free parameters in the coupled system. These parameters characterize skin friction in the borehole, and diffusion, transmissivity and turbulent transport in the subglacial flow layer. Dimensionless results show that, under most circumstances, the effects of skin friction in the borehole are negligible. Diffusion, transmissivity and especially turbulent transport in the basal layer are found to influence subglacial water flow strongly. We use our model to predict fluctuations of borehole-water levels that result from different types of disturbances. We show how this framework can be used to estimate subglacial hydraulic properties by comparing model results with data collected during field experiments on Trapridge Glacier, Yukon Territory, Canada in 1989 and 1990.

1. Introduction

It is widely recognized that subglacial water-flow systems play a key role in regulating the motion of glaciers and ice streams. However, understanding of the morphologies and hydraulic properties that characterize subglacial flow systems remains incomplete. Perhaps nowhere is the influence of the basal drainage system on glacier motion better demonstrated than in the case of a glacier surge. This is embodied in the review by Raymond (1987): “A pivotal question in the surge mechanism concerns the cause of buildup of stored water and high basal water pressure.. .. Major questions concerning how water flows in a distributed system of basal cavities or other passages and how this water affects sliding need to be addressed.”

We seek to characterize subglacial flow conditions by estimating hydraulic parameters that regulate water flow at the bed. In what follows, we develop a theoretical model of water motion in a coupled borehole-subglacial flow system. Our approach does not follow the traditional view that clean ice overlies rigid bedrock. Instead, we present the model with the idea that glacier ice can rest on unlithified permeable sediments. This picture of the glacier bed is suggested by the observations of Boulton and Jones (1979), Clarke and others (1984), Alley and others (1986), Blankenship and others (1986), Boulton and Hindmarsh (1987) and Engelhardt and others (1990b). The description we present is cast in terms of ground-water flow through a saturated porous medium (Boulton and Jones, 1979; Shoemaker, 1986; Clarke, 1987); however, our mathematical characterization also allows consideration of other distributed flow regimes: flow as a sheet or thin film between ice and bedrock (Weertman, 1957, 1964; Lliboutry, 1968; Kamb, 1970), flow through interconnected water-Filled cavities (Walder, 1986; Kamb, 1987) and distributed channelized flow, either through an ice-bedrock network (Weertman, 1972) or over basal sediments (Boulton and Hindmarsh, 1987; Walder and Fowler, 1989). For interconnected cavities or distributed channel networks, physical descriptions of the actual flow systems are not given by our model; instead, these systems are represented, in a general way, as permeable units having equivalent hydraulic characteristics. Although the theory is relevant to many subglacial flow regimes, our model does not apply to certain drainage configurations. Flow through a single channel incised upward into basal ice (Röthlisberger, 1972; Shreve, 1972) or downward into bedrock (Nye, 1973) is not described by our model because such configurations do not constitute widespread interconnected sets of flow paths.

To focus our discussion, basal water flow will be represented as occurring beneath a glacier that rests on a saturated substrate through which a significant amount of water is transmitted. In hydrogeologic terminology, an aquifer is a layer, formation or group of formations of geologic material, saturated with water, and having a high degree of permeability (Marsily, 1986). Thus, the glacier bed we consider will be referred to as a subglacial aquifer.

In ground-water studies, well-response tests are commonly used to evaluate aquifer properties. These tests consist of disturbing the equilibrium hydraulic head in the aquifer by changing the amount of water in a well and observing the equilibrium recovery either in the same well or in nearby wells. In the simplest case, a parcel of water is suddenly removed from an equilibrated well and water level in the well is monitored until the predisturbed value is regained. In an analogous technique, boreholes through a glacier that penetrate a confined basal aquifer can be used to observe fluctuations of hydraulic head within the confined layer. Surprisingly, application of this technique has received little attention in glaciological studies. Hodge (1976) measured water-level drops in boreholes while drilling through South Cascade Glacier, Washington. He also reported inducing damped oscillations by suddenly displacing water in the boreholes with the drill tip. Engelhardt (1978) described borehole water-level fluctuations in Blue Glacier, Washington, that were induced by pumping additional water into the borehole. C. Smart (personal communication) has measured borehole-drainage rates and the responses of boreholes following episodes of pressurization and release. Engelhardt and others (1990a) have observed the rates of water-level lowering as boreholes reached the base of Ice Stream B, Antarctica.

We begin by developing a theoretical framework for borehole-response tests. Following this, we consider important physical aspects of the model as highlighted by dimensional analysis. We use a dimensionless formulation to predict the responses of coupled borehole-subglacial water-flow systems having different hydraulic characteristics. Next, we demonstrate how the theoretical description is used to estimate hydraulic parameters by comparing model results with data collected from field experiments on Trapridge Glacier, Yukon Territory, in 1989 and 1990. We conclude with a discussion of the model, its limitations and generalization to a wide range of subglacial flow regimes.

Theory

Types of Disturbance

In the following sections we present a model that simulates the response of a coupled borehole-aquifer system to three different types of disturbance, corresponding to field observations that we have made. One type of response occurs when a water-filled borehole is suddenly opened to the basal aquifer. Such a situation arises, for instance, when the bed is reached by hot-water drilling and the borehole becomes connected to the subglacial flow system. We refer to this process as a connection–drainage disturbance. Another type of response is observed when the water level in an open borehole connected to the basal aquifer is displaced from its equilibrium position and allowed to recover. The simplest way of inducing this type of response is to lower a sealed pipe into the borehole, wait a sufficient amount of time for the water level to re-equilibrate, and then quickly remove the pipe. This procedure is widely used in ground-water studies and is commonly referred to as a slug test. A slightly more complicated test involves pressurizing the air in a borehole that is sealed at the top and observing the response when the pressure is suddenly released. Because it is difficult to seal the borehole perfectly, in practice the pressure is usually released before the water level has stabilized. We refer to this type of disturbance as a packer test.

High flow velocities are rarely encountered in most sub-surface hydrologic applications. Thus, energy losses in the well and the effects of turbulent transport in the aquifer are typically neglected in standard ground-water models. During a connection-drainage disturbance, however, a substantial volume of water can drain from the borehole in a short period of time. Under these conditions, water flow may be turbulent, both in the borehole and in the subglacial aquifer near the borehole. Because flow velocities can be significant in some of the situations that we will consider, turbulent effects are included in the following development.

Motion of Water in the Borehole

Water flow in a borehole of radius r_w is described by the parameters shown in Figure 1. Following standard procedures for well-response analyses (e.g. Cooper and others, 1965; Kamp, 1976; Kipp, 1985), a well screen or filter of radius r_f is included as part of the geometrical description. Such filters are not actually used in our field studies, but disturbances at the bottom of the borehole due to drilling likely result in excavations that are conveniently represented by a filter. To facilitate treatment of water flow between the borehole and the subglacial aquifer, it is assumed that the filter fully penetrates the permeable layer and that flow into or out of the borehole is uniformly distributed across the entire aquifer thickness.

Fig. 1. Model geometry showing parameters and variables that regulate flow in the vicinity of a borehole. A basal layer of uniform thickness b and hydraulic conductivity K₁ is assumed to rest on an aquitard having conductivity K₂ ≪ K₁. Water from a borehole of radius r_w is introduced into the flow layer at the filter radius r_f. Disturbances, causing water to enter or leave the basal aquifer, are produced by suddenly changing the borehole water-column height h_w or by changing the pressure acting at the top of the borehole. In the model, the pressure change is expressed as the height ht of a water column that would produce an equivalent fluid pressure at its base. In the absence of flow, the piezometric surface represents the hydraulic head h_B in the basal aquifer.

For packer tests, water-level fluctuations are induced by sealing the borehole, pressurizing it, then suddenly releasing the pressure. The pressure rise p during these tests can be expressed as the height h_T of a water column that would produce an equivalent fluid pressure at its base: p = ρ_wgh_T, where ρ_w is water density and g is the acceleration due to gravity. Thus, for packer tests, the disturbance can be represented as a downward-directed surface force acting on the top of the water column:

(1)

For connection-drainage disturbances and slug tests, there are no additional pressure forcings; hence, h_T = 0 in these cases.

When water in the basal aquifer is stationary, the borehole-water level represents the static or piezometric head in the immediate vicinity of the borehole bottom. Under these conditions, piezometric head is equivalent to hydraulic head (Marsily, 1986, p. 51). When water flows in the basal aquifer, hydraulic head increases as kinetic energy is gained by the fluid; the increased hydraulic head is manifested as a rise in the borehole-water level. In either case, hydraulic head h_B(r,t) in the basal aquifer acts to regulate the borehole-water level. Assuming that the water-column height is much greater than the aquifer thickness b, and that h_B(r,t) is uniform over the distance of the filter radius, hence h_B(0,t) = h_B(r_f,t), the upward-directed surface force supporting the water column is

(2)

If τ₀ is the frictional shear stress acting at the wall of the borehole, the force exerted by the borehole wall on the water column is

(3)

where h_w is the height of water above a point at the bottom of the borehole. In fluid mechanics, it is customary to define the skin-friction coefficient C_f in a long pipe by

(4)

where ν is the mean fluid velocity in the pipe (Kay and Nedderman, 1985, p. 170). Thus, the frictional force acting on the water column can be expressed as

(5)

where sgn (x) is the algebraic sign function (Bracewell, 1978, p. 61); this is included to account for the fact that the frictional force acts in a direction opposite to the mean velocity. The downward-directed gravitational body force acting on water in the borehole is

(6)

If we assume that ν represents a uniform water velocity over the borehole cross-section, and that water compressibility is a negligible component of the rate of change of water-column height, we have

(7)

where h(t) = h_w(t) + b. In this case, the momentum Ρ of water in the column can be written as

(8)

The total force acting on the water column will be equal to the rate of momentum outflow across the bottom surface of the borehole plus the rate of change of momentum in its interior. Thus,

(9)

where terms involving water-column velocity have been rewritten in accordance with Equation (7). The first term is the time rate of change of momentum of the water column, the second is the flux of momentum across the borehole base and the righthand-side terms correspond to the sum of forces expressed by Equations (1), (2), (5) and (6).

In general, the skin-friction coefficient in Equation (9) depends on both the Reynolds number Re and the wall roughness (Prandtl, 1952). For a pipe of internal diameter d, the Reynolds number expresses the ratio of inertial to viscous forces: where n is the dynamic viscosity of the fluid and ν is its mean velocity. Because our present concern is only with straight vertical conduits in ice, we will assume that wall roughness can be neglected. With this assumption, the relationship between the skin-friction coefficient and Reynolds number is of the form

(10)

where a and γ are positive constants. For Reynolds numbers ≲2000, the flow in a smooth pipe will be laminar with a = 16 and γ = 1; otherwise, for values of Re up to about 10⁵, the flow will be turbulent with a = 0.079 and γ = 0.25 (Kay and Nedderman, 1985). Boreholes through glaciers usually have small radii compared to their lengths; a typical borehole through the ice of Trapridge Glacier has a radius of 0.05 m and a length of 70 m. Thus, we consider a borehole to be a long smooth pipe. If water velocity in the borehole is such that Re ≲2000, the laminar form of Equation (10), combined with the momentum conservation expression (9), leads to

(11)

Equation (11)is the differential equation that describes the height of water h in the borehole at any rime t. In writing the final form of this expression, we have made use of two previous assumptions; namely, ν = dh/dt and h_w, ≫ b, so that h ≈ h_w. Also note that sgn(dh/dt) does not appear in the frictional term, because direction is given by the velocity, which now appears to the first power. The skin-friction coefficient will be slightly underestimated by Equation (11) when 2000 ≲ Re ≲ 10⁵. As we will show, however, skin friction is a minor component of the borehole-aquifer system. Thus, switching from the laminar to the turbulent form of Equation (10) represents a small correction to an insignificant term. For these reasons, we will simplify our model and not consider a separate friction coefficient for turbulent flow in the borehole.

Expressions describing displacement of the water level from an initial position of equilibrium in a coupled well-aquifer system have been derived by Cooper and others (1965) and by Kipp (1985). Like Equation (11), these expressions were based on conservation of momentum. Kabala and others (1985) used the framework of Cooper and others (1965) in numerical modelling of the responses of well-aquifer systems to sudden changes of water levels. In their developments, Cooper and others (1965) and Kipp (1985) neglected loss of momentum due to skin fricüon. It was shown by Kamp (1976) that this term is important only in cases where the well radius is very small, or when the oscillations are slowly damped, or if the initial displacements are large compared to the well radius. Because some of our observations include large initial displacements, we will retain all the non-linear terms in Equation (11).

Motion of Water in the Aquifer

Movement of water in the basal aquifer is assumed to obey the following balance equation:

(12)

where q_j is the fluid-volume-flux vector (specific discharge), n is porosity and α and β are constant compressibility coefficients for the porous medium and the fluid, respectively. Equation (12) is based on conservation of fluid and solid mass, and is a standard expression appearing in many developments of the equation of transient ground-water flow (Jacob, 1940; Cooper, 1966; Wiest, 1966; Gambolati and Freeze, 1973; Bear and Verruijt, 1987). The assumptions made in the development of Equation (12) are as follows: (i) displacements of solid grains occur only in the vertical direction; (ii) thickness and density of the overlying material are constants and atmospheric pressure fluctuations are negligible; (iii) temporal changes in fluid pressure are much greater than the rate at which pressure gradients are advected by motion of the solid skeleton; (iv) the fractional change in the fluid volume flux is much greater than the fractional change in fluid density. We recognize that the first assumption is valid only if horizontal deformations of the sediments can be disregarded with respect to the vertical deformations. Verruijt (1969) has explored this assumption and has shown that in some cases, though not all, the errors introduced by this assumption will be negligible. The second assumption is justified by the observation that time-scales of over-burden-pressure variation are much greater than the duration of a response test. Assumption (iii) will be of questionable validity in some cases, as mentioned by Gambolati and Freeze (1973). However, for the situation with which we shall be concerned — purely radial fluid flow—assumption (iii) will automatically be satisfied if assumption (i) is true; under these conditions, the fluid and solid skeleton velocities are orthogonal and there will be no advection of pressure gradients due to skeleton displacements. The final assumption also seems reasonable because the fluid we are concerned with is water, which is only very slightly compressible; thus, we expect that the fractional change in fluid density will be extremely small.

For borehole-response tests, the time-scales over which observations take place are usually much smaller than the normal time-scales over which hydraulic head in a subglacial aquifer varies. In the case of Trapridge Glacier, subglacial water pressure typically varies over a period of hours, whereas the duration of a borehole-response test is, at most, a few minutes (Fig. 2). Furthermore, because our disturbances are small compared to natural pressure variations, it is likely that response tests influence only a local part of the glacier bed in the vicinity of the borehole. Thus, response tests represent small brief perturbations. For these reasons, we suggest that the following simplification is reasonable: the region immediately surrounding the borehole can be treated as a horizontal and homogeneous aquifer in which the pressure gradient will be independent of azimuth for the duration of a response test. If we also assume that the aquifer is isotropic, we can restrict our attention to one-dimensional radial flow in the region surrounding the borehole.

Fig. 2. Variation of water pressure beneath Trapridge Glacier during summer 1990 as recorded by two sensors, 90 P01 (solid line) and 90 P02 (dashed line), 18.8 m apart at the glacier bed. Day 200 corresponds to 19 July. Data were collected by these two sensors at 2 min intervals for the entire period shown. Arrow tips delimit the duration of 16 individual response tests in the same borehole. Natural pressure fluctuations typically occur over time-scales that are large compared to the duration of a single response test. Note that sensors 90 P01 and 90 P02—separated from the hole in which these tests were being performed by, respectively, 19.0 m and 10.4 m—recorded nearly simultaneous responses.

We wish to combine Equation (12) with a constitutive relation that is applicable to a wide range of flow velocities; the small discharges for which Darcy’s law is valid limit its usefulness to linear laminar-flow regimes. Thus, we adopt an expression, suggested by Ergun and Orning (1949, equation 6), that facilitates a smooth transition between laminar and turbulent flow in a porous medium. After conversion to our notation, and modification to allow for radial flow direction, the constitutive relation can be written as

(13)

where q is the radial component of the volume-flux vector, S₀ is the surface-to-volume ratio of solids, and A and Β are positive constants that control partitioning between the two righthand-side terms. With appropriate choices of A and B, the head change will be dominated by the first term on the right side of Equation (13) in laminar flow; the second term will dominate when the flow is turbulent. By inspection, we see that both terms on the right side of Equation (13) are negative when ∂h_B /∂r < 0 and q > 0, corresponding to flow away from the borehole. Also, when ∂h_B /∂r > 0 and q < 0, both terms on the right side of Equation (13) are positive; in this case, flow is directed back towards the borehole.

Equation (13) can be simplified if we assume that the Kozeny-Carman relation is applicable in a subglacial environment. This relationship between permeability k and porosity is given by

(14)

(Carman, 1956). If we combine Equation (14) with the usual definition of hydraulic conductivity Κ where

(15)

(e.g. Bear, 1972, p. 109), then the constitutive relation can be rewritten in terms of hydraulic conductivity. In this case, Equation (13) becomes

(16)

As stated, it is readily apparent that the last expression reduces to the one-dimensional form of Darcy’s law in cylindrical coordinates, q = –K(∂h_B /∂r), if A = 1 and Β = 0. Therefore, without loss of generality, we will set A = 1 for the remainder of our development. With A fixed, a suitable choice for Β still permits relative proportioning between head loss in laminar and turbulent-flow regimes.

Before obtaining the final flow equation, we must solve the quadratic expression in Equation (16) for q. The two roots of Equation (16) are

(17a)

and

(17b)

(Press and others, 1986, p. 145), where we have simplified our notation by defining the constant C₁ = BS₀(1 − n)/8gn³. Because we require that a negative head gradient produce a positive volume flux, the second root must be chosen. Noting that the magnitude of the head gradient can be written as

(18)

simplification of Equation (17b) leads to

(19)

where C₂ = 4K²C₁.

We are now ready to combine the fluid-flow Equation (12) with the constitutive relation in Equation (19). Since we are considering only radial flow, Equation (12) can be written in terms of the radial component of q_j in cylindrical coordinates as

(20)

where we have introduced specific storage S ₈ = ρ _w g(α + nβ). Physically, S₈ represents the volume of water released from storage when a unit decline in hydraulic head occurs in a unit volume of aquifer. For confined aquifers of constant thickness b, it is customary to define aquifer storativity S and transmissivity Τ as follows:

(21)

and

(22)

(e.g. Bear, 1972, p. 215). Using these definitions, and substituting Equation (19) into Equation (20), we obtain

(23)

With Equation (18), the last equation can be expanded and re-arranged to give the partial differential equation governing water flow in the aquifer:

(24)

As a consistency check, we see that a laminar-flow regime corresponds to a weighting coefficient value of Β = 0; in this case, C₂ = 0 and Equation (24) reduces to the standard radial diffusion equation

(25)

(e.g. Marsily, 1986, p. 162), that arises from Darcy’s law.

Before considering borehole-aquifer coupling, we note that values for the weighting coefficient Β can be calculated based on porosity and a critical Reynolds number Re′ for the aquifer. For porous media, a flow transition occurs when 10 ≤ Re ≤ 100 (Marsily, 1986, p. 74); during this transition, the flow regime changes from laminar to turbulent. Ergun and Orning (1949) showed that for Re ≈ 60, with n = 0.35, the two terms on the righthand side in Equation (13) have nearly equal effects upon pressure drop. If we choose Re′ to be the critical Reynolds number at which the terms are of equal magnitude, and simplify our picture of the porous medium by assuming spherical solid grains, we obtain the following expression for the weighting coefficient: Β = 240(1 − n)/Re′.

Coupling Between Borehole and Aquifer Water Motion

The rate at which water flows into or out of the borehole must be equal to the rate at which it leaves or enters the aquifer, if there is no storage within the filter. This is consistent with our previous assumptions; namely, that water compressibility is a negligible component of the rate of change of water-column height and that discharge is uniform across the filter. Under these conditions, continuity of water volume requires that

(26)

where q(r_f, t) is the radial volume flux across the filter. Solving the last expression for the discharge across the filter and combining the result with Equation (16) leads to

(27)

where we have used the fact that ∂h_B /∂r and dh/dt have the same algebraic sign. In effect, Equation (27) couples water flow in the borehole and water flow in the basal aquifer: terms involving borehole-water velocity dh/dt are evaluated from the solution of Equation (11); the radial head gradient at the filter ∂h_B(r _f, t)/∂r is obtained from the solution of Equation (24).

Non-Dimensionalization

In the previous section, we developed a model of water flow in a combined borehole-subglacial aquifer system. The model requires input of material constants, as well as geometric and hydraulic parameters. Unfortunately, model inputs tend to be combined in ways that do not permit straightforward assessment of their individual contributions. (Transmissivity Τ = Kb, for instance, is an important term containing two parameters of interest: hydraulic conductivity and aquifer thickness.) To obtain insight into inherently non-unique parts of the model, we turn to dimensional analysis. This approach provides an efficient way to examine parameter sensitivities and highlights key physical aspects of the model.

Dimensionless Formulation

We start by defining dimensionless variables as follows: time t* = t/t₀; volume flux q* = q/q₀; radial distance r* = r/r₀; water-column height h* = h/h₀; hydraulic head surface forcing The characteristic constants t₀, q₀, r₀ and h₀ are arbitrary, but reasonable choices should involve time, flux and length scales that are representative of the actual physical system. With this consideration in mind, we set the characteristic water-column height equal to the height of water in an open, undisturbed and equilibrated borehole. (As will be subsequently discussed, this value of h₀ also corresponds to that used in dimensional simulations.) With h₀ fixed, the remaining characteristic values follow naturally:

(28)

(29)

(30)

For a porous medium of hydraulic conductivity K, q₀ represents the specific discharge under a constant-head gradient h₀/r₀ As defined, the time constant t₀ approximates the theoretical period of oscillation expected for a vertical U-tube manometer (Prandtl, 1952, p. 50).

The final dimensionless quantities we define are the laminar-flow skin-friction parameter

(31)

the diffusivity parameter

(32)

the transmissivity parameter

(33)

and the Ergun parameter

(34)

The relative importance of skin friction in controlling the rate of laminar water flow in the borehole is characterized by the number ζ. A frictionless borehole wall corresponds to ζ = 0. The number X characterizes the relative importance of diffusion in regulating flow in the subglacial aquifer. As X → ∞, disturbances in the borehole tend to alter the hydraulic head in the aquifer immediately. Conversely, X ≈ 0 means that the head distribution in the surrounding aquifer will not be affected by borehole disturbances. Similarly, the number r indicates the importance of advection in the vicinity of the borehole. In this case, as r → ∞, water tends to be transmitted instantaneously from the borehole into the aquifer, while r → 0 leads to a perfectly unconnected borehole (i.e. one for which there is absolutely no leakage into the glacier bed). Finally, the Ergun parameter ξ characterizes the importance of turbulent transport in subglacial water flow. For ξ = 0, the flow is purely laminar. On the other hand, large values of ξ indicate that the dominant flow regime is turbulent and the deviation from Darcy’s law becomes apparent.

With these definitions, the governing Equations (11), (19), (20) and (26) can be written in dimensionless form:

(35)

(36)

(37)

(38)

Motion of water in the borehole is now described by Equation (35), while Equation (37) governs water flow in the subglacial aquifer. The constitutive relation is expressed by Equation (36) and the input boundary condition, coupling Equations (35) and (37), is given by Equation (38).

Recasting the problem in nondimensional form has simplified the mathematical description; we are left with four dimensionless parameters upon which model solutions depend. These parameters highlight the fundamental physics in the model and allow examination of model sensitivities in a straightforward way.

Simulations

Solution Procedure

A set of finite-difference expressions corresponding to Equation (23) or (37) was solved simultaneously with a system of two first-order differential equations equivalent to the second-order Equations (11) or (35) using an implicit, fifth-order, Runge-Kutta scheme with adaptive time-stepping (Hairer and Wanner, 1991). Computational efficiency was enhanced by using two staggered spatial grids upon which nodes were logarithmically spaced according to the transformation R = ln(r/r′) (e.g. Jarvis and Clarke, 1974), where r′ = 1.0m is a non-dimensionalizing constant that does not rescale the problem. Staggered grids allow specification of hydraulic head on one grid and convenient calculation of head gradient and volume flux on the other. This procedure eliminates the necessity of evaluating derivatives higher than first-order and eases implementation of the boundary conditions. The logarithmic transformation decreases nodal spacing in the vicinity of the borehole, where head changes are the greatest, and reduces the required number of nodes. The resulting algorithm is sensitive to spatial step size. Through detailed analysis, we have determined that consistent solution results are obtained for constant logarithmic step sizes ΔR ≤ 0.22.

Boundary and Initial Conditions

Using staggered spatial grids, the input boundary condition is conveniently given by Equation (26) or (38), with r_f corresponding to the first node on the flux grid. We allow two possibilities for the outer-boundary condition: a boundary of prescribed head or a boundary of prescribed head gradient. In particular, we have used special cases of these general boundary conditions. We consider the subglacial flow layer to be “open” to water flow if, at some distance from the borehole, the preexisting head value h₀ remains undisturbed. This situation imposes a constant-head boundary condition h_Β(r_max,t) = h₀ and, for suitably large values of approximates an infinite aquifer. Alternatively, we consider the system to be “closed” to water flow if, at some distance from the borehole, the hydraulic-head gradient is zero. This condition describes a no-flux boundary, in which case q(r_max,t) = 0. For dimensionless simulations, the outer-boundary conditions are or corresponding to constant head and zero gradient boundaries respectively. Regardless of which boundary condition is used, it is reasonable to expect that, at a sufficiently large radius, disturbances in the borehole will not be sensed — physically or numerically. In testing this limit numerically, our simulations have shown that, for slug and packer tests, the model is insensitive to the prescribed outer-boundary condition when r_max ≥ 20 m.

A “closed” system, such as the one we have described, does not represent a typical aquifer because it fails to transmit a significant amount of water. Nevertheless, we have included this possibility because there are times when large parts of the bed beneath Trapridge Glacier appear not to drain. At these times, subglacial water is ponded but water can still be moved about within these regions. (For instance, water pumped down a borehole sometimes causes flow out of neighbouring boreholes.) Response tests performed under these conditions still allow estimation of subglacial hydraulic properties.

In the case of a connection-drainage simulation, the initial conditions are as follows: the head is uniform within the aquifer so that all nodal positions represent the same hydraulic potential h_Β(r,0) = h₀ ; the borehole is full of water so that the height of the water column is equal to the ice thickness h(0) = h_i ; water in the borehole is stationary dh/dt = 0. For dimensionless simulations, the corresponding initial conditions are , and For slug and packer tests, we again assume that the hydraulic head is initially uniform throughout the aquifer and that the borehole-water level is stationary when the tests begin. Before the system is disturbed, water level in the open borehole represents an equilibrium with the basal aquifer. In the case of a slug-test simulation, the initial water-column height is set equal to the uniform head within the aquifer, minus the height of water h_s displaced by the slug: h(0) = h_B(r, 0) − h_s. In the case of a packer-test simulation, the initial water-column height is simply h(0) = h_B(r, 0). For dimensionless slug-test simulations, the initial conditions are and dh*/dt* = 0. For dimensionless packer-test simulations, the initial conditions are and dh*/dt* = 0.

Sensitivity Analysis

Connection-drainages involve the sudden opening of water-filled boreholes to the basal aquifer. Such distur-bances can result in a wide range of responses, making them particularly well suited for a sensitivity analysis. The effects of the skin-friction parameter ζ, diffusivity parameter X, transmissivity parameter r and Ergun parameter ξ on the character of connection-drainage disturbances are shown in Figure 3. By independently varying these parameters, we can assess their individual contributions to the overall response and also predict what the responses would be for flow systems with vastly different hydraulic characteristics.

Fig. 3. Results of dimensional analysis showing model sensitivity to the four free parameters, a. Wide variations of the skin-friction parameter ζ have virtually no influence on simulation results. Trapridge Glacier boreholes correspond to ζ ≈0.01, indicating that frictional energy loss at the borehole wall is an insignificant process in the coupled borehole-subglacial aquifer system, b. Large values of the diffusivity parameter X suggest that diffusion is an important process in the basal flow layer. However, solution character is insendtive to X as evidenced by the small changes that result from variation of X over several orders of magnitude, c. Small variations of the transmisdvity parameter r have a strong influence on the rate at which water flows out of the borehole and into the aquifer, and vice versa. Despite this sensitivity, r is small in magnitude compared with X. Such a comparison indicates that advection in the vicinity of the borehole is less important than diffusion for the overall system, d. The Ergun parameter ξ characterizes the importance of turbulent transport in the basal flow layer. For ξ = 0 the flow regime is purely laminar. Larger values of ξ correspond to increasing deviations from Darcy’s law. Typical values of the Ergun parameter indicate that flow in the basal layer is strongly regulated by the effects of turbulence. Furthermore, the character of simulated solutions is also sensitive to ξ. For the simulations shown, a transition between underdamped and overdamped responses occurs somewhere in the range 10² < ξ < 10³.

Figure 3a shows the importance of skin friction in controlling the rate of change of water-column height, assuming flow in the borehole is laminar. If h₀ = 50 m, then ζ = 1 requires a borehole radius of less than 0.006 m, an unrealistically small value. Boreholes having reasonable radii, say 0.05 m, lead to smaller values of the skin-friction parameter, corresponding to diminishing importance of skin friction. As can be seen in Figure 3a, for ζ ≤ 0.1 the solution results are nearly identical. This result supports the previously discussed conclusion of Kamp (1976); namely, that under most conditions skin friction can be neglected.

The influence of diffusion on the character of connection-drainage disturbances is illustrated in Figure 3b. Larger parameter values, corresponding to increasing importance of water storage, result in slower drainage rates. It is readily apparent that a large change (five orders of magnitude) in the difiusivity-parameter value results in only small changes in the drainage response; for 10⁴ ≤ X ≤ 10⁷ the differences between solution results are barely perceptible. Thus, connection-drainage solutions appear to be relatively insensitive to diffusion. However, the large values of X calculated from model inputs indicate that diffusion is still an important process governing the flow.

Figure 3c shows the influence of transmissivity in the immediate region surrounding the borehole. The drainage response following a connection is seen to be extremely sensitive to transmissivity-parameter values; small variations of r (less than one order of magnitude) give rise to a wide range of responses. Parameter values in the range r ≤ 8 produce overdamped responses, while for r ≥ 10 the solutions are underdamped. Although solution character appears to be strongly dependent on advection in the vicinity of the borehole, typical values of the transmissivity parameter are small in comparison with the diffusivity and Ergun parameters. This suggests that advection near the borehole plays an overall less important role than either diffusion or turbulence in the basal aquifer.

Figure 3d demonstrates the importance of turbulent transport in the aquifer. For ξ = 0, turbulence is neglected and the flow is purely laminar. Ergun-parameter values 0 ≤ ξ ≤ 100 are seen to result in underdamped solutions, while for ξ ≥ 1000 the system is overdamped. Thus, the character of the response is also sensitive to ξ. Values of ξ in the range 10³ ≤ ξ ≤ 10⁵, calculated from model inputs, indicate that turbulent flow in the aquifer is important in regulating borehole drainage after a connection.

The underdamped responses that are predicted for certain dimensionless parameter values indicate that the borehole-aquifer system is capable of sustaining force-free oscillation. Such oscillations resemble those of the classic spring-mass system: the column of water in the borehole plus some part of the water in the aquifer constitute the mass of the system; the restoring force is provided by the difference between the pressure head in the aquifer and the non-equilibrium water level in the borehole; the damping force arises from the friction that accompanies water flow through the borehole and the aquifer. Transition between overdamped and underdamped responses depends on the mass of water in motion and the system’s ability to transmit this water. Hence, the degree of damping depends on h₀, r_w, S and Τ — geometric and hydraulic quantities that are contained in the dimensionless variables X and r. In the case of a connection-drainage, the inertial force is significant due to the large mass of water involved. The Ergun parameter ξ plays an important role in this case, because the frictional losses in turbulent flow are an important component of the damping force. As Figure 3 illustrates, an underdamped response is predicted when turbulent effects are ignored.

Results

As a demonstration of the theoretical model we compare connection-drainage, slug-test and packer-test simulations with field observations in Figures 4, 5 and 6, respectively. In Table 1 we have listed the model parameter values that were used to obtain the simulated solutions. These solutions were achieved through trial-and-error forward modelling; parameters were adjusted to obtain the best visual fits between the data and the simulated solutions. In this paper, we compare simulated solutions with field observations to demonstrate model usage; thus, the model inputs in Table 1 should not be construed as final results.

Fig. 4. Borehole-connection observation 90CD38 (solid line) from summer 1990 and modelled drainage response (dashed line). Divergence between observed and modelled results becomes apparent near the end of the connection when the head gradient is small and might be due to water input from the drill, a contribution that is ignored in the model.

Fig. 5. Slug test 90ST38-A (solid line) from summer 1990 and simulated slug-test response (dashed line). While the model predicts oscillations about the predisturbed water level, observed oscillations take place about a lower level. Note that this particular test was performed just before noon on day 206, corresponding to 25 July 1990, during a rapidly rising pressure limb (see Fig. 2).

Fig. 6. Packer tests (a) 89PT68-B and (b) 89PT68-C (solid lines) from summer 1989, together with simulated packer-test results (dashed lines). As with slug tests, the model predicts oscillations about the predisturbed water level, whereas actual fluctuations are about a somewhat lower level. These two tests were performed in the same hole on two successive days — the borehole was re-opened in the morning of the second day.

Table 1. Parameters for borehole-response simulations*

A connection-drainage observation in Borehole 38 during summer 1990 (designated 90CD38) is shown in Figure 4 along with a modelled drainage response. During this observation the borehole water level dropped approximately 20 m in 20 s; thus, the mean water velocity was about 1.0 ms⁻¹ during the initial moments of connection. This velocity, together with the simulation input values, implies a mean specific discharge of roughly 0.3 ms⁻¹ at a distance of 0.1m from the borehole center. Such high water-transfer rates motivated our use of a non-linear constitutive relation to characterize the subglacial aquifer. As the borehole water level approaches the equilibrium head value, simulated and observed results diverge. This model predicts a more rapid return to predisturbed conditions than actually occurs. For this particular connection, modelled and observed results converge to the background head value approximately 120 s after the connection.

Figure 5 shows the first of a series of slug tests performed in Borehole 38 during summer 1990 (designated 90ST38-A) and a simulated slug-test response. The negative initial water-column displacement corresponds to removal of the slug. After the slug is withdrawn, this model predicts an oscillatory response about the predisturbed water level. The measured response is also oscillatory — however, it does not oscillate about the predisturbed water level but rather about a lower level. Similar responses have been observed during several different series of slug tests; repeated tests result in lower re-equilibrated water levels. We shall discuss this behaviour subsequently.

Two packer tests performed in Borehole 68 during summer 1989, together with simulated responses, are shown in Figure 6. The first, 89PT68-B (Fig. 6a), was performed around 1820 h on day 209. The second, 89PT68-C (Fig. 6b), was performed around 1055 h the following day. Because the borehole had frozen shut overnight, it was re-opened in the morning of the second day. As with slug tests, this model predicts oscillations about the predisturbed water level, whereas actual fluctuations are about a somewhat lower level. The oscillations recorded during packer test 89PT68-B are more rapidly damped and are of longer period relative to those observed during test 89PT68-C. Comparison of model input parameters (Table 1) reveals that smaller head and larger conductivity values result in higher-frequency oscillations that are less quickly damped.

Discussion

Generally good agreement between observed and modelled results demonstrates that borehole-response tests, together with the theoretical framework we have developed, can be used to estimate subglacial hydraulic properties. As previously mentioned, the simulations shown in Figures 4, 5 and 6 are intended to demonstrate model usage. In these simulations, parameter variations were intentionally limited to highlight those to which we have found the model is most sensitive. It is important to realize that focusing on sensitive parameters is a natural tendency with forward modelling, and that this approach can fail to produce a fully consistent set of hydraulic parameters. For example, differences in the modelled responses shown in Figure 6 are due only to changes in initial head and hydraulic conductivity values. For these simulations, all other parameters were held constant. However, according to Equations (14) and (15), changes in hydraulic conductivity must be accompanied by changes in porosity and/or the surface-to-volume ratio of solids. Such inconsistencies must be resolved when seeking actual hydraulic parameter estimates, again emphasizing that the model inputs listed in Table 1 should not be construed as final results.

Significant temporal changes in hydraulic conductivity, suggested by model inputs for the simulations shown in Figure 6, are not expected in conventional groundwater flow systems. However, such changes are not unreasonable in a subglacial sediment layer. One possibility is that water discharge from the drill causes fine sediments to be entrained and transported away from the borehole. Alternatively, basal sediments might be deforming at a non-uniform rate, thereby changing the porosity. In either case, hydraulic conductivity would be altered.

Discrepancies between Observed and Modelled Results

In comparing measured and simulated responses, several discrepancies have been noted, suggesting the possibility that there are additional properties or processes that are not included in our model. Based on initial simulations, we have discovered the following discrepancies between observed and modelled results: (i) simulated borehole connections recover to the equilibrium head value more rapidly than observed connections; (ii) slug and packer tests result in lower re-equilibrated water levels, whereas the models that we have used predict that pre-disturbed levels will be regained; (iii) initial head values specified for slug- or packer-test simulations are typically lower than those that actually existed at the times of the response tests.

Subsequent inverse modelling, using the non-dimensional formulation, has shown that these discrepancies can be eliminated with judicious choices of input parameters; modifications to the model are not required. In a separate paper, we will present better estimates of the subglacial hydraulic properties of Trapridge Glacier and we will demonstrate the importance of using inversion as part of the overall procedure.

An interesting aspect of slug and packer tests is that re-equilibrated water-column heights are consistently observed to be lower than the predisturbed levels. This situation is schematically illustrated in Figure 7. We have considered a number of possible explanations for this behaviour: changes in background pressure, changes in sensor position, poor testing procedure, temporary water storage in the aquifer and water expulsion from the immediate flow region. Basal water storage and/or expulsion from the flow layer are the most likely causes of this effect (Stone, 1993).

Fig. 7. Conceptual slug test showing reduced water levels following response tests. Part a shows slug insertion and removal, along with the associated displacements of borehole-water level. Temporal changes in water-column height h, as measured by a pressure sensor in the borehole, are sketched in part b. Dotted lines indicate initial water levels, and the dash lengths in part b represent approximately 1 min of observation. The slug is inserted into the borehole at point A and the disturbance produces an immediate rise in water level, corresponding to the A–B transition. The insertion creates a pressure imbalance at the bottom of the borehole that causes water to flow from the borehole into the subglacial aquifer, as shown between Β and C. At C, the borehole-water level has returned to its original position and the slug is removed, again creating a pressure imbalance. In response to this imbalance, water flows back into the borehole, as indicated between D and E. The water level does not, however, recover to its initial value within the observation time window.

Generalizations to Other Flow Systems

We have previously stated that our model can be applied to distributed sheet and channelized flow regimes as well as to ground-water flow. For drainage through a very thin layer, the representation provided by our model is directly analogous to sheet-like flow as discussed, for example, by Weertman (1972). To illustrate this connection, consider the relations governing flow between infinite parallel walls separated by a distance b. For one-dimensional laminar flow, the analytic expression relating specific discharge and hydraulic head gradient is

(39)

(Todd, 1959, p. 315; Bear, 1972, p. 692), and for turbulent flow, standard empirical formulas give

(40)

(Francis, 1975, p. 218–19; Henderson, 1966, p. 91–101) where f₀, n′ and C are respectively the Darcy friction factor, the Manning roughness parameter and the Chézy coefficient. Note that in presenting the formulas in (40) we have used the hydraulic radius for an infinite sheet of thickness b and also made the sign corrections that are necessary to account for flow direction. The analogy between flow through porous media and Weertman-like sheet flow is immediately obvious if we express Equation (16) in terms of x derivatives and set A = 1:

(41)

For low values of specific discharge, Equation (41) gives q = –K(∂h_B /∂x) and suggests the correspondence Κ = ρgb²/12n between Darcian flow and laminar sheet flow. For large values of specific discharge, Equation (41), in accordance with Equation (18), gives

(42)

and suggests the correspondence

(43)

Thus, for both laminar and turbulent sheet flow, relationships between specific discharge and hydraulic head gradient are represented by analogous expressions in our model. Furthermore, it is worth noting that the constitutive relation in Equation (41) does more than just represent laminar and turbulent flows — it also provides a continuous solution for the transition between these regimes. Analyses of linked cavity configurations (Walder, 1986; Kamb, 1987) yield expressions that, like Equations (40), involve but the coefficients involve numerous geometrical variables that would be challenging to disentangle. For cavity configurations, networks of arborescent or braided channels, and other distributed systems, our model provides estimates of the hydraulic properties of permeable layers that “effectively” characterize the actual flow systems. In these cases our model suggests alternate — but hydraulically equivalent— representations of the real drainage configurations.

Concluding Remarks

We have presented a physical framework that allows estimation of subglacial hydraulic properties, when combined with field observations of boreholes responding to induced changes in basal water pressure. In general, agreement between responses predicted by the model and those we have observed is satisfactory. To this extent, our theory seems to provide a reasonable description of the coupled borehole-subglacial flow-layer system. We have also shown that our model can be applied to a variety of distributed drainage systems. As such, it is potentially useful for many wet-based glaciers.

Because response tests are restricted to a small region of influence near the borehole, application of this model to many different borehole response tests — separated both in space and time — is a means by which the distribution of hydraulic parameters can be quantified. This approach is useful for understanding the heterogeneous properties of subglacial drainage systems.