Introduction
Fluctuations in thickness or lateral extent of a glacier may be caused by climatic forcing or by internal instabilities. It is becoming increasingly evident that glaciers resting on deformable beds can exhibit internal instabilities and large amplifications of small climatic forcing. An understanding of the glacier-bed system will be necessary to allow interpretation of past and ongoing glacier changes in terms of climatic forcing, and to predict future ice-sheet behavior forced by anthropogenic climatic change. Here I report progress of an ongoing program of coupled glacier-bed modeling to address this problem.
Model
The model used here is based on several observations of glaciers supported by results from theory (Reference WeertmanWeertman, 1972; Reference Engelhardt, Harrison and KambEngelhardt and others, 1978; Reference Boulton and HindmarshBoulton and Hindmarsh, 1987; Reference AlleyAlley, 1989):
– the upper surface of a typical glacier till is rough to a glacier and causes basal sliding to be slow;
– deforming basal till with high water pressure is orders of magnitude softer than basal ice;
– deforming till will creep into low-pressure regions at the glacier bed;
– a connected water system at the ice-bed interface will have high water pressure in the absence of well-developed channel systems;
– natural glacier beds range from unconsolidated with low cohesion (will deform under glaciogenic stresses if water pressure sufficiently high) to strong bedrock (not deformable under glaciogenic stresses, and difficult to erode).
Taken together, these conditions suggest that incipient drainage channels fed by basal water beneath a wet-bedded glacier that does not lose its water to subglacial aquifers and that generates basal till will be closed by till creep. A glacier resting on any till then will exhibit rapid basal velocities if, and only if, the bed is easily eroded and produces a thick basal till.
A model for cold ice based on these concepts can be stated in words as follows (equations are given in the appendix):
– a glacier frozen to its bed moves entirely by internal ice shearing (deforming permafrost at the bed is defined as glacial ice for convenience); if a glacier bed is at the pressure-melting point, then the ice velocity is the sum of the basal and ice-deformational velocities;
– water produced at the bed first saturates any porosity produced by ongoing erosion, and then flows through subglacial aquifers under a potential gradient less than or equal to that from the ice-air surface slope and aquifer slope. If the water supply exceeds these two sinks, then an ice-contact basal water system develops;
– in the absence of a basal water system, basal velocity is limited to slow basal sliding without cavitation, and is assumed negligible here;
– water produced by basal melting supplies a distributed, high-pressure basal water system (Reference WeertmanWeertman, 1972; Reference AlleyAlley, 1989). With increasing water supply, this system occupies an increasing fraction of the bed, “drowning” some roughness elements, and speeds sliding. “Drowning” an obstacle and increasing the areal fraction of water at the bed requires that the water pressure increase above the local ice-rock contact pressure;
– an increase in water pressure in the basal system decreases the viscosity and yield strength of subglacial sediments, allowing or increasing sediment deformation;
– till generation increases with clast velocity and clast force on substrate (Reference HalletHallet, 1979; assuming that abrasion is a major source of till), and depends on substrate resistance.
The modeling procedure adopted is relatively simple. Surface temperature, snow accumulation rate, bed elevation, and geothermal flux are specified at each grid point along a flow line. Finite-difference integration is begun for specified ice-surface elevation at the divide and proceeds outward. For each distance step, an iterative procedure is used to find the steady surface slope and thus ice thickness. Each surface slope tested for a distance step yields an ice thickness, balance velocity, and basal shear stress. The Reference RobinRobin (1955) steady-temperature model is used to calculate ice temperature and basal temperature. Ice velocity from internal shearing is estimated using a simple, shear-stress-only, depth-averaged model allowing for temperature dependence of creep (after Reference AlleyAlley, 1984).
If the basal temperature equals the pressure-melting point, then the basal melt rate is calculated. Water loss to saturate porosity produced by erosion and to enter subglacial aquifers is subtracted from the melt rate. If excess water remains, it is added to any excess water from upstream to yield the basal flux. This basal flux is used to calculate the average thickness and pressure of the basal water system. The water thickness, basal shear stress and bed roughness (which varies with water thickness and thus with water pressure) yield the sliding velocity, and the water pressure, basal shear stress, and bed properties and thickness yield the velocity from deformation of subglacial sediments.
Generation of till is assumed proportional to the strain rate and effective pressure in till, and till continuity is used to calculate till thickness. The dynamic velocity then is the sum of the sliding, till-deformational, and depth-averaged ice-shearing velocities. The surface slope is varied until the balance and dynamic velocities match. This procedure is repeated down-stream, yielding a steady ice-sheet profile.
Results
A variety of simulations has been run, using constants and boundary conditions appropriate for the West Antarctic flow line through Ice Stream Β (see Table I). The model flow lines show two fundamental patterns, one dominated by ice shearing and the other dominated by basal processes. All model runs find ice shearing to be dominant near the divide, producing the classic convex-up ice sheet with surface slope and basal shear stress increasing along flow to discharge the increasing ice flux from surface accumulation along flow. If the geothermal flux is sufficiently high (> ≈ 50 m W m−2) and the subglacial aquifer is sufficiently impermeable (< 0 (10−5 m2 s−1)), then a transition occurs downstream to reduced shear stress, slow ice shearing and fast basal velocity. Such behavior is shown in Figure I. Downstream of this transition, basal shear stress remains nearly constant at a value somewhat greater than the till yield strength; however, slow increase or decrease in basal shear stress may occur along flow, depending on parameters chosen.
Fig. 1. Ice thickness, hi in km, deforming bed thickness, hb , in m, surface velocity, usfc , in m a−1, and bed velocity, ubed , in m a−1for flow-line simulation using values listed in Table I, no lateral convergence, and program that chooses steeper slope when two possible slopes exist.
Table I. Constants used in Simulations, Unless Otherwise Specified. see Alley (1989), Reference Alley, Blankenship, Rooney, Bentley and GiachiAlley and Others (1989) for Justification of most Choices
I do not model longitudinal stresses in the ice, so the transition from ice shearing to basal velocity can be quite abrupt. This causes the modeled surface deceleration along flow in Figure 1. Upstream of the transition, ice shearing dominates and surface velocity exceeds depth-averaged velocity; however, downstream of the transition surface velocity decreases essentially to the depth-averaged velocity.
(Depth-averaged velocity increases monotonically downstream.)
In many areas both the ice-shearing and the basal-velocity configurations are steady. In these areas, a steep surface slope causes rapid ice shearing equal to the balance velocity. The steep surface slope also causes a large pressure gradient in subglacial aquifers, allowing all of the water produced to be drained through porous aquifer flow; thus, an ice-contact basal water system is not developed and basal sliding and till deformation are slow or absent. Under the same conditions, a gradual surface slope causes ice shearing to be slow compared to the balance velocity but reduces the water flux through subglacial aquifers, allowing a basal water system to form and to lubricate till deformation and sliding. Both steady states are stable.
It is probable that a real glacier would “choose” one of these two states based on historical accident or two-dimensional bed irregularities. During initial growth, a glacier would evolve into one or other of these states. However, a sufficiently large perturbation in the appropriate direction would cause a switch to the other steady state. The implications of these results for ice-stream initiation are interesting and merit further study.
Figure 2 shows results of a sensitivity study based on Figure I. The distance from the ice divide to the grid point farthest downstream with convex-up surface slope (ice-shearing-dominated) is plotted for varying geothermal flux and aquifer capacity. One series of numerical experiments was conducted in which the program chose the steeper steady surface slope at each step, and a second series chose the more gradual steady slope. The horizontal separation between these curves in Figure 2 is the width of the zone (envelope) where two steady states exist. I arbitrarily truncate the solution domain 600 km from the ice divide. In the case of very high aquifer capacity and low heat flow (no basal motion), the ice front is reached at about 510 km. To allow a longer ice sheet, the divide thickness would need to be increased above the 3000 m chosen here. The geothermal fluxes and aquifer capacities used in Figure 2 fall within the range of possible natural conditions.
Fig. 2. Sensitivity of distance to grid point farthest downstream with convex-up surface slope in simulation shown in Figure 1 for runs that choose steeper slope (STEEP) and gentler slope (FLAT), a. Sensitivity to variation in aquifer capacity, K. b. Sensitivity to geothermal flux, G.
Figure 2 assumes a readily erodible (or semi-consolidated) bed. Figure 3 shows sensitivity to bed erodibility, assumed constant along flow, using values for other parameters listed in Table I. The basal shear stress and till thickness are plotted at x = 600 km from the
Fig. 3. Sensitivity of simulations of Figure 1, for runs that choose steeper slope (STEEP) and gentler slope (FLAT), to variations in bed erodibility, Kt . Plotted are values at x = 600 km of the basal shear stress, Tb , and the till thickness, hb. STEEP and FLAT plot sufficiently close together for hb that only one curve is shown.
divide. For these simulations, the initiation of basal motion occurs between about x = 200 km and x = 400 km, well upstream. Deforming-bed thickness proves to be only weakly dependent on the point of initiation of basal motion, causing both cases to plot within the width of the curve shown on Figure 3. Basal shear stress at the downstream end shows a stronger dependence on point of initiation of basal motion for resistant beds, although very erodible beds reduce the shear stress close to the till yield strength regardless of the choice upstream.
Figure 4 shows an attempt to simulate more closely the flow line through Ice Stream B, West Antarctica. Lateral convergence and downstream shallowing of the bed are specified to be similar to observed values (e.g. Reference AlleyShabtaie and others, 1988). Bed erodibility is specified to increase downstream through the convergence zone, equivalent to ice streams occurring on and following highly erodible, soft sediments. The match to the real ice stream is not exact (nor was an exact match attempted), but the similarities are encouraging. (The reader should remember that the ice stream is not in steady state (e.g. Reference Shabtaie, Bentley, Bindschadler and MacAyealShabtaie and others, 1988) and that there is no unequivocal evidence of bed deformation beneath Ice Stream B.)
Fig. 4. Attempt to simulate longitudinal profile of Ice Stream B. a. Model (continuous curve) and observed (x; Reference Shabtaie, Bentley, Bindschadler and MacAyealShabtaie and Bentley, 1988) flow-band width, with width at divide taken as unity, and specified bed erodibility, Kt in m P a−1, b. Resulting basal velocity, ubed , and surface, usfc velocity (continuous curves), for program choosing steeper surface slope whenever possible, and observed surface velocities (x; Reference Shabtaie, Bentley, Bindschadler and MacAyealShabtaie and Bentley, 1988), all in m a−1, c. As in b, but model (continuous curve) and observed (x; Reference Shabtaie, Bentley, Bindschadler and MacAyealShabtaie and Bentley, 1988) surface elevation, zs , and bed elevation, zb , in km, and model bed thickness, hb, in m. Note that zb is specified, zs is a model output, and hb is plotted as thickness about zb.
Discussion and Conclusion
The results here support the finding of Reference Boulton and JonesBoulton and Jones (1979) that the subglacial aquifer system exerts an important control on ice-sheet profile. Given the numerous differences between their model and mine, this result should be considered robust. Efficient basal aquifers extending beyond the ice-sheet margin can lower basal water pressures and reduce basal velocities; inefficient aquifers promote fast basal motion.
Imbalance between heat supplied to the bed (geothermal flux + heats of sliding and deformation) and conducted into the ice is critically important. Ice frozen to its bed or with only slow basal melting lacks lubricating water and fast basal velocities. In this respect, both the high geothermal flux and the low bed elevation of West Antarctica promote
fast basal motion; for given ice thickness, lower bed elevation places the surface at lower, warmer elevation in the atmosphere, reducing heat conduction through the ice and promoting basal melting.
The resistance of the bed to abrasion also proves to control the ice-surface profile, as shown in Figure 3. The values used for bed hardness in Figure 3 were chosen in a somewhat arbitrary manner, but probably span the range of real materials. It is difficult to conceive of sub-till abrasion of a fresh granite being rapid enough to generate a meters-thick till layer, and it is equally difficult to conceive of active ice exerting a significant shear stress on unconsolidated basal sediments under high water pressure without deforming them to a significant depth. Both fresh granite and unconsolidated sediments are possible natural substrates, and the values in Figure 3 were chosen empirically to span this range of behaviors.
The existence of two steady states under a range of conditions is an exciting result readily understandable in terms of glacier physics. It introduces the possibility of a bifurcation, with a small perturbation causing a large switch between stable states.
Only two steady states appeared in the current simulations, but the equations in the appendix suggest the possibility of a third. If the bed slopes upstream, it will tend to trap basal water beneath a sufficiently flat surface slope. Under such circumstances, a greatly thickened water layer might “drown” common roughness elements (if roughness elements ≥O(10 mm) high are scarce or absent; Reference Weertman and BirchfieldWeertman and Birchfieid, 1982) causing the sliding velocity to approach the balance velocity for a basal shear stress reduced close to or below the till yield strength (≈ cohesion) so that bed deformation and ice shearing are slow. This is one proposed explanation for the existence of ice plains with very low surface slope at the mouth of Ice Stream Β (Reference Alley, Blankenship, Bentley and RooneyAlley and others, 1987); it deserves further scrutiny. (The bed assumed here is sufficiently rough to suppress this behavior.)
The model here would be improved by better knowledge of some physical systems (especially bed physics and roughness) and by changes in some boundary conditions (for example, adding elevation-dependent surface temperature and accumulation rate). A sufficiently strong elevation dependence of the surface boundary condition might cause one of the steady states to disappear, for example. Two-dimensional, time-dependent modeling eventually will be needed.
Acknowledgements
I thank D.R. MacAyeal for helpful suggestions. This research was supported in part by the U.S. National Science Foundation, Division of Polar Programs under grant DPP-8716016.
Appendix
The bed model calculates basal temperatures following Robin’s (1955) steady state treatment in which dissipative heating is assumed to occur at the bed. The basal melt rate, m, is caused by the sum of the geothermal flux and the dissipative heating less the conduction into the ice. The water supply, q, then is:
where x is the distance coordinate along flow with x − 0 at the ice divide. This water is lost into porosity generated by erosion at a rate
The maximum flux through subglacial aquifers is:
where ρi is the ice density, ρw is the water density, g is the gravitational acceleration, is the ice-air surface slope, αb is the ice-bed surface slope, and the constant Κ is the aquifer capacity and depends on the aquifer thickness h a and the aquifer hydraulic conductivity Ka. This assumes that the aquifer overlies an aquitard and that the aquifer thickness is small compared to its length so that one-dimensional flow is a good approximation. Typically, the term containing the bed slope in Equation (3b) is small compared to that containing the surface slope, and can be ignored. The model basal water system (at the ice-bed interface) drains water qw given by:
The model approximates the basal water system as a film of average thickness d given by:
where μ = 1.8 × 10−3 Pa s is the viscosity of water (Reference WeertmanWeertman, 1972). The model basal sliding velocity, u s, is:
where Ks is a constant related to the bed roughness divided by the controlling-obstacle height (Reference Weertman and BirchfieldWeertman and Birchfield, 1982) and where the basal shear stress, τ, is:
with hi the ice thickness.
At the ice-bed interface, the model calculates the effective pressure (the difference between the íce-overburden pressure and the water pressure), N, as (Reference AlleyAlley, 1989):
where β is a measure of bed roughness and f is the fraction of the bed area where a thickened water system exists. For a till bed, f is estimated as (Reference AlleyAlley, 1989):
where d is in meters; other relations can be written for beds with different surface roughnesses.
Bed deformation is approximated by (Reference AlleyAlley and others, 1989):
where ub is the horizontal velocity at the top of the deforming till, hb is the deforming thickness, and Kb is a till-softness coefficient. Here τ * is the till yield strength:
with tan ϕ the internal friction and C the cohesion. This treatment of bed deformation assumes no variation in effective pressure with depth in the till and no resulting yield-strength limitation on deforming thickness. It also ignores possible ploughing or discrete shearing. Equation (10) is based on a best-fit result for Ice Stream Β calculated with these assumptions. More – complete formulations should be tested in the future, but are unlikely to change results qualitatively.
The till thickness, hb, is calculated from continuity and a till-gcneration relation. I use:
where Kt is an abrasion softness, ub/hb is the average till strain-rate (and thus measures the velocity of clasts moving over a rigid substrate), and N is the effective pressure (and thus measures the pressure of a moving clast on a rigid substrate). The physical basis for this relation is weak. However, bed erodibility is likely to vary over many orders of magnitude whereas variations in other parameters affecting till generation are likely to be smaller; thus, it may be more important to know whether a bed is unconsolidated or highly consolidated than to know the exact form of the erosion relation.
The depth-averaged ice-deformation velocity, ui , is calculated for an ice sheet without longitudinal stresses following Alley (1984). The depth-averaged dynamic ice velocity, ud , then is:
and is compared to the balance velocity to determine steady state. Values of constants used in model runs are taken from earlier model attempts to match Ice Stream Β itself, but not the entire flowline from the ice divide (Reference AlleyAlley and others, 1989). These are listed in Table I.
