Seasonal acceleration of inland ice via longitudinal coupling to marginal ice

Abstract

We use an ice-flow model to demonstrate how flow variations initiated in the marginal zone of an ice sheet affect flow farther inland through longitudinal (along-flow) coupling. Our findings allow for an alternate interpretation of seasonal accelerations observed near the equilibrium line of the Greenland ice sheet (Zwally and others, 2002). We demonstrate that these observations can be explained by accelerations initiated up to 12 km closer to the margin where the ice is ∼40% thinner, is heavily crevassed, experiences a seasonal doubling of velocity, and where the ablation rate, surface meltwater flux and ice temperature are likely higher. Our modeling and observations suggest that conditions and processes normally found near ice-sheet margins are adequate for explaining the observations of Zwally and others (2002). This and considerations of the likely subglacial hydrology in the marginal zone lead us to suggest that seasonal accelerations may have limited impact on ice-sheet mass balance even in the face of climate warming.

Introduction

A correlation between surface meltwater generation and short-lived (≤1 month), seasonal accelerations in surface velocity of 10–20% was observed at Swiss Camp (Zwally and others, 2002; hereafter ‘Z02’), located near the equilibrium-line altitude (ELA), 40 km inland from the western margin of the Greenland ice sheet (GIS) (Fig. 1a). These observations have led to renewed interest in the process of hydrologically driven ice fracture (hereafter, ‘hydrofracture’) whereby pre-existing, near-surface fractures fill with water from surface streams and lakes and, due to the density difference between ice and water, propagate through the ice, ultimately connecting the supra- and subglacial water systems (e.g. Boon and Sharp, 2003). Theoretical work suggests that this process may apply at the flanks of polar ice sheets, where the ice is thick (>1000 m) and cold (subfreezing) (Weertman, 1973; Van der Veen, 1998, 2007; Alley and others, 2005a), and the process is believed to be of fundamental importance for explaining the observations of Z02 (e.g. Alley and others, 2005b).

Fig. 1. (a) Landsat Enhanced Thematic Mapper (ETM) image from 1 August 2001, showing the location of Swiss Camp (SC), the radio-echo sounding (RES) profile (dotted line), the model flowline (dotted and dashed lines) and the heavily crevassed region downstream from SC (dark band running from lower left to upper right near image center), marked ‘CR’ in (b) (study area shown at inset). Surface elevation contours are from Bamber and others (2001). (b) 2 MHz ice-penetrating radar profile along the dotted curve shown in (a). Data in the upper ∼100 m (heavy black band) are not resolved by the radar system. The bed is indicated by the bright, continuous reflector at an elevation of ∼0 m near the righthand side of the image. The heavily crevassed region discussed in the text is labeled ‘CR’. Detailed RES surveys and analysis suggest that the vertical diffractor ∼1 km downstream from SC represents a vertical conduit that failed to connect with the subglacial drainage system (J.S. Greenbaum and others, unpublished information). The vertical diffractor ∼4 km downstream from SC likely coincides with a ∼600 m instrumented borehole from fieldwork that took place during the 1990s (Thomsen and others, 1991).

By connecting the subglacial hydrology and the ice-sheet surface, hydrofracture allows for the possibility of a fast link between increased surface melting (triggered by increased atmospheric temperature) and increased ice-sheet motion (through increased lubrication at the ice–bed interface). Assuming that hydrofracture can operate near the ELA of a polar ice sheet, it has been proposed as a ‘mechanism for rapid, large-scale, dynamic responses of ice sheets to climate warming’ (Z02). Below, we refer to this as the ‘ELAhydrofracture hypothesis’. Parizek and Alley (2004) investigated the effects of this hypothesis by using the observations of Z02 to link the rate of sliding and the volume of up-glacier surface melt on annual timescales in a predictive model of the GIS. They found that the link leads to a positive feedback: melt-induced sliding on the ice-sheet flank promotes thinning and flattening there, enlarges the area experiencing surface melting and sliding, and allows the ELA and the onset of sliding to propagate increasingly far inland over time. Simply put, a warming climate exposes a greater fraction of the ice sheet to surface-melt-induced sliding over time. By including this mechanism, Parizek and Alley (2004) predict a more negative mass balance for the GIS, sooner, than if the mechanism is neglected.

