Subglacial hydrology plays a key role in many glaciological processes, including ice dynamics via the modulation of basal sliding. Owing to the lack of an overarching theory, however, a variety of model approximations exist to represent the subglacial drainage system. The Subglacial Hydrology Model Intercomparison Project (SHMIP) provides a set of synthetic experiments to compare existing and future models. We present the results from 13 participating models with a focus on effective pressure and discharge. For many applications (e.g. steady states and annual variations, low input scenarios) a simple model, such as an inefficient-system-only model, a flowline or lumped model, or a porous-layer model provides results comparable to those of more complex models. However, when studying short term (e.g. diurnal) variations of the water pressure, the use of a two-dimensional model incorporating physical representations of both efficient and inefficient drainage systems yields results that are significantly different from those of simpler models and should be preferentially applied. The results also emphasise the role of water storage in the response of water pressure to transient recharge. Finally, we find that the localisation of moulins has a limited impact except in regions of sparse moulin density.

The most frequently used relation between ice deformation rate, έ, and stress, ፐ, is the power law, commonly called Glen’s flow law, έ = Aፐn , in which A is an ice stiffness parameter and n is an empirical constant. A can be estimated from the simple exponential relation where A 0 is a constant independent of temperature; E, commonly called the enhancement factor, depends on ice crystal orientation, impurity content and other factors; Q is the activation energy for creep; R is the universal gas constant; and T is the absolute temperature. Laboratory experiments yield values of A 0 = 9.514MPa–3a–1 for secondary creep. Typical borehole closure experiments then give E ≈ 0.16. This low value probably results from the fact that, when deforming into a borehole, ice is subject to stresses that are inconsistent with the preferred orientation of c axes that has developed over many years under a stress configuration with no borehole present. Closure data from Vostok hole 3G yield E ≈ 0.7. This higher value may reflect a unique stress environment yielding fabrics that are somewhat better oriented for borehole closure.

Annual-layer thickness data, spanning AD 1534–2001, from an ice core from East Rongbuk Col on Qomolangma (Mount Everest, Himalaya) yield an age–depth profile that deviates systematically from a constant accumulation-rate analytical model. The profile clearly shows that the mean accumulation rate has changed every 50–100 years. A numerical model was developed to determine the magnitude of these multi-decadal-scale rates. The model was used to obtain a time series of annual accumulation. The mean annual accumulation rate decreased from ∼0.8 m ice equivalent in the 1500s to ∼0.3 m in the mid-1800s. From ∼1880 to ∼1970 the rate increased. However, it has decreased since ∼1970. Comparison with six other records from the Himalaya and the Tibetan Plateau shows that the changes in accumulation in East Rongbuk Col are broadly consistent with a regional pattern over much of the Plateau. This suggests that there may be an overarching mechanism controlling precipitation and mass balance over this area. However, a record from Dasuopu, only 125 km northwest of Qomolangma and 700 m higher than East Rongbuk Col, shows a maximum in accumulation during the 1800s, a time during which the East Rongbuk Col and Tibetan Plateau ice-core and tree-ring records show a minimum. This asynchroneity may be due to altitudinal or seasonal differences in monsoon versus westerly moisture sources or complex mountain meteorology.

At Engabreen, Norway, an instrumented panel containing a decimetric obstacle was mounted flush with the bed surface beneath 210 m of ice. Simultaneous measurements of normal and shear stresses, ice velocity and temperature were obtained as dirty basal ice flowed past the obstacle. Our measurements were broadly consistent with ice thickness, flow conditions and bedrock topography near the site of the experiment. Ice speed 0.45 m above the bed was about 130 mm d−1, much less than the surface velocity of 800 mm d−1. Average normal stress on the panel was 1.0–1.6 MPa, smaller than the expected ice overburden pressure. Normal stress was larger and temperature was lower on the stoss side than on the lee side, in accord with flow dynamics and equilibrium thermodynamics. Annual differences in normal stresses were correlated with changes in sliding speed and ice-flow direction. These temporal variations may have been caused by changes in ice rheology associated with changes in sediment concentration, water content or both. Temperature and normal stress were generally correlated, except when clasts presumably collided with the panel. Temperature gradients in the obstacle indicated that regelation was negligible, consistent with the obstacle size. Melt rate was about 10% of the sliding speed. Despite high sliding speed, no significant ice/bed separation was observed in the lee of the obstacle. Frictional forces between sediment particles in the ice and the panel, estimated from Hallet’s (1981) model, indicated that friction accounted for about 5% of the measured bed-parallel force. This value is uncertain, as friction theories are largely untested. Some of these findings agree with sliding theories, others do not.

