7.1. Approach

Before attempting a fuller explanation of the data, we develop a theory of the operation of the hydraulic system in which the flux of meltwater and sediment transmissivity determine the water pressure in the sediments for a given ice load. There are strong diurnal variations in surface meltwater production, which are transmitted through the glacier and till into the aquifer and produce variations in water pressure and sediment effective pressure. These variations in basal water pressure and the average water-pressure profile depend on the rheology of the subglacial sediments under the varying ice load caused by the advance and subsequent retreat of the glacier. Although it is clear that the tills at the site underwent shear deformation during the mini-surge, we believe that our experimental data only record pressures in the sediments beneath the deforming horizon, because of the fall in water pressure that occurs immediately before transducers cease recording, which we interpret as reflecting the beginning of sediment dilation produced by shear deformation. We therefore ignore shear deformation in the theory and modelling presented below.

We use poroelasticity theory to account for the spatial and time-dependent distribution of basal water pressure during the period of rapid ice loading as the glacier advances. The objective is to use as many of the available field data as possible as input to a model and to restrict the intrinsic model parameters to the minimum possible number, producing a model capable of predicting water-pressure profiles in the till and aquifer in response to an arbitrary time-varying recharge in which consolidation of the subglacial till layer is taken into account. Physical parameters were restricted to the permeability and volume compressibility of subglacial sediments (responsible for the compaction/dilation of porous material) and an additional parameter related to the glacier’s ability to deliver water from the ice surface to the subglacial till. This latter was introduced for the general case, but proved to be redundant for our experimental site because of the rapid transmission of surface water to the bed.

The irregularity of surface elevation, till thickness, till permeability and, most importantly, the evolving upper boundary condition of a spatially distributed ice load moving upslope, rendered a finite-difference or finite-element approach inappropriate. The chaotic nature of the surface water production data (Fig. 9), containing frequency components from 6 hours to almost half a year, presents a further complication. To deal with these complexities we have developed a mathematical model similar in operation to an analogue computer, where the entire area of interest is represented by a one-dimensional array of interconnected vertical cells (Fig. 12). Each cell represents a hydraulically interconnected thin vertical slice of ice, till and aquifer. The advantage of this approach is that the hydrodynamics of water in the cell can be analytically resolved for a given water production rate at the ice surface (at the top of the cell) with the incoming and outgoing water flux to the left and right of the cell accounting for the horizontal water dynamics. The total water-pressure field can be computed for given Fourier components of the runoff data and ice load through the numerical solution of a complex matrix equation using a standard Fortran library subroutine. The final result is given as a sum of the inverse Fourier components. The model was tested against analytic solutions for some simple approximations and produced extremely accurate results for spatial resolution, even for a relatively small number of vertical cell elements. The results are consistent with the assumption that subglacial groundwater flow lies in a flow-parallel vertical plane, apart from at the end of the experiment, when we suggest that local transverse channels developed some distance back from the maximum extent of the advance.

Fig. 12. The cell model used for hydraulic analysis. The assumed longitudinal flow system is represented by a one-dimensional array of cells representing ice (i superscript), till (t) and aquifer (a). and Pa are water-pressure potentials at the base of the glacier, in the middle of the till and at the top of the aquifer, respectively, with a varying thickness, , of each column. The water exchange between different cells is assumed to take place in the aquifer only. Apart from the first and last, cells have a horizontal extent of 1m along the transect.

7.2. Water flow in consolidating subglacial sediments

7.2.1. Introduction

In the theory of consolidation, soils react differently to the same effective stress increment depending on their deformation history. Figure 13 shows a typical stress/strain relationship for a soil under uniaxial loading. The virgin consolidation curve, representing the initial compression of the soil during sedimentation and subsequent loading, shows a steady irreversible consolidation along the plastic yield envelope (OAB), along a ‘normal consolidation line’ (NCL) with a fairly high coefficient of compressibility. If the soil is unloaded then reloaded, it follows an expansion (swelling)–recompression loop (BCD). In this region it is ‘overconsolidated’. A soil can store its history of prior, higher consolidation (‘pre-consolidation’), with a lower void ratio than it would have had it undergone virgin consolidation up to its current load. Reference Boulton and DobbieBoulton and Dobbie (1993) show that vertical profiles of pre-consolidation can be preserved in fine-grained sediments with gradients significantly greater than normal gravitational gradients. They conclude that these were a memory of very high water potential gradients in the till created by strong subglacial rates of water flow through it.

Fig. 13. Consolidation states for soils/sediments. Void ratio (e) = 1/(1 – n) where n is porosity. NCL is normal consolidation line, CSL is critical state line, p is the effective pressure and p₀ the initial effective pressure. BCD illustrates the pre-consolidation process (BC is expansion, CD recompression).

Although it is clear that sediments (soils) react nonlinearly to varying applied (effective) stresses, we can nevertheless reasonably apply poroelastic theory to describe the hydraulic response of subglacial sediments using the following observations of specific ice loading and runoff history. The varying ice load and groundwater-pressure data shown in Figure 10 can be decomposed into three principal components:

a slowly increasing spatially distributed glacier load followed by a stationary state,

a slowly increasing average water pressure followed by a slow decline,

a fast (diurnal or near-diurnal) fluctuation in water pressure.

We assume that the difference between the glacier overload, σ₁(t), and the averaged groundwater-pressure component, , is effectively responsible for the irreversible consolidation of till and aquifer sediments (the OAB part of strain/stress curve along NCL in Fig. 13) with some coefficient of volume compressibility, M_V , whilst the fast-fluctuating pressure component, , causes a reversible response in the till and aquifer sediments, characterized by the coefficient of volume compressibility, m_V.

Our decomposition of the pressure field into fast and slow components relies heavily on the observed groundwater-pressure data, where fast diurnal fluctuations are clearly superimposed on an otherwise slowly changing pressure profile. Strictly speaking, this decomposition is valid only when the two following conditions are satisfied:

(1)

These conditions are not restrictive during a major part of seasonal groundwater evolution, with several obvious exceptions:

during the first few days after glacier movement over a site (e.g. days 7–10 at 125m along the transect; see Fig. 10),

during some extremely strong peaks of groundwater pressure (e.g. day 140 at 85m along the transect; see Fig. 10),

during the fastest stage of glacial advance when all records show that water pressures increase at a rate that is faster than the rate of increase in ice load (e.g. between days 90 and 105 at 85m along the transect; see Fig. 10).

This decrease of effective pressure will cause expansion (swelling) of sediments along curve BC in Figure 13, and can bring porous materials to the critical state (CSL line) and permit shearing of till. This shearing can substantially increase the local transmissivity of the till, and when it occurs through the whole thickness of till, can make a very fast hydraulic connection between the glacier sole and the aquifer. Here, we address only the hydraulic response of subglacial sediments before the onset of active shear deformation in them.