Here, we propose a modified form of the ELA-hydrofracture hypothesis to explain the observations of Z02. Assume that at some distance, D, downstream from Swiss Camp (SC) surface melt reaches and lubricates the ice-sheet bed, resulting in seasonal acceleration by a factor A. The acceleration at and downstream from D also causes an instantaneous acceleration at SC through the effects of longitudinal (along-flow) coupling. If coupling is weak for a reasonable value of A, such that D must be small in order to match the observations of Z02, surface melt must be accessing the bed and affecting sliding locally (at SC). If coupling is strong and D is significant for a reasonable value of A, a local change in lubrication at SC is not required to explain the observations of Z02. Below, we use an appropriate ice-flow model and recent field observations to explore this alternate hypothesis and its broader implications.

Flow Modeling

We use a two-dimensional flowline model to explore how temporal changes in velocity downstream from SC could affect velocities at SC. The model, which solves the full two-dimensional stress equilibrium equations (i.e. plane strain) using the finite-volume method, is discussed in detail in Price and others (2007) and is applied similarly in Price and others (2008).

In two-dimensional Cartesian coordinates, conservation of momentum for a viscous fluid in a low Reynolds number flow is expressed by

(1)

where x and z are the along-flow and vertical coordinates, respectively, and repeat indices imply summation. The first term on the lefthand side of Equation (1) is the body force, the product of ice density, ρ, and the acceleration due to gravity in the i direction, g_i . The second term is the stress divergence with the full stress tensor, σ_ij , given by the deviatoric stress, τ_ij , minus the pressure, P:

(2)

where δ_ij is the Krönecker delta (or identity matrix). The constitutive relation linking deviatoric stress and strain rate is given by

(3)

where is the strain-rate tensor,

(4)

and η is the effective viscosity,

(5)

In Equation (4), u_i represents the components of the velocity vector u (in the x direction, i) or w (in the z direction, j). In Equation (5), B(T) is the temperature-dependent rate factor described by an Arrhenius relation (e.g. Paterson, 1994, p. 86), n is the power-law exponent (taken equal to 3) and is the effective strain rate given by

(6)

Ice is assumed to be incompressible such that

(7)

The finite-volume grid, shown as thin gray lines in Figure 2, has a horizontal grid spacing of ∼850 m. The ice thickness is discretized by 20 finite volumes, using ‘sigma coordinates’, with a variable spacing so resolution increases near the ice–bed interface. Further details of the model and solution methods are given in Price and others (2007).

Fig. 2. (a) Mean-annual horizontal velocity field for isothermal model. Black dotted curves show the actual surface and bed topography. (b) Instantaneous horizontal velocity field after the basal resistance has been reduced from D km downstream of SC to the margin. In Figure 3a, this particular model would plot at D = 15 km, A = 1.25, and would be associated with a ∼2% acceleration at SC. Thin gray curves represent the computational grid.

Figure 2 shows a portion of the model domain near SC, which is based on surface and bed elevation profiles obtained through ground-based global positioning system (GPS) and radio-echo sounding (RES) surveys (data shown in Fig. 1b). The profile is oriented along flow and covers the region within ±20 km of SC. Surface elevations along the profile have been fitted to a smooth parabola, and bed elevations have been smoothed with a high-order polynomial fit. From ±20 km downstream of SC to the ice-sheet margin we have extrapolated a smooth bed slope over 20 km so the surface and bed profiles meet in a small ice cliff at the ice-sheet margin, 40 km downstream from SC. The full model domain extends several hundred kilometers upstream from the study area shown in Figure 2, where we apply a zero-flux boundary condition at an ice divide. We note that a limited number of experiments in which we account for the actual bedrock topography (shown as a dotted line in Fig. 2a) confirm that the primary model results discussed below (Fig. 3) are not significantly affected by the use of smoothed bedrock topography. Here, we use smoothed bed topography for computational efficiency.

Fig. 3. Contours of percent increase in surface velocity at SC over initial value for (a) the isothermal model (low-sliding fraction) and (b) the polythermal model (high-sliding fraction) as a function of the parameters D (horizontal axis) and A (vertical axis) discussed in the text. Crosses indicate the values of D and A used to construct contours. Dotted lines denote the range of D for which A = 2 satisfies the observations of Z02.

