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. Reference Cooper, Bredehoeft, Papadopulos and BennettCooper and others, 1965; Reference KambKamp, 1976; Reference KippKipp, 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:

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 (Reference MarsilyMarsily, 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

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

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

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

where sgn (x) is the algebraic sign function (Reference BracewellBracewell, 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

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

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

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,

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 (Reference PrandtlPrandtl, 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

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 (Reference Kay and NeddermanKay 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

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 Reference Cooper, Bredehoeft, Papadopulos and BennettCooper and others (1965) and by Reference KippKipp (1985). Like Equation (11), these expressions were based on conservation of momentum. Reference Kabala, Pinder and MillyKabala and others (1985) used the framework of Reference Cooper, Bredehoeft, Papadopulos and BennettCooper and others (1965) in numerical modelling of the responses of well-aquifer systems to sudden changes of water levels. In their developments, Reference Cooper, Bredehoeft, Papadopulos and BennettCooper and others (1965) and Reference KippKipp (1985) neglected loss of momentum due to skin fricüon. It was shown by Reference KambKamp (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:

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 (Reference JacobJacob, 1940; Reference CooperCooper, 1966; Reference WiestWiest, 1966; Reference Gambolati and FreezeGambolati and Freeze, 1973; Reference Bear and VerruijtBear 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. Reference Verruijt and WiestVerruijt (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 Reference Gambolati and FreezeGambolati 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 Reference Ergun and OrningErgun 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

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

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

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

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

and

(Reference Press, Flannery, Teukolsky and VetterlingPress 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

simplification of Equation (17b) leads to

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

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:

and

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

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

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

(e.g. Reference MarsilyMarsily, 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 (Reference MarsilyMarsily, 1986, p. 74); during this transition, the flow regime changes from laminar to turbulent. Reference Ergun and OrningErgun 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

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

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

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:

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 (Reference PrandtlPrandtl, 1952, p. 50).

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

the diffusivity parameter

the transmissivity parameter

and the Ergun parameter

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:

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.

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 (Reference Hairer and WannerHairer 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. Reference Jarvis and ClarkeJarvis 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 Reference KambKamp (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.

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 (Reference StoneStone, 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 Reference WeertmanWeertman (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

(Reference ToddTodd, 1959, p. 315; Reference BearBear, 1972, p. 692), and for turbulent flow, standard empirical formulas give

(Reference FrancisFrancis, 1975, p. 218–19; Reference HendersonHenderson, 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:

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

and suggests the correspondence

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 (Reference WalderWalder, 1986; Reference KambKamb, 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.