7.2.2. Model description

In Appendix A, we derive a closed equation describing the evolution of water pressure in the compressible sediments under applied (external) load, σ(t), from the first principles of the theory of poroelasticity. In its final form this equation reads:

(2)

where ϕ and k are the porosity and specific permeability of sediments, η is the dynamic viscosity of water, c_W and c_s are the compressibility coefficients of water and solid grains, and

(3)

is the bulk compressibility of dry sediments (with density p_T = p_S(1 — ϕ), where p_s is the density of mineral grains).

For weakly consolidated sediments such as the subglacial tills at Breiðamerkurjökull, the bulk compressibility, c_T (≈10^–8 to 10^–5Pa^–1), is much greater then the compressibility of mineral grains, c_s (≤10^–11 Pa^–1), so we can neglect terms with c_s. In this case Equation (2) reduces to the Terzaghi consolidation equation with an additional correction caused by water compressibility:

(4)

The coefficient of hydraulic diffusivity (or Terzaghi’s coefficient of consolidation), c_V, in front of the Laplacian is, in our notation, given by:

(5)

where K = p_Wgk/p is the hydraulic conductivity and g is the acceleration due to gravity. We can identify our (poroelastic) bulk compressibility coefficient, c_T, with the coefficient of volume compressibility, m_v ≡ c_T, traditionally used in soil mechanics.

While we can reliably disregard the contribution of solid grain compressibility on water-pressure dynamics in Equation (4), the effect of finite water compressibility can be important and, in the case of high gas saturation, may be dominant. The compressibility of water with volume gas content, v_G (volume of gas per unit volume of gas-saturated water), can be evaluated as:

(6)

With the coefficient of water compressibility c_W ≈ 4.4 × 10^–10Pa^–1 and compressibility of gas c_G ≈ 0.76 × 10^–5Pa^–1, the compressibility of saturated water is completely dominated by gas for v_G ≥ 0.6 × 10^–4. As a result, whenever ϕv_Gc_G ≈ 10^–5 ϕv_G becomes of the order of the sediment coefficient of volume compressibility, m_v (~10^–8 to 10^–5Pa^–1), the effect of water compressibility on water-pressure dynamics in Equation (4) becomes essential or dominant (we currently ignore this, but it can result in a serious overestimate of permeability).

Application of a Fourier transform to Equation (4), with respect to time, yields an equation for water pressure in the frequency domain:

(7)

We can formally rewrite Equation (7) in the form:

(8)

introducing a frequency-dependent coefficient of compressibility, c_Tω), where c_Tω) = m_v for high-frequency (diurnal) fluctuations in water pressure, and c_T(ω) = M_v for low frequencies, corresponding to slow changes in stress (overload of the advancing glacier) and to the general build-up of water pressure, following the discussion in the previous section. Formally, this means that we create an approximate constitutive relationship for the sediments satisfying a causal relationship:

(9)

where Tre _ij is the absolute dilatation of sediments and T is the integration variable. Equation (9) assumes a viscoelastic response of the sediments to an arbitrary history of (effective) stress, σ(f), and, in general, is not applicable to the description of deformation, including irreversible consolidation. However, in the case of an advancing glacier, the slow component of effective stress shows a general increase with time, following the ‘consolidating’ direction along line NCL in Figure 13. Provided that we do not consider glacier retreat, the rheology given by the constitutive Equation (9) with a frequency-dependent compressibility coefficient, c_T, is a reasonable model for sediment consolidation. When the compressibility of water can be neglected, Equation (8) reduces to:

(10)

with Terzaghi’s law of effective stress clearly assumed for all frequencies.

We stress that we introduce a frequency-dependent compressibility as a technical device to describe essentially non-linear behaviour of subglacial sediments by a linear system of equations, as given by the consolidation curve in Figure 13. The logic that renders this approach physically reasonable was discussed in detail in section 7.2.1. In general there are four, initially independent, material parameters in standard poroelasticity theory (Reference BiotBiot, 1941, 1962) such that any modernization of the constitutive relationship in a way similar to Equation (9) inevitably results in a more complicated rheology than simple poroelasticity with a viscous frame, particularly with respect to the shearing properties of water-saturated sediments.

Thus far we have ignored the gravity-induced hydrostatic component of water pressure. This can be accounted for by adding a constant term, –pwg, to the vertical hydraulic pressure gradient.

7.3. Water flow in the till

Consider a permeable, porous layer of thickness d representing a till sandwiched between the glacier sole and a more permeable aquifer of thickness D. If the specific permeability, K, and the thickness of the till are much less than the corresponding values for the aquifer (kd c k_AD), water flow in the till can reasonably be considered as a onedimensional vertical flow. The one-dimensional solution of differential Equation (10) with specified boundary conditions at the top and bottom of the till is derived in Appendix B in the frequency domain and is given by:

(11)

where z is the vertical (depth) coordinate (0 ≤ z ≤ d), sinh is the hyperbolic sine of a complex variable, and the complex parameter, A, is defined as:

(12)

with (±) corresponding to the positive or negative frequencies in the Fourier transform, and subscript ‘T’ referring to the physical parameters of the till. The solution in Equation (11) gives a distribution of water pressure across a layer, as a function of water-pressure values at the ice-till interface, pi, and at the till-aquifer interface, p ₁. The second term on the righthand side of Equation (11) is responsible for an additional pressure induced by compression of the till layer under applied stress σ.

There are two useful, closely related parameters that define characteristic temporal and spatial scales of the process. The characteristic response time of water pressure in the till to pressure fluctuations on its boundaries is:

(13)

and the characteristic frequency-dependent diffusion length (see Appendix B for details):

(14)

This determines how far into a till a pressure perturbation of given frequency, J, on its boundary can penetrate, and determines two types of behaviour:

1. 1. Undrained loading. For short-term perturbations, when ωt _t ≪ 1, the diffusion distance is less then the till thickness: δ(ω) ≫ d. In this case, the water pressure in the till layer is relatively little affected by the water-pressure variations p ₁ and p _A on its boundaries, and will be restricted to a small thickness ≈δω) near the top and bottom boundaries. In addition, however, increases in applied stress, σ(ω), such as that due to increased loading during glacier advance, will be borne by the interstitial water with induced water pressure p(ω) ≅ σ(ω), which will be unable to drain and thereby to transfer the increased load to interparticle contacts. Such undrained loading can readily cause failure through almost the whole thickness of the till.

2. 2. Drained loading. For longer-term perturbations of pressure, when ωT _T ≪ 1, the diffusion length is greater than till thickness, and the runoff-induced pressure (hydraulic potential) changes almost linearly with depth in the till, resulting in a quasi-stationary (Darcy) flow driven by the pressure difference at the top and bottom of the till layer. Increments of load are readily borne by interparticle contacts because of rapid drainage of the till.