Thermomechanical-flow modeling studies in the region (Funk and others, 1994; Lüthi and others, 2002; Wang and others, 2002) indicate that the ice–bed interface in the region is at the pressure-melting point. This indicates that some fraction of the observed surface velocity in the vicinity of SC results from motion at, or very near to, the glacier bed, either as a result of sliding or through enhanced deformation in temperate, soft basal ice (Lüthi and others, 2002). In our model, basal motion occurs through deformation within a several-meters-thick layer, beneath the ice that is discretized by several layers of gridcells (Price and others, 2008). The sliding speed is defined as the horizontal velocity at the top of this basal-fluid layer, for which the rate of deformation is given by

(8)

In Equation (8), the effective stress, τ_e, is analogous to the effective strain rate defined in Equation (6). We specify the rate factor, C, and the power-law exponent, p, in order to control the fraction of the surface velocity that results from motion near the bed. In all experiments discussed below, we specify p = 1, in which case our sliding law is analogous to viscous basal-slip laws commonly employed in other flow-modeling studies (e.g. MacAyeal, 1989). The value of C is initially held steady and uniform in order to match observed surface velocities along the model flowline (C is taken as 2.5 × 10⁻⁵ and 6.7 × 10⁻⁵ Pa⁻¹ a⁻¹,for the isothermal and polythermal models, respectively, equivalent to viscosities, C^{− 1}, of 4.0 × 10⁴ and 1.5 × 10⁴ Pa a, respectively). For the models discussed below, the initial surface velocity in the region of SC is 35–40 cm d⁻¹, similar to the mean-annual value reported by Z02. At a location ∼15 km downstream from SC, model surface velocities are ∼25 cm d⁻¹, equal to measurements of the mean-annual surface velocity at this location (Rumrill and others, 2006; Neumann and others, 2007). We take the initial, horizontal velocity field to be representative of the mean-annual horizontal velocity in the region near SC. (A similar ‘tuning’ of the sliding velocity for the case of p > 1 would require smaller C, a thinner basal-fluid layer, or some combination of the two.)

Over the short timescales considered here (<1 year), changes in ice thickness due to non-steady flow have a minimal effect on the domain geometry and on the modeled stress and velocity fields. This is confirmed by time-dependent model runs (discussed further below) for which we apply a simplified but reasonable pattern of surface mass balance (i.e. accumulation increasing with distance upstream from the modern-day ELA, near SC, and ablation increasing with distance downstream from the ELA). Thus, we simplify the modeling by comparing two instantaneous velocity fields, one representing the unperturbed, mean-annual flow field (tuned to the observed mean-annual velocities in the region, as discussed above) and one representing the perturbed, seasonally accelerated flow field. There is no requirement that the initial flow field be in steady state and, by this approach, we avoid the detrimental effects of a standard model ‘spin-up’, which would lead to significant changes in the initially prescribed domain geometry. Uncertainties in the initial temperature field (and thus B(T) in Equation (5)) are not treated through any explicit tuning and are discussed further below.

For a portion of the model domain extending from D km downstream of SC to the ice-sheet margin, we assume that surface-generated meltwater can reach the glacier bed, causing an increase in lubrication and accelerated basal motion over that region relative to the mean-annual motion (Fig. 2) (the implementation of this process is discussed in more detail below). In the present context, a logical measure of this acceleration is the fractional increase in surface velocity that takes place at D, the location of the transition from relatively slow to relatively fast basal motion (i.e. the transition from relatively less to relatively more basal lubrication). We define this fractional increase as A. The ELA-hydrofracture hypothesis implies that accelerations originate local to SC, in which case D ≈ 0 km and A = 1.1–1.2; all of the 10–20% acceleration observed at SC is attributed to a change in basal lubrication and basal motion at, or very close to, SC.

We use our flow model to explore other reasonable combinations of D and A that result in accelerations of 10–20% at SC. We limit our exploration of A to a range of 1–3, based on observations of seasonal acceleration in northeast Greenland (Mohr and others, 1998), similar observations on other high-latitude, polythermal glaciers (Iken, 1974; Bingham and others, 2003; Copland and others, 2003) and recent observations downstream from SC (Rumrill and others, 2006; Neumann and others, 2007). We limit our exploration of D to between 2 and 30 km; a value of D ≤ 2 km is essentially ‘local’ to SC, and model results indicate that for D ≥ 30 km the range of A noted above has little or no influence on the velocity at SC.

To obtain a particular value of A, which is a model output, we increase the amount of basal motion by increasing the rate factor, C (Equation (8)); at and downstream from D, Cis increased as a step function by a factor, f, of its initial value, C ₀. In reality, a change in basal motion resulting from a change in meltwater input will be time-dependent so C = f(t)C ₀, which generates A = A(t). We used time-dependent model runs with time-steps ∼1 month to test the model sensitivity to different functional forms of f(t) (and the resulting A(t)) over the period of a year, and found that at SC the maximum increase in surface velocity as a function of time is similar, regardless of the shape of f(t); the largest increase in surface velocity at SC occurs nearly coincident with the greatest increase in f(t) (and thus the largest value of A(t)). These observations confirm that thickness changes (i.e. thinning) resulting from the increase in basal motion have a minor effect on the velocity field over the period of a melt season. Thus, we simplify the modeling by removing time dependence; over the domain of interest, the maximum difference between the mean-annual horizontal velocity field and that resulting from the perturbation C = f(t)C ₀ is well represented by the instantaneous response to the perturbation C = f _max C ₀.

Results

Model results are summarized in Figure 3. Because we have incomplete information on (1) the temperature of the ice, and therefore the value of the flow-law rate factor, and (2) the fraction of the surface velocity arising from motion at the glacier bed, we explore the parameter space for two end-member sets of models. First, we consider isothermal models (temperate ice), with a uniform flow-law rate factor, for which the fraction of surface velocity, u _sfc, arising from basal motion, u _bed, is small in the region of SC (u _bed/u _sfc ∼ 0.3) (Fig. 3a). Second, we consider polythermal models for which u _bed/u _sfc is large in the region of SC (u _bed/u _sfc ∼ 0.8) (Fig. 3b). For the polythermal models, we use the ‘ice-sheet’ temperature profile modeled by Funk and others (1994; similar to the profile at SC calculated by Wang and others, 2002) to describe a depth-varying but overall stiffer rate factor within the ice column. (We make this simplification rather than conducting a detailed calculation of the temperature field due to (1) uncertainties that would be introduced into the temperature field by assuming it is in steady state with modern-day variables or, in the case of a time-dependent calculation, as a result of uncertain climate and ice dynamics histories, and (2) the lack of a measured temperature profile in the region to use as a modeling target.) Because we expect that cold, stiff ice that is predominantly sliding will result in more effective longitudinal coupling than warm, soft ice that is primarily deforming (Kamb and Echelmeyer, 1986), the results for these two sets of models place additional bounds on the parameter space.

Discussion

The flow model is able to reproduce the seasonal velocity accelerations observed at SC for multiple combinations of D and A. Values consistent with the ELA-hydrofracture hypothesis plot in the lower-lefthand corners of Figure 3a and b. For D > 0, an increasingly large value of A is required. Rumrill and others (2006) and Neumann and others (2007) report a seasonal acceleration of twice the mean-annual velocity (i.e. A = 2) at a location ∼15 km downstream from SC (GPS1 and 2; Fig. 1a). Figure 3 indicates that, for A = 2, the observations at SC can be explained by increased basal motion originating from ∼5–12 km downstream from SC. Based on borehole measurements and thermomechanical modeling studies in the region (Funk and others, 1994; Lüthi and others, 2002; Wang and others, 2002), the polythermal model results (Fig. 3b) are more likely representative of the longitudinal coupling length near SC. Hence, the modeling conducted here supports the contention that a ∼2 × seasonal acceleration initiated ∼12 km downstream from SC can result in a 10–20% acceleration at SC.

Other factors support the contention that large (≥2×), surface-meltwater induced, seasonal accelerations may be preferentially initiated in the region 12–15 km downstream from SC, rather than at SC itself. The surface in this region is heavily crevassed, as observed in the field by G.C. and T.N. and as indicated by multiple, chaotic diffractors observed in RES profiles (‘CR’; Fig. 1b). Crevassing is likely the result of extending flow* over a bedrock bump (x = 15 km; Fig. 1b). Heavily crevassed, relatively thinner ice (Fig. 1b) at relatively lower elevation (Fig. 1a) is likely to allow surface meltwater relatively easier access to the ice–bed interface.