Assuming a channelized drainage system in steady state, we investigate the influence of enhanced surface melting on the water pressure in subglacial channels, compared to that of changes in conduit geometry, ice rheology and catchment variations. The analysis is carried out for a specific part of the western Greenland ice-sheet margin between 66° N and 66°30′N using new high-resolution digital elevation models of the subglacial topography and the ice-sheet surface, based on an airborne ice-penetrating radar survey in 2003 and satellite repeat-track interferometric synthetic aperture radar analysis of European Remote-sensing Satellite 1 and 2 (ERS-1/-2) imagery, respectively. The water pressure is calculated up-glacier along a likely subglacial channel at distances of 1, 5 and 9 km from the outlet at the ice margin, using a modified version of Röthlisberger’s equation. Our results show that for the margin of the western Greenland ice sheet, the water pressure in subglacial channels is not sensitive to realistic variations in catchment size and mean surface water input compared to small changes in conduit geometry and ice rheology.

The University of Maine Ice Sheet Model was used to study basal conditions during retreat of the Laurentide ice sheet in Maine. Within 150 km of the margin, basal melt rates average ∼5 mm a−1 during retreat. They decline over the next 100 km, so areas of frozen bed develop in northern Maine during retreat. By integrating the melt rate over the drainage area typically subtended by an esker, we obtained a discharge at the margin of ~1.2 m3 s-1. While such a discharge could have moved the material in the Katahdin esker, it was likely too low to build the esker in the time available. Additional water from the glacier surface was required. Temperature gradients in the basal ice increase rapidly with distance from the margin. By conducting upward into the ice all of the additional viscous heat produced by any perturbation that increases the depth of flow in a flat conduit in a distributed drainage system, these gradients inhibit the formation of sharply arched conduits in which an esker can form. This may explain why eskers commonly seem to form near the margin and are typically segmented, with later segments overlapping onto earlier ones.

Empirical data suggest that the rate of calving of grounded glaciers terminating in water is directly proportional to the water depth. Important controls on calving may be the extent to which a calving face tends to become oversteepened by differential flow within the ice and the extent to which bending moments promote extrusion and bottom crevassing at the base of a calving face. Numerical modelling suggests that the tendency to become oversteepened increases roughly linearly with water depth. In addition, extending longitudinal deviatoric stresses at the base of a calving face increase with water depth. These processes provide a possible physical explanation for the observed calving-rate/water-depth relation.

The temperature distribution in a polar glacier is described by the equation of heat conduction,

1

where K is the thermal diffusivity of ice, Q is the internal heat generation, p is the ice density, and C is the heat capacity. To obtain a solution to this equation, boundary conditions at the surface and bed must be known. The boundary condition at the bed is generally taken to be the temperature gradient in the ice required to conduct the geothermal heat upward into the glacier, with certain modifications where the pressure melting temperature is reached. The boundary condition at the surface is the ice temperature, which is usually assumed to be equal to the mean annual atmospheric temperature. This assumption is incorrect in the ablation area and in the percolation and saturation zones of the accumulation area. In this paper I examine the reasons for the break down of this assumption, and attempt to indicate the magnitude of the error introduced.

The atmospheric temperature at a glacier surface changes seasonally; thus measurements of the “surface” temperature for use with Equation (1) are generally made at some depth, z 0, in the glacier below which the effect of these seasonal variations is negligible. If the seasonal variation can be represented by a sinusoidal function, this depth is given by: 2