From Equation (13), we see that the characteristic response time depends on the permeability and compressibility of sediments. For a 1m thick silty-sandy till layer with hydraulic conductivity of about 6 × 10^–7ms^–1 and skeleton compressibility (typical for soils) ~2 × 10^–7Pa, the characteristic response time is of the order of 1 hour. Given that the characteristic response time is proportional to the square of till thickness, a 3m layer of the same till will show undrained behaviour for pressure fluctuations with semidiurnal period, but drained behaviour over longer periods. A silty-clay till with a hydraulic conductivity less than 10^–10ms^–1 would have a response time of about 240 days. A glacier advance similar to that monitored by our experiment, in which the rate of loading is relatively high, would ensure that the till remained in an undrained state even without any fluctuations in recharge from the overlying glacier.

The water flux through the till-aquifer interface is given by the derivative of the solution (11) over depth z:

(15)

For very high frequencies ωT _T ≫ 1) the water flux through the till-aquifer interface is insensitive to the conditions at the top of the till (the pressure term p₁ in Equation (15) is suppressed by an exponential factor sinh^–1(λd) ≅ exp . High-frequency pressure variations in runoff therefore cannot penetrate through the till layer (see Appendix B for more details).

At the low-frequency (drained) limit, Equation (15) yields:

(16)

The first term on the righthand side of Equation (16) refers to standard Darcy flow caused by the pressure difference across the layer. The second term gives an extra influx into the aquifer caused by till consolidation or compression/decompression under variable effective stress.

To resolve the pressure distribution in the aquifer we need an expression connecting the flux, Q _TA, across the till-aquifer interface with the influx from the glacial ice into the till layer, Q₁. This flux is given by an equation similar to Equation (15):

(17)

Expressing p₁ in this equation as a function of Q₁, p_A and stress, σ, and substituting into Equation (15), gives:

(18)

This is a remarkable result. if the flux, Q₁, supplied by the glacier to the subglacial hydraulic system is known (in our case we presume that it is equal to the water production rate at the glacier surface) then Equation (18) gives the water flow exchange between the till and aquifer sediments (including consolidation-induced correction) as a function of the pressure distribution in the aquifer only. The complex influence of the till layer is accounted for by the hydraulic conductivity, k_T, and a single dimensionless parameter λd (with λ defined by Equation (12)).

7.4. Hydraulic coupling between ice and till

if we assume that during periods of strong surface melting most discharge from the base of the glacier into the till represents surface meltwater that has penetrated through the glacier, the simplest conceivable situation is one in which the water pressure, p’, at the top of the till is determined by the water head, H, in the glacier:

(19)

The head variation is determined by the difference between the discharge rate from the surface, R, and the water flux, Q₁, into the till:

(20)

where ψ is the water volume per unit volume of glacial ice. We do not specify the type of conduit system, assuming only that it is a well-connected system of conduits and that an additional pressure gradient caused by water flow through the glacial ice is much less than the hydrostatic gradient p_Wg. Substituting for H from Equation (19) into (20) yields:

(21)

or, applying the Fourier transform,

(22)

Excluding pressure p₁ from Equations (22) and (17), we have the flux, Q₁, crossing the ice-till interface as a function of R and pressure, p_A, at the top of the aquifer:

(23)

where

(24)

and another characteristic time parameter, T, is given by ice and till physical parameters:

(25)

Note that even in the case when till permeability is high enough to provide an easy hydraulic connection with the aquifer, i.e. when ωt _t ≪ 1, parameter ωT can be high enough to suppress flow through the ice-till interface. For a till of hydraulic conductivity K ≈ 6 × 10^–7 m s^–1, ωT ≈ 10² ψ for diurnal water-pressure fluctuations. Consequently, where the water content of glacier ice is several per cent, ωT ≥ 1 and the runoff rate in the first term of Equation (23) is multiplied by a small factor, 1/ωT. Another consequence is a phase shift, exp (—iπ/2), in the flux, Q₁, with respect to the runoff rate, R, caused by a large value of ωT in the denominator of Equation (23). The second term on the righthand side of Equation (23) represents an upward water flow component, caused by fast till compaction/consoli-dation or by excessive aquifer pressure. Equation (23) ignores the fact that the local ice thickness will determine an upper limit to the water pressure.

More careful considerations, including adding effects caused by viscous water flow through glacial conduits, show that further corrections to Equation (22) contain non-linear terms involving water-flux and water-pressure products (see section 11).

In conclusion we give an expression for water pressure at the ice-till interface as a function of runoff, R, and water pressure at the till-aquifer interface. Substituting Equation (23) into Equation (17) and resolving it with respect to the ice-till interface pressure, p ₁, we have

(26)

if the water pressure, p_A, at the top of the aquifer is known, Equation (26) gives water pressure at the top of the till layer and, as a consequence, the pressure distribution through the till (given by Equation (11)). Note that, apart from the obvious dimensional parameter, dη/k_T, providing a characteristic value of pressure drop in pascals through the till layer caused by runoff of magnitude R = 1 ms^–1, only two nondimensional parameters, ωT and λd, are involved in the important relationship (26).

7.5. Hydraulic coupling between till and aquifer

Consider the vertical cell through the glacier, till and aquifer, schematically shown in Figure 12. If the aquifer has a hydraulic conductivity and thickness much greater than the till layer above it, there will be almost vertical water flow through the till layer that supplies a generally subhorizontal water flow in the aquifer (cf. Reference Boulton and DobbieBoulton and Dobbie, 1993). Numerical restrictions on physical parameters which justify this scenario are given below. An analysis of water mass balance in the ith aquifer cell of vertical thickness D and horizontal extent L, where index i specifies the discrete spatial coordinate of the aquifer slice along the y axis, yields:

(27)

where c_A is the volume compressibility coefficient of aquifer, analogous to c_T used for the compressibility coefficient of till. The concept developed in section 7.2 is applied to the volume compressibility of the aquifer. The aquifer responds with compressibility c_A = M_VA to a steady (low-frequency) ice loading and general build-up in water pressure, and with a different (smaller) compressibility,

c_A = m_VA, to the high-frequency diurnal fluctuations in water pressure. The lefthand side of Equation (27) describes the water-content change in cell volume caused by variations in effective stress, σ _i – p,-. The last two terms on the righthand side of the equation account for water exchange with the adjacent cells (Fig. 12), i.e. for the horizontal water flow in the aquifer. The easiest way to solve Equation (27) is to replace horizontal water fluxes qi+1, i q _i, i–1 and with approximated water fluxes between the cells:

(28)

where k _A is the specific permeability of the aquifer. Following our previous discussion, this approximation applies only to the drained regime, when the diffusion length (for the highest available frequency) is greater than the cell thickness in the direction of flow. In the case of a till layer of thickness d, the relationship between the diffusion length δ(ω) and d is not known a priori, and high-frequency fluctuations in water discharge can lead to an undrained condition. The case of an aquifer is, however, different. We can always choose a cell thickness, L, small enough to satisfy the drained condition, even for the highest (Nyquist) frequency. In this case, the system of Equation (27) can be written as

(29)

where

(30)

is the hydraulic diffusivity (consolidation coefficient) of the aquifer. The Fourier transform version of Equation (29) yields a linear system of equations for pressure p_i in the aquifer cells which can be readily solved.

This approach demands a selection of small values of L and, consequently, an unnecessarily large number of aquifer cells. A better approach is to apply solution (B8) of Appendix B to the horizontal flow through the slab of aquifer with an arbitrary cell thickness. This solution can be written as:

(31)

This equation gives a distribution of water pressure in the ith cell in the horizontal direction along the dominant flow, when pressure values in the adjacent cells p_{i –1} and p _i+1 are considered as given boundary conditions. The last term on the righthand side of Equation (31) describes the effect of the vertical influx, Q_TA, from the overlying till layer. A possible consolidation of the aquifer sediments is accounted for by the term with applied stress, σ. The pressure value in the centre of the ith cell is equal to p_i ≡ p_i (y = L/2). If the thickness of till is variable and the applied stress (σ = σ(y) = σ_i ) is non-uniformly distributed along the glacier bed, we choose the size of a discrete cell, L, small enough to account for the observed variation in till thickness, glacier bed elevation profile and ice surface profile.

Substitution of the influx Q_TA = Q_TA(R, σ, p), given by expressions (18) and (23), into Equation (29) results in a closed matrix equation for the unknown array of pressure values p_i in the centres of aquifer cells with i = 1, N:

(32)

where f₁ and f₂ are hyperbolic functions defined in Equation (24), and

(33)

The only dimensional coefficient,

(34)

on the righthand side of Equation (32) determines the amplitude of aquifer pressure response in pascals to the runoff influx, R, in ms^–1. Two new dimensionless coefficients appear in Equation (32):

(35)

and

(36)

Equation (32) can be formally rewritten as a matrix equation:

(37)

with the explicitly known matrix coefficients Z_ji and array A_j and given Fourier components of runoff, R, and ice load σ_j . This matrix equation can be resolved with the help of any standard numerical algorithm for complex matrix inversion. Application of an inverse Fourier transform yields the water-pressure distribution in the aquifer for an arbitrary runoff and ice-loading history. The pressure of water in the till layer is given by Equations (26) and (11) if the water pressure in the aquifer is known.

The spatial distribution of till (and/or aquifer) thickness and physical properties such as hydraulic conductivity and compressibility can be accommodated prescribing individual characteristic times to different cells. By recognizing that the contrast in conductivity between the till and underlying aquifer will result in vertical and horizontal flows, respectively (Reference Boulton and DobbieBoulton and Dobbie, 1993), we are able to solve a two-dimensional hydraulic problem with an inhomogeneous distribution of physical parameters and boundaries (including moving boundary conditions; see section 7.6). This robust approximation allows us to dramatically simplify the calculation, so that the vertical dynamics can be completely resolved analytically, and the horizontal spatial distribution is found by complex matrix inversion. The approach is similar in spirit to an analogue computer method, but has been implemented by means of a numerical algorithm using standard subroutines for fast Fourier transform and complex matrix inversion.

7.6. Complications caused by the moving glacier boundary

In Equation (37) we prescribed a spatial index, i, to the source on the righthand side of the matrix equation to describe a spatial distribution of till (and/or aquifer) physical parameters. This also allows us to account for the glacier advance as a moving boundary condition on the top of the till layer. We initially assume that recharge from the glacier surface to the bed is spatially uniform over the area of data collection or, more generally, as far as the water divide. Before they are reached by the advancing glacier, transducers do not show significant pressure variations. This is taken into account by multiplying our runoff dataset by a Heaviside function H(t — t_i ), where t_i is the time of glacier advance to a location specified by spatial index i.

The two last terms on the righthand side of Equation (36) represent the effect of consolidation in the till and aquifer sediments. The terms are functions of the ice load, σ_i , in each cell element and are different for different cells. This requires interpolation of glacier surface elevations before undertaking a Fourier transform of the time-dependent ice load. Both problems are caused by the moving boundary condition on the top of the till layer. It complicates data preparation without affecting our computational algorithm. Nevertheless, the evolving upper boundary conditions essentially prevent us from estimating aquifer hydraulic parameters directly in the Fourier domain.

7.7. The subglacial model

The subglacial model is represented by 127 cells, where all cells have the vertical structure shown in Figure 12 and, with the exception of the first and the last, have the same horizontal extent △ = 1m. The first cell is introduced to account for the proglacial area, which has not been overrun by the glacier advance, though it may be influenced hydrologically, and the last, extended, cell represents the rest of the subglacial system from the location of the 125m transducer to the water divide. The thickness of till in each cell is taken into account. The effect of elevation is accounted for by means of an additional gravity-induced hydraulic gradient. The horizontal thickness of the last cell, representing the whole area from the location of the 125m transducer up to the unknown water divide, is the only free parameter in the model. It can be readily determined by matching the modelled water-pressure spatial profile in the aquifer along the transect with the experimental results.

The pressure histories of transducers are all very similar (Fig. 10). They show five successive generalized pressure regimes:

1. 1. The aquifer transducer lies above the water table until the glacier margin advances to within 5–10m of the site.

2. 2. As the glacier advances to within about 5–10m of the site, the water table rises. Pressures rise first in the aquifer, and then in the till. Figure 11a shows the 65m transducers, where an upward pressure gradient is maintained for about 5 days, and for 3–4 days after glacier overriding of the site, when the site lies about 4m behind the advancing ice margin. This reflects upward expulsion of water from the aquifer towards the proglacial zone.

3. 3. Shortly after the glacier has overridden the site and the establishment of a downward pressure gradient (Fig. 11b), the pressure at the top of the aquifer initially follows the ice pressure at the 12, 30, 65 and 85m sites (Fig. 10), before steadily falling below the ice pressure as the glacier advance slows. The 125m site is an exception to this. Water pressures remain well below the ice pressures after overriding until about day 100, after which the trend follows that of other sites.

4. 4. After about day 150, the glacier reaches its maximum extent, and a quasi-stationary subglacial hydraulic regime is established (Fig. 10).

5. 5. After about day 190, there is a steady decline in aquifer pressure interrupted by two large water-pressure pulses, which appear to be caused by periods of heavy rain (Fig. 10).

The model formalism in section 7 represents a preferred analysis of the glacier hydraulic system in the Fourier domain. However, several issues prevent us from proceeding straightforwardly in this way.

In the following sections, we derive hydraulic properties of the system from the field data by the use of the model (section 8), use the model to explore the implications for basal processes of variations in till conductivity and compressibility (section 9) and of recharge transit times through the glacier (section 10), suggest processes which might account for mismatches between the field data and the model (section 11) and suggest how seismic approaches could be developed further to explore active hydraulic processes at the base of a glacier (section 12).

8.1. Aquifer conductivity, compressibility and dynamics

We assume that the discharge in the aquifer is derived from water that has drained vertically downwards through the till between the local flow divide and the terminus (Fig. 3). The water pressure in the aquifer is determined by the magnitude of this discharge and aquifer transmissivity. This water pressure sets the base pressure for water pressures in the till. The efficiency of this drainage pathway is the major control on till water pressures and effective pressures, and therefore the state of consolidation of the till and the frictional resistance that it offers to glacier movement.

Aquifer properties can most readily be estimated by first considering the quasi-stationary stage. Average daily water pressures in the aquifer are relatively stable, determined primarily by aquifer hydraulic conductivity and by the distance from the glacial terminus to the water divide. Fortunately, these two parameters have different effects on the absolute and relative values of water pressure at different locations. Fitting the predicted quasi-stationary pressure values to measured data from different transducers allows us to recover the distance to the water divide and the conductivity of the aquifer, provided that aquifer thickness (assumed to be 50 m) can be independently estimated. This process indicates a distance to the water divide of ~750 ± 150m from the 125m transducer compared with the measured distance of about 1 km to the transient flow divide on the glacier surface during the mini-surge.

The overall shape of the aquifer pressure curve depends strongly on the hydraulic diffusivity of the aquifer, . If we assume that aquifer properties are constant, we find that we get a good model fit either for the quasi-stationary phase, or for the earlier period when the ice load increases rapidly, but not for both. Figure 14 shows how model 1, which assumes an incompressible aquifer, and has been fitted to the quasistationary stage, predicts aquifer pressures in the rapid loading stage that are much lower than measured. If, however, a consolidation/compressibility effect is included by changing the ratio between two volume compressibility coefficients M_V/m_V (see section 7.2) a surprisingly good model fit is achieved both for the period of rapidly increasing ice pressure, as a result of consolidation of aquifer sediment, and the quasi-stationary phase. Model 2 in Figure 14a–d shows the best-fit results for the 30, 65, 85 and 125m transducers compared with experimental data. These best-fit properties for a 50m thick aquifer, derived using all aquifer transducer data, are shown in Table 1.

Fig. 14. Fitting the models (green and red curves) to aquifer pressure (transducer) data (black curve), (a) for the 125m site, (b) for the 85m site, (c) for the 65m site and (d) for the 30m site. Model 1 predictions (red) ignore aquifer compressibility. Model 2 predictions (green) include a compressibility value of m_V = 2 × 10^–7Pa^–1. The discrepancies between measured and modelled values reveal three important features. Firstly, there is a large mismatch between modelled and measured pressure values during the period of rapid ice loading unless aquifer compressibility is taken into account, and a single compressibility value is a good fit during this phase and the succeeding relatively stable phase. Secondly, after about day 200, the modelled pressure prediction is systematically higher than the measured value, which we interpret as caused by the development of a transverse, low-pressure subglacial channel, up-glacier of the 125m site, that draws down groundwater pressures in the aquifer. Thirdly, the three anomalously high aquifer pressure peaks between about days 160 and 190 at 30m reflect hydrofracturing events in the till very close to this site, which connect the local aquifer very directly with recharge sites at the base of the glacier without the buffering effect of the till.

Table 1. Aquifer characteristics derived by fitting the model to the experimental data. The analysis calculates transmissivity, the product of thickness and conductivity. Conductivity is calculated assuming aquifer thickness to be 50 m. The estimated value of the hydraulic diffusivity (consolidation coefficient, ) and compressibility, m_V, depends on this assumption

Prior to the experiment, we had expected the aquifer to behave passively as a medium of constant transmissibility in discharging the groundwater flux, and did not expect aquifer compressibility to play such an important role in determining aquifer water pressures during ice loading. It is clear, however, that the hydraulic pressure build-up in the aquifer during the first 40 days after glacier overriding is strongly affected by aquifer compressibility, resulting in a cumulative compression of ~1 to 2 m. This enhanced water pressure during ice build-up inevitably leads to increased water pressures in the overlying till, thereby reducing its resistance to glacier movement and potentially enhancing the rate of ice advance. We presume that the very high water pressures during this phase, which can exceed ice pressures at the 65 and 85m transducers, must reflect the absence of easy hydraulic connections above these sites and the glacier surface, which would have released water pressures. This is consistent with the fact that the most rapid phase of loading is associated with the passage of an uncrevassed compressional wave front (Fig. 5a).

It is also particularly important to note that the ratio of the plastic (consolidation) to poroelastic (compressibility) components (M_V/m_V) of 2.6 is very similar to the values observed in laboratory experiments on soil compaction, where M_V corresponds to the virgin consolidation (NCL) part of the deformation curve in Figure 13 (e.g. Reference CraigCraig, 1997). It is normally presumed that there is a severe problem of ‘upscaling’ from laboratory measurements to the field, but the large-scale natural consolidation experiment that we have monitored suggests that, in this particular case, there is no significant up-scaling error.

8.2. Recovering the diurnal component of runoff

There is a strong contrast between the hydraulic signatures in the till and the aquifer. The aquifer naturally averages (filters, suppresses) the high-frequency components of influx whilst till transducers show a highly fluctuating signal at the highest sampling frequency. However, it is difficult to extract information about the dynamics of diurnal pressure fluctuation in the till from our model because we are limited to measurements of surface water production every 24 hours (Fig. 9). We have therefore attempted to recover the diurnal component of runoff using the recorded differences between aquifer and till pressures, given by Equation (26). Knowing the pressure difference across the till and the loading history for the whole dataset allows us to find Fourier components at all frequencies and to resolve Equation (26) with respect to runoff, R. The main difficulty is having to combine experimental results from different transducers, which makes inversion of runoff data non-unique. Two attempts were made to invert runoff data from the 125m transducer for days 0–96.75, from the 85m transducer for days 97–132.75, and from the 30m transducer for the rest of the observation period. It proved impossible, however, to produce a single inversion that matched the runoff data for the whole period. It was possible to obtain a good match between inverted runoff and experimental data for the period after day 137 when the mean value of pressure difference across the till is conserved. We suggest that for the earlier period the movement of crevasses over the site caused local variations in the intensity of runoff (e.g. Fig. 10a, insert). Better results were obtained using synthetic runoff data, where the diurnal fluctuations are added as a simple harmonic signal with maximum at 1800 h:

(38)

where R is the measured (24 hour) water production rate.

8.3. Till conductivity, compressibility and dynamics

A relationship between runoff and water pressures in the till permits us to infer till properties such as hydraulic conductivity, compressibility, the penetration depth of the diurnal runoff pulse and the characteristic timescale of response. The best data for this exercise are those from the 30m transducer for which there is a complete set of till and aquifer pressures. The best match between model and field data is shown in Figure 15. There is a very satisfying agreement for aquifer and till pressures, with discrepancies (1) around day 140, possibly caused by variations in ice profile and velocity that were not picked up by our periodic observations, and (2) between days 165 and 190, possibly the consequence of hydrofracturing events in the aftermath of high-water-pressure events, as well as (3) the general decline in actual compared with modelled pressures after day 200, that could be a consequence of local channel development (see section 11.2). The data do not show a phase delay caused by till filtering, implying that the characteristic time for the till is less than 6 hours, the pressure sampling period. Figure 16 shows similar results for modelling pressure variations in the till at the 65m (mid-till), 85m (mid-till) and 125m upper transducers. There are no obvious phase delays in these data that allow us to put upper limits on the volume compressibility coefficients, m_V, and, consequently, lower limits on the hydraulic diffusivity (consolidation coefficient, c_V). We cannot therefore directly infer the compressibility of the till from these data, but we have assumed that it is similar to that which can be inferred, using our theory, from the data presented by Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others (2001a) where the water-pressure monitoring frequency is sufficient to show depth-dependent phase lags. Figure 17 shows these results, together with a theoretically expected diurnal water-pressure variation and a well-constrained coefficient of volume compressibility, m_V ≈ (7.5 ± 2.5) × 10^–7 Pa^–1. The best fit shown in Figures 15–17 permits us to infer the till properties shown in Table 2.

Fig. 15. (a) Simulated pressures (red curve) for the till transducer at 30m compared with measured pressures (black curve). (b) Simulated and measured pressure drop across the till, showing the maximum gradients at times of maximum water influx (Fig. 9). The till transducer failed to register for about 60 hours after day 220.

Fig. 16. Comparison between mid-till transducer data (black curve) and modelled mid-till water pressures (red curve) at (a) 125 m, (b) 85m and (c) 65 m.

Fig. 17. High-frequency measurements made by till and aquifer transducers reported by Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others (2001a), showing the phase lag between transducers that the mini-surge data cannot show because of the low frequency of measurement. Modelled pressures for the mid-till transducer are also shown.

Table 2. Inferred characteristics of the till at transect sites. The final row shows values calculated from the data of Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others (2001a). The values in parentheses are what would be expected if the compressibility of till were the same as in their article. The inequalities follow from the fact that our measurements are insufficiently frequent to reveal a phase shift at the depths of our transducers. It is noticeable that at 65 and 85 m, the till was sufficiently thick (see Fig. 7) to prevent the diurnal cycle from penetrating to the base of the till

We are surprised by the small range of deduced values of hydraulic conductivity, which we suggest reflects the facts that (i) the till matrix is generally massive and homogeneous, apart from large boulders within it, and lacks the large-scale granulometric variability of stratified sediments, and (ii) the well-mixed sandy/silty matrix favours a narrow range of conductivity compared with clay-rich tills whose clay content often varies considerably over short distances. We are also surprised that trench excavation and refilling appear not to have introduced permeability variations during remoulding, suggesting that trench excavation and refilling did not produce significant sorting and that glacial consolidation and deformation overprinted any consolidation during trench refilling.

Laboratory consolidation tests conducted on a basal till sample collected from a trench dug into the margin of Storglaciaren, northern Sweden, by Reference Baker and HooyerBaker and Hooyer (1996) gave estimates of compressibility coefficients of 1–5 × 10^–7 Pa^–1 in the effective pressure range 100–300kPa. This agrees well with our mean in situ average of ≈7.5 × 10^–7 Pa^–1 under effective pressures between 0 and 50 kPa.

Using the till thickness along the transect from Figure 7, interpolating between the values of conductivity shown in Table 2 for the transducer sites and using a constant value of compressibility (m_V = 7.5 × 10^–7Pa^–1), we have modelled the till response for every 5m interval along the transect. Figure 18 shows the results for hydraulic diffusivity and characteristic time. It shows the effect of the finer-grained till characteristic of the inner part of the transect (Fig. 8) in slowing the response time. In principle, the till in the inner part of the transect will show a greater tendency to undrained loading and instability than that in the outer part.

Fig. 18. Inferred till properties along the transect. Conductivity and compressibility were deduced by theory from the water-pressure and water-production data; these were used in turn to infer the penetration depth of the diurnal wave and the characteristic penetration time. Properties in (b-e) reflect the finer grain size of the till in the inner part of the transect (see Fig. 8).

11.1. Self-adjustment of aquifer transmissivity through hydrofracturing

The transducers at 125m exhibit intriguing behaviour that we are unable to explain using constant values of aquifer properties. Even though our model 2 (section 8.1 and Fig. 14) fitted the quasi-stationary phase, which incorporates aquifer compressibility, and predicts aquifer pressures well for the period of rapid ice-load increase between about days 85 and 150, there remains a poor fit between about days 37 and 85. During this period, the aquifer transducer shows a very significantly higher water pressure than predicted by either model 1 or 2 (Fig. 14a). The period prior to day 37 is complicated by the upward flow of water near the margin (see section 7.7). We can only interpret this anomaly as reflecting a strong increase in transmissivity between about days 75 and 85.

This becomes plausible if we consider the field observation that fine-grained silty beds commonly occur between gravel lenses within the sub-till outwash sequence (Fig. 4). The change in apparent aquifer properties occurs at about the time of several large runoff-generated pressure spikes. If the effective transmissivity of the aquifer unit at about 125m was limited because of silt aquitards between gravel aquifer compartments, a major water-pressure spike at the ice-bed interface would pressurize the immediately underlying aquifer. This would generate a large potential gradient across the silt bed separating it from the underlying unpressurized aquifer. Consequent hydrofracturing across the silt bed would create new hydraulic pathways, immediately increasing the transmissivity of the aquifer. This might explain the change at about days 75–85, during which newly developing pathways progressively changed aquifer properties. We find that we are able to model the water-pressure behaviour of the aquifer before and after this period if we assume an effective aquifer thickness of 15m beforehand and 50m afterwards (the latter being the thickness inferred from drilling and seismology). We suggest that the apparent increase in effective hydraulic thickness results from the hydrofracturing of conductivity barriers that allow hydraulic compartments to be added to the aquifer.

11.2. Groundwater pressure drawdown due to subglacial channel formation

During days 192–242 (Fig. 14) there is disagreement between model and experiment for the 12, 30 and 65m aquifer transducers (the 125 and 85m transducers were no longer operational). The model predicts a slower general trend of pressure decrease than those clearly seen at the 30 and 65m transducers (Fig. 14c–d) and a higher response to prominent peaks in runoff data. It is impossible to satisfy both data trends because the only physical parameter that strongly affects the aquifer response is the characteristic response time of the aquifer. Decreasing this parameter to account for a rapid seasonal decline in pressures will automatically increase the sensitivity of the aquifer to (i) short influx pulses and (ii) the build-up of water pressure during rapid initial ice loading, thereby increasing the amplitudes of the two prominent peaks.

We suggest that this discrepancy can be explained by local shifts in the mode of subglacial drainage. The retreat of the glacier after the mini-surge ultimately led, in the early 1990s, to a direct connection between Stemmulon lake to the northeast of the experimental site (Fig. 1) and Jökulsárlón, into which it drained. We suggest that across-flow drainage into Jökulsárlón started during this stage, as the glacier stagnated. The anomalously low water pressures during this period could thus be explained by a subglacial channel developing, that ran across the glacier flowline, upglacier of the monitored transect, and intercepted a significant proportion of the aquifer flow, thereby drawing down heads. This would both produce a lower than expected head in the aquifer and a smaller aquifer response to the major rainfall events of days 216 and 234.

11.3. Transient upward water-pressure gradients in the till

Upward head gradients were measured far from the glacier margin during brief diurnal periods. Figure 11b shows the normal state with a downward head gradient and downward flow through the till, with an insignificant diurnal pressure fluctuation in the aquifer. Figure 11c shows a period of strong variations in aquifer pressure, when aquifer pressure peaks temporarily exceed those in the till for periods between 6 and 12 hours in the latter part of the day. In general, the aquifer is recharged by a downward water flux derived from fractures, moulins and crevasses in the glacier. We suggest, however, that where there are no fractures or crevasses in overlying basal ice to provide proximal local recharge at the top of the till, strong increases in aquifer pressure produced by strong recharge from nearby but not overlying crevasses create a transient upward gradient, which decays as aquifer pressures fall again.

11.4. Till hydrofracturing during strong recharge/high-pressure events

The water-pressure gradient across the till normally falls within the range 5–60 kPam^–1. However, during some periods of high recharge from the surface, generally associated with heavy rainfall, this can increase to over 100 kPa m^–1. On a number of such occasions pressure peaks in the till then decline sharply to values very similar to those in the aquifer (e.g. Fig 11d). They are then sustained at such levels for several days, although showing some fluctuations, before previous values of pressure gradient in the till are progressively re-established.

As has been argued for silt beds in the aquifer in section 8.2, we suggest that the large pressure differences between the till and aquifer, produced by strong recharge from the surface, create hydrofractures in the till. These produce efficient hydraulic pathways between the ice-bed interface and aquifer, which rapidly draw down pressures on the till surface and therefore in the till. This locally enhances the rate of recharge to the aquifer, increasing aquifer pressures. If hydrofracturing is widespread, a general increase of water pressure in the aquifer would result, giving enhanced peaks of diurnal aquifer pressure. If till deformation continued during this phase, the break in till cover produced by hydrofracturing would be progressively sealed, re-establishing the previous hydrological regime. In general, therefore, we might expect hydrofracturing of thin till to be a means of inhibiting large water-pressure differences across the till.

11.5. Water-pressure effect of shearing in the till

Till transducer records tend to show a strong fall in water pressure immediately before they cease recording (Fig. 10). We suggest this reflects dilation associated with shearing in the till, which increases conductivity and reduces water pressure. We believe that transducers cease recording shortly after the till in which they are embedded begins to suffer significant shear deformation. Our records, then, will relate only to non-deforming conditions. This assumption underlies our approach. When deformation of till occurs, we suggest that the potential gradient in the deforming horizons becomes relatively small as a consequence of increased conductivity due to dilation (Reference Boulton, Menzies and RoseBoulton and Hindmarsh, 1987; Reference Boulton and DobbieBoulton and Dobbie, 1993), thereby reducing the effective thickness of the till to which our analysis applies.

When the uppermost of two till transducers fails, pressure at the upper transducer falls below that of the lower. We suggest that as dilation increases and water pressures fall in the upper part of the till as a consequence of shear deformation, pressures in the undeforming body of the till remain high, so that water drains upwards towards the glacier sole. A reversal of drainage that we suggest reflects this process is seen in Figure 17 late on days 253 and 254, when contemporary monitoring of strain (Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others, 2001b) confirms that the upper part of the till is deforming.

We have attempted to infer the properties of the till along the whole of the transect by considering changes in till thickness and surface elevation, and by assuming that hydraulic conductivity and compressibility can be interpolated between the values at measured sites (Fig. 18a-c). We have then determined the penetration depth for the diurnal wave and the characteristic response time for the till (Fig. 18d–e). Penetration depth shows a pattern that is consistent with grain-size variation along the transect (Fig. 8), with the finer-grained till in the up-glacier part of the transect having a smaller penetration depth and higher characteristic time than the outer part. This part of the till will also be more prone to undrained loading and consequent potential instability.

13.1. Analysis in the frequency domain

1. 1. In order to understand the frequency response of the system to changes in water recharge and ice loading, the hydraulic pathway from the glacier surface to the bed and through subglacial sediments is modelled in the Fourier domain within a framework of three subsystems (ice-till-aquifer) of different linear poroelastic properties, coupled at the ice-till and till-aquifer interface, and with one-dimensional Darcy diffusion of water; it is solved analytically for arbitrary runoff and loading histories. This allows processes to be reformulated in the Fourier frequency domain and resolved analytically at each given location for arbitrary runoff and ice-loading histories. It permits the response of system properties to a broadband signal and changes inboundary conditions to be readily and directly determined, in contrast to the merely ad hoc testing of system properties that can be achieved by finite-element or -difference schemes. The model permits:

   the deduction of in situ values of both compressibility and hydraulic conductivity of the till and aquifer directly from the transducer data,

   the determination of the ratio of plastic consolidation poroelastic compressibility, which is shown to be similar to that measured in laboratory tests, indicating absence of the normally assumed ‘up-scaling’ problem.

2. 2. Major complications are caused by the non-linearity of sediment consolidation (Fig. 13) and by evolving ice-till boundary conditions. We model these non-linear responses of till and aquifer to effective stress with an essentially linear model by considering the frequency-dependent compressibility of sediments (section 7.2.1). The problem of the advancing glacier is resolved by ‘slicing’ the system into a number of vertical cells connected only by horizontal water flow in the aquifer and with independently evolving boundary conditions at the top of the till, caused by the changing glacier load.

3. 3. It is unfortunate that the characteristic frequencies of the till and glacier ice were outside (higher than) the 6 hour data frequency available from most transducer signals. This is because of the relatively high conductivity and low compressibility of the coarse-grained till at Breiðamerkurjökull. However, the Fourier approach did demonstrate the existence of a quite unexpected aquifer-scale compressibility that created an important water-pressure wave beneath the advancing glacier terminus.

4. 4. We have used our model to explore how fine-grained tills of low hydraulic conductivity, high compressibility and/or great thickness would behave in response to variations in recharge from the glacier. Lower conductivity and higher compressibility inhibit penetration into the till of water-pressure fluctuations so that the signal attenuates at shallow depth, reducing the amplitude of oscillation, increasing the penetration time, increasing the phase lag of the pressure signal with depth and creating the potential for mechanical instability associated with undrained loading, including during rapid increases in loading due to advance of a steep glacier front. The rapid rate of loading associated with the advance of a steep glacier front produces anomalously high water pressures in aquifers as they consolidate, thereby driving up the water pressure in any overlying till. This can be a source of positive feedback to the forward movement of a glacier, and may in part account for the mini-surge. Highly compressible tills may also show transient upward water-pressure gradients as the recharge rate on their upper surface falls. An oscillating water pressure on their upper surface produces periodic, though phase-lagged, reversal of the water-pressure gradient, and could explain periodic oscillations of the décollement surface in tills (e.g. Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others, 2001a).

5. 5. The model is also used to explore the impact of drainage rate of surface meltwater through a glacier on subglacial hydraulics and dynamics. At the Breiðamerkurjökull site, the drainage time is less than the 6 hour frequency of most transducer records, so we are unable to capture its influence on the frequency response of the subglacial hydraulic system. We assume that the drainage time increases under thicker ice, and model how high-frequency components of change in the rate of surface water production will be progressively attenuated under thicker ice. In a near-margin zone, diurnal surface melt fluctuations are reflected in the subglacial hydraulic system. Under increasingly thick ice, higher-frequency fluctuations were progressively attenuated in the subglacial response, so that diurnal responses are progressively replaced by dominantly seasonal, annual and multi-year frequencies. We presume that at some ice thickness, or because of cold ice conditions, there will be negligible surface water penetration to the bed and negligible melt-driven subglacial hydraulic effects, such that water pressures will only be determined by basal melt rates and the drainage pathway. In such a system, the time-scales of water-pressure response will be major determinants of basal dynamic processes. We might therefore expect higher-frequency dynamic effects in marginal zones and lower-frequency effects in inner zones of the same ice sheet. Understanding the whole hydraulic pathway from surface to bed and the effective conductivity of glacier ice is a key issue in understanding the way in which ice sheets may respond dynamically to a warming climate through processes similar to those monitored at Breiðamerkurjökull.

13.2. Analysis in the physical domain

Characteristics of the water-pressure record, and the extent to which they diverge from expectations of the model, are evidence of other important hydraulic processes.

1. 1. The surface water production rate, transmissivity of the drainage pathway and fluctuations of ice loading control the state of consolidation of subglacial sediments. The time-dependent variation of ice and water pressure determine the changing state of consolidation of subglacial sediments. The complex hydraulic history at Breiðamerkurjökull produces a complex consolidation history, with previous preconsolidation histories being removed by shear remoulding if rising water pressures exceed a critical value. We deduce from these data that heavy pre-consolidation in a till does not preclude the till having undergone shear deformation at some time during its history.

2. 2. It is deduced from the model that a shift in the water-pressure regime during the experiment can only be explained by a major increase in aquifer transmissivity. It is suggested that this occurs because of hydrofracturing in a sequence in which individual lenticular aquifers, reflecting the geometry of glaciofluvial fans, are separated by finegrained, less permeable beds. Rapid loading produced transient strong pressurization of immediately subglacial aquifer lenses, which generated hydrofractures in the finergrained, less permeable beds which separate the upper aquifers from lower, relatively unpressurized aquifers. Hydraulic connection between the aquifers equalized the water pressures between them, thereby increasing the connected transmissivity of the subglacial hydraulic system and producing general water-pressure drawdown.

3. 3. Convergence of the pressures in the till and the top of the aquifer occurred immediately after high till pressure peaks at several sites on several occasions, even though diurnal runoff peaks remained high. It is suggested that rapid water-pressure increase on the till top created hydrofractures in the till, making easy hydraulic connections between the glacier sole and the underlying aquifer, thus locally lowering till pressures so that they come close to aquifer pressures. Subsequent shear deformation of till sealed these interconnections, permitting till-top water pressures again to rise well above aquifer pressures.

4. 4. Upwelling occurs immediately beyond the margin of the advancing glacier, with water rising through the till from the underlying aquifer driven by a locally strong upward head gradient of up to 25 kPam^–1. Observations elsewhere at the margin show that upwelling can be locally strong enough to cause liquefaction and buoyancy in near-surface sediments. Such upwelling, and the reduction in proglacial sediment strength that it causes, may also have facilitated the creation and easy gliding of a push moraine at the advancing glacier terminus. Van der Reference Van der Meer, Kjaer and KrugerMeer and others (1999) have described ‘water escape structures’ in a formerly subglacial aquifer from the glacier Slettjökull in southern Iceland, which we suggest may be former locations of hydrofractures created by pressurized sub-till aquifers.

5. 5. The measured water pressures are much lower than those predicted by the model at a late stage in the experiment; this is suggested to be the result of formation of a subglacial drainage channel near the experimental site, which led to general drawdown of groundwater pressures.

6. 6. The Breiðamerkurjökull tills are similar in granulometry to most tills lying on the ancient shield rocks in the central areas of the European and Laurentide ice sheets. Our model shows that such tills of high conductivity and low compressibility adjust rapidly to changes in local pressure, such that instabilities are rapidly attenuated, unless the till is very thick. Tills overlying soft sedimentary rocks in the outer zones of the mid-latitude ice sheets show different behaviour, of low conductivity and high compressibility. Modelling the behaviour of such tills demonstrates their potential for undrained loading, long response times and sustained instability. The highly lobate form of ice margins in these zones (Reference Boulton, Hagdorn and HuttonBoulton and others, 2003) may reflect such instabilities.

7. 7. The poroelastic theory used in the model to analyze hydraulic effects can also be used to predict seismic responses. We suggest that, building on existing work (e.g. Reference SmithSmith, 1997; Reference AnandakrishnanAnandakrishnan, 2003; Reference Vaughan, Smith, Nath and Le MeurVaughan and others, 2003), the theory could also be used to infer spatial and temporal changes in sediments at the ice-bed interface, including porosity, compressibility and possibly conductivity.