We hypothesize that the area from CR to the ice-sheet margin is the source for significant seasonal accelerations in the study area. Further, we hypothesize that the basal hydrology in this region is similar to that beneath an alpine or High Arctic valley glacier. That is, we suppose that a distributed basal drainage system, with limited capacity, exists near CR during the winter and early in the melt season. As the melt season progresses, and as the flux of meltwater from the surface increases, the capacity of the drainage system also increases by becoming channelized (e.g. Bingham and others, 2003, 2005); the subglacial water system serves as a buffer against continually increasing velocity with increasing surface melt. This is analogous to behavior observed on John Evans Glacier, a High Arctic polythermal valley glacier on Ellesmere Island, Canada. Bingham and others (2003) observed that spring accelerations there were smaller and of shorter duration during the relatively warm summer of 2000 than during the following, relatively cooler summer. Dye-tracing experiments explained the paradoxical result; relatively more surface melt in 2000 led to a faster switch to a channelized basal water system, and thus to relatively lower basal water pressures and sliding velocities. Truffer and others (2005) discuss a similar example based on observations from several Alaskan glaciers.

The differences between the ELA-hydrofracture hypothesis and the hypotheses favored here have implications for the large-scale response of the GIS to future climate warming. In our hypothesis, significant seasonal accelerations are initiated at highly crevassed regions where there is copious meltwater. In the ELA-hydrofracture hypothesis, significant seasonal accelerations can initiate anywhere there are preexisting fractures at the surface and where surface meltwater can collect. Our hypothesis suggests that increased sliding as a function of increased surface melt will be limited by evolution of the subglacial water system; increases in ice flux may be related to increases in surface melt, but will also depend on the type and evolution of the subglacial water system (e.g. modeling of Arnold and Sharp, 2002). To date, the ELA-hydrofracture hypothesis has generally been interpreted such that the effect of seasonal accelerations on annual ice flux (and thus annual mass balance) parallels the quantity of annual surface melt; increases in ice flux are taken as proportional to increases in surface melt (e.g. modeling of Parizek and Alley, 2004). Our hypothesis implies that the effects of meltwater-induced acceleration on inland ice, through longitudinal coupling, will be relatively smaller than implied by the ELA-hydrofracture hypothesis; for a given longitudinal-stress gradient at CR (versus at SC) the force ‘pulling’ on upstream ice is initiated half as far inland and is ∼2 × less (due to the difference in ice thickness). Finally, we note that acceleration on the ice-sheet flank implies thinning, which will cause a given elevation to retreat inland over time. If acceleration and thinning occur nearer to the margin, where the slope is relatively steeper, than nearer to the ELA, where the slope is shallower, the rate of inland retreat will be relatively slower. Thus, our hypothesis implies a smaller rate of ablation-zone widening over time than does the ELAhydrofracture hypothesis. Overall, our hypothesis suggests that the effects of meltwater-induced acceleration in a warming climate will have a smaller impact on the mass balance of the GIS than suggested by the ELA-hydrofracture hypothesis.

Conclusions

We find that seasonal accelerations observed near the ELA of the GIS (Z02) can be explained by larger accelerations, initiated closer to the ice-sheet margin, that affect velocities upstream through longitudinal coupling. While our favored mechanism for the accelerations (seasonal lubrication of the ice–bed interface via surface meltwater) is identical to that proposed previously, the non-local origin of the accelerations, implied here, is important. We propose that accelerations are initiated in highly crevassed, relatively thinner ice, nearer to the ice-sheet margin, where large seasonal accelerations are known to occur. If the glacier hydrology in this region is similar to that beneath other highly crevassed glaciers, the magnitude and duration of seasonal accelerations, their inland extent and their effect on ice-sheet mass balance may be limited, even in the face of climate warming.

It is unlikely that the observations of Z02 represent an isolated phenomenon. Tsai and Ekstrom (2007) report on seismic observations that suggest widespread, seasonal changes in sliding velocity along the flanks of the GIS and note that the seasonality and increasing frequency of events imply a link to temperature or temperature-related variables. While surface meltwater is the most obvious link between changes in atmospheric temperature and changes in basal sliding, it is less obvious how this link might affect large-scale ice-sheet dynamics. Z02 conclude that their observations indicate ‘a mechanism for rapid response of the ice sheets to climate change’. Recent reviews in glaciology, scientific editorials and cross-disciplinary studies (see Appendices A–C) suggest that this hypothesis has been largely accepted by the scientific community. As a consequence, similar statements have been communicated in the popular press (see Appendix D) with the appearance of scientific consensus. Based on the work presented here, we suggest that further study and a more nuanced interpretation are necessary before these particular conclusions of Z02 are fully warranted.