(Carslaw and Jaeger, 1959, p. 65) where w is the period of the fluctuations (in this case 2π/year), θr is the temperature range from winter minimum to summer maximum, and A is the maximum acceptable change in temperature at depth z 0. For example, in the dry zone of the accumulation area if we take K = 16 m2/year, a value appropriate for unpacked snow, θ r = 30 deg, and Δ = 0.4 deg, we obtain z 0 = 10 m. This is the basis for the common assumption that the 10 m temperature is approximately equal to the mean annual temperature.

In the superimposed ice zone superimposed ice occurs immediately beneath the winter snow cover, and K for ice at — 10°C is about 38 m2/year. Furthermore, the temperature fluctuation at the ice—snow interface cannot be approximated by a sinusoidal function because the accumulating snow cover insulates the ice during the winter, and because the temperature rises rapidly in the late spring when percolating melt water reaches the interface. Equation (2) can still be used to calculate an approximate value for z 0 of about 15 m, but due to the non-sinusoidal temperature variation at the ice-snow interface, the temperature at depth is commonly a few degrees above the mean annual atmospheric temperature.

The magnitude of this difference, which we will call,Δθ can be calculated from Equation (1) if it is assumed that convection, internal heat generation, and transverse and longitudinal conduction are negligible, and if the proper boundary condition at the ice-snow interface is known. The equation to be solved is:

3

The boundary condition at depth z0 is taken to be the temperature gradient below this depth. The problem is thus reduced to one of determining the temperature of the ice-snow interface as a function of time,θs(t.

In the soaked zone of the snow cover rests on permeable firn rather than on superimposed ice, and melt water percolating down through the snow pack can penetrate some distance into the firn. Upon refreezing, this water releases the heat of fusion, thus warming the firn. Mathematically, this can be represented by adding an internal heat-production term to Equation (I), thus:

4

Due to this internal heat production, Δθ may be substantially larger in the soaked zone than in the superimposed ice zone or ablation area.

In order to determine θs(t) and Q(z), six 30 m bore holes were drilled on the south dome of the Barnes Ice Cap. Five of these holes were along a flow line extending from the divide to the margin. Temperature measurements were made in each hole in mid-July 1973, and three times in June and July 1974. Two finite-difference calculations were carried out with the use of these data. In one, the first of the 1974 measurements was used as an initial condition, and Equation (1) was integrated, using an assumed form of θs(t) as a boundary condition. The form of θS(t) was varied until reasonable agreement was obtained between calculated and measured profiles at the times of the second and third measurements in 1974. This calculation thus permits an estimate of θs(t) during the critical period of melt from early June to mid-July.

In the second calculation, the July 1973 temperature profile was used as an initial condition, and Equation (1) was again integrated, this time simulating a time period of one year. In this model, snow was allowed to accumulate on the ice surface. The accumulation pattern of the snow cover as a function of time was determined from climatic records from nearby weather stations. The temperature variation at the snow-air interface was also determined from these records, using measured lapse rates to correct for differences in altitude (private communication from R. Barry in 1974). Calculations using this model suggest that three factors have a significant influence on the temperature distribution in the ice. One is the fact that most accumulation occurs in the fall and spring, with relatively little accumulation during the winter. The second is that thermal diffusivity of the snow apparently increases during the early spring as atmospheric temperatures begin to rise. The third and most important factor is that the snow-ice interface remains at 0°C late into the fall due to the presence of melt water that percolated down through the accumulating snow cover in the early fall, and the temperature of this interface rises to 0°C early in the spring, again due to percolation of melt water.

The temperature measurements and climatic records were also used to estimate values for Δθ on the Barnes Ice Cap. In the ablation area and superimposed ice zone Δθis 2 to 4 deg. Then there is a rapid increase to about 5 deg in the lower part of the soaked zone. This jump occurs over a distance of only a couple of kilometers across the boundary between the two zones.

Perhaps the most significant consequence of this increase in surface temperature in the superimposed ice and soaked zones is that the temperature throughout this part of a glacier will also increase by approximately the same amount. Thus if the base of a glacier were not already at the pressure melting temperature up-glacier from the soaked zone, it is likely that it would rise to the pressure melting temperature somewhere beneath the soaked or superimposed ice zones. The consequences of this for basal erosion and the entrainment of morainal material need to be examined.

The structure and flow field in the margin of the Barnes Ice Cap was determined through observations on the ice-cap surface, in four bore holes, and in a 125 m ice tunnel. A band of fine bubbly white ice with a single maximum fabric appears at the glacier surface about 160 m from the margin. This band is overlain by coarse blue ice with a four-maximum fabric, and underlain by alternating bands of fine ice with a single-maximum fabric and moderately coarse ice with a two or three-maximum fabric.

The effective strain rate was determined from the bore-hole and tunnel deformation data, and possible variations in the other three parameters in Glen’s flow law, , were studied. It appears that τ xy is independent of depth near the surface, and that relative to the coarse blue ice, A is 40 to 50% lower in the white ice and possibly 10% lower in the fine blue ice.

Dips of foliation planes decrease rapidly with increasing depth and distance from the margin. This foliation is assumed to have developed near and parallel to the bed some distance from the margin. An analysis based on this assumption predicts the observed change in dip, but suggests that it did not develop under the present flow field. The ice cap was probably thicker a few tens of years ago, and the observed foliation pattern may be a relict from that time.

Three types of glacier margin are found along the edge of the Greenland ice sheet near Thule: ice cliffs, ramps and ice-cored moraines. Where the glacier margin is perpendicular to prevailing katabatic winds, drifting snow accumulates along it in stagnant wind-drill ice wedges. Upward flow of active ice behind these wedges causes ice originally near the base of the glacier to rise to the surface. Where this basal ice is free of debris, a gently sloping ramp develops. However, where the basal ice contains sufficient debris, a layer of till accumulates on the glacier surface. Ice beneath the till is insulated and a debris-capped ice ridge or ice-cored moraine forms, Ice cliffs occur where the ice-sheet margin is parallel to prevailing winds and is thus swept clear of drifting snow. Although the ice sheet in the Thule area appears to have had a negative mass balance for many years, all three types of glacier margin are believed to be equilibrium forms that can develop and persist on a glacier with a balanced mass budget.

Foliation in wind-drift ice wedges generally dips down-glacier but foliation in active ice dips up-glacier. It is inferred that foliation in the wedges was once sedimentary stratification that has been tipped upward and locally overturned.

Numerical methods based on quadrilateral finite elements have been developed for calculating distributions of velocity and temperature in polar ice sheets in which horizontal gradients transverse to the flow direction are negligible. The calculation of the velocity field is based on a variational principle equivalent to the differential equations governing incompressible creeping flow. Glen’s flow law relating effective strain-rate ε̇ and shear stress τ by ε̇ = (τ/B) n is assumed, with the flow law parameter B varying from element to element depending on temperature and structure. As boundary conditions, stress may be specified on part of the boundary, in practice usually the upper free surface, and velocity on the rest. For calculation of the steady-state temperature distribution we use Galerkin’s method to develop an integral condition from the differential equations. The calculation includes all contributions from vertical and horizontal conduction and advection and from internal heat generation. Imposed boundary conditions are the temperature distribution on the upper surface and the heat flux elsewhere

For certain simple geometries, the flow calculation has been tested against the analytical solution of Nye (1957), and the temperature calculation against analytical solutions of Robin (1955) and Budd (1969), with excellent results.

The programs have been used to calculate velocity and temperature distributions in parts of the Barnes Ice Cap where extensive surface and bore-hole surveys provide information on actual values. The predicted velocities are in good agreement with measured velocities if the flow-law parameter B is assumed to decrease down-glacier from the divide to a point about 2 km above the equilibrium line, and then remain constant nearly to the margin. These variations are consistent with observed and inferred changes in fabric from fine ice with random c-axis orientations to coarser ice with single- or multiple-maximum fabrics. In the wedge of fine-grained deformed superimposed ice at the margin, B increases again.

Calculated and measured temperature distributions do not agree well if measured velocities and surface temperatures are used in the model. The measured temperature profiles apparently reflect a recent climatic warming which is not incorporated into the finite-element model. These profiles also appear to be adjusted to a vertical velocity distribution which is more consistent with that required for a steady-state profile than the present vertical velocity distribution.

Laboratory studies suggest that neither bubbles nor dirt particles migrate rapidly enough in glacier ice to be responsible for the alternating layers of bubbly and clear ice or dirty and clean ice which constitute foliation. We therefore suggest that these variations in bubble or dirt content are inherited from primary inhomogeneities such as may occur in sedimentary stratification in the accumulation region, in crevasse fillings, or during debris entrainment at the base of the glacier: the appearance of these inhomo-geneities is later modified by strain during flow to produce foliation. We consider six types of inhomogeneity, or components of foliation, and show that, at the very large total strains expected in glaciers, all are eventually flattened, stretched out, and rotated to form a layered structure roughly perpendicular to the direction of maximum total shortening. Most characteristics of observed foliation can be explained by this hypothesis. For example, in the marginal zones of polar ice sheets the rapid decrease in dip of foliation with depth and with distance up-glacier from the margin can he explained by a model in which the foliation is assumed to be nearly parallel to the base of the glacier some distance from the margin, and is deformed passively with the ice thereafter. However, some observations of cross-cutting foliations may require localized inhomogeneous shear parallel to the "new" foliation.

Oxygen-isotope ratios indicate that a distinctive band of white ice along the margin of the Barnes Ice Cap is of Pleistocene age. It is estimated freom a flow mogerl that beneath the center of the ice cap the thickness of the band shuuld be abont 0.6 times its thickness at the margin, or about 8 m. However, an ingerpengernt estimate, based on calculated vertical strain-rates and explicitly assuming no basal melting, predicts a thickness of about 22 m beneath the center of the ice cap. The discrepancy between the two thickness estimates is interpreted as indicating that basal melting has occurred. Calculated basal temperatures support this conclusion.

Cylindrical samples of ice with 0.0 to 0.35 volume fraction fine sand were tested in unconfined uniaxial compression at stresses between 5.3 and 6.4 bar and at temperatures between −7.4 and −9.4° C. Secondary creep rates were obtained from the slope of the total strain vs. time curve and were normalized to 5.6 bar and −9.1° C. Creep rates in ice with low sand concentrations were in some cases higher and in other cases lower than in clean ice. However at higher sand concentrations the creep rate decreases exponentially with increasing volume fraction sand. The latter results are in general agreement with theories developed to explain dispersion hardening of metals, and suggest that each sand grain is surrounded by a tangled network of secondary dislocations which impede passage of primary glide dislocations.

Between 3 June 1982 and 8 July 1985, a stake net consisting of up to 32 stakes covering the greater part of Storglaciären was surveyed 70 times, yielding roughly 2000 separate determinations of vertical and horizontal velocity. The time interval between surveys averaged about 1 week during the summer and 2 months during the winter.

Horizontal velocities were normally highest during periods of high daily temperature or heavy rain early in the melt season. Comparable or sometimes higher temperatures or rainfalls later in the season usually had less effect, though minor velocity peaks were often present in August and early September. During periods for which bore-hole water-level measurements are available, velocity peaks generally coincided with periods of high basal water pressure, but not all periods of high water pressure resulted in velocity peaks. Despite increasing basal water pressures, velocity decreased gradually during the winter.

Vertical velocities also vary seasonally. Beneath the upper part of the ablation area the glacier bed is overdeepened. Vertical velocities here are ˜3 mm/d higher during the summer. Down-glacier from the overdeepening, vertical velocities are ˜1 mm/d lower during the summer. These and other characteristics of the vertical velocity pattern are best explained by appealing to: (1) a decrease in strain-rate with depth, and (2) seasonal variations in this depth-dependence.

Five periods of high velocity lasting from 3 to 11d were studied in detail. In an area where the bed is overdeepened, force-balance calculations suggest that basal drag decreased between 16 and 40% during these high-velocity events. This resulted in a decrease in compressive strain-rate at the up-glacier end of the overdeepening, an increase at the down-glacier end, and a slight increase in lateral shear strain-rates. Down-glacier from the overdeepening, basal drag increased during two events owing to an increased push from up-glacier and pull from down-glacier. Lateral shear strain-rates increased sharply here.