Considering thermal-viscous collapse of the Greenland ice sheet.

We explore potential changes in Greenland ice sheet form and flow associated with increasing ice temperatures and relaxing effective ice viscosities. We define "thermal-viscous collapse" as a transition from the polythermal ice sheet temperature distribution characteristic of the Holocene to temperate ice at the pressure melting point and associated lower viscosities. The conceptual model of thermal-viscous collapse we present is dependent on: (1) sufficient energy available in future meltwater runoff, (2) routing of meltwater to the bed of the ice sheet interior, and (3) efficient energy transfer from meltwater to the ice. Although we do not attempt to constrain the probability of thermal-viscous collapse, it appears thermodynamically plausible to warm the deepest 15% of the ice sheet, where the majority of deformational shear occurs, to the pressure melting point within four centuries. First-order numerical modeling of an end-member scenario, in which prescribed ice temperatures are warmed at an imposed rate of 0.05 K/a, infers a decrease in ice sheet volume of 5 ± 2% within five centuries of initiating collapse. This is equivalent to a cumulative sea-level rise contribution of 33 ± 18 cm. The vast majority of the sea-level rise contribution associated with thermal-viscous collapse, however, would likely be realized over subsequent millennia.

Introduction

The transfer of terrestrial ice into the sea accounted for over half the 3.3 ± 0.4 mm/a sea‐level rise observed during the 1993–2007 period [Cazenave and Llovel, 2010]. While the Greenland ice sheet is believed to have been in near balance under the relatively cool pre‐1990 climate, recent climate warming has resulted in an increasingly negative ice sheet mass balance [Rignot et al., 2008]. The Greenland ice sheet is now the single largest cryospheric source of sea‐level rise. Between 2005 and 2010, Greenland lost 263 ± 30 Gt/a of ice, equivalent to a sea‐level rise contribution of 0.73 ± 0.08 mm/a [Shepherd et al., 2012]. This ice loss stems from both surface mass balance processes, which increase meltwater runoff, and ice dynamic processes, which increase iceberg discharge. Although these two mass loss components were approximately equal in magnitude over the 2000–2008 period [van den Broeke et al., 2009], the surface mass balance component of Greenland mass loss appears to have accelerated faster than the ice dynamic component since ca. 2005 [Enderlin et al., 2014].

Under the representative concentration pathway (RCP) 2.6–8.5 scenario range, global mean sea‐level rise by year 2100 will likely be between 0.40 and 0.63 m, with the majority of sea‐level rise due to the thermal expansion of sea water [Church et al., 2013]. This RCP range, corresponding to +2.6 and +8.5 W/m2 of radiation forcing in year 2100 relative to pre‐industrial times, spans the complete range considered by the Intergovernmental Panel on Climate Change Assessment Report 5 (IPCC AR5). Under this range of climate forcing, Greenland is projected to contribute between 0.06 and 0.12 m of this sea‐level rise [Church et al., 2013]. Due to both committed climate change and anticipated further climate forcing, however, the post‐2100 sea‐level rise will be more substantial than pre‐2100 sea‐level rise [Collins et al., 2013]. Under a medium greenhouse gas forcing scenario of 500–700 ppm CO2, which corresponds approximately to RCP 4.5, the IPCC projects a total Greenland ice sheet sea‐level contribution of up to 0.90 m by the year 2500 [Church et al., 2013]. While the ice dynamic contribution to Greenland's current mass loss is presently less than the surface mass balance contribution, there is substantial uncertainty associated with the evolution of ice dynamics over longer timescales.

The uncertainty associated with projecting ice dynamics is encapsulated by an end‐member kinematic constraint scenario that postulates 0.80 m as the “plausible” upper limit of the year 2100 sea‐level rise contribution due to accelerated ice dynamics throughout the entire cryosphere [Pfeffer et al., 2008]. The uncertainty associated with projecting ice dynamics ultimately stems from the nonlinear ice physics governing glacier form and flow. Cryo‐hydrologic warming, whereby glacier ice is warmed by latent heat released by refreezing water, has been postulated as one such mechanism for the rapid thermal response of ice sheets [Phillips et al., 2010]. Effective ice viscosity, and thus ice deformation and flow, is nonlinearly dependent on ice temperature, which means relatively small changes in ice temperature can result in relatively large changes in ice deformation (Figure 1). Cryo‐hydrologic warming has been recognized by IPCC AR5 as one of several processes that suggest changes in ice sheet form and flow can happen more rapidly than was previously recognized in IPCC AR4 [Vaughan et al., 2013]. The complex processes and feedbacks influencing ice sheet behavior yield different effects, and the overall net response is some combination of many driving factors. For example, the progressive replacement of relatively low‐viscosity Wisconsin ice by relatively high‐viscosity Holocene ice contributes to subtle transient ice sheet slowing and thickening [Reeh, 1985; Colgan et al., 2015]. Increases in ice temperature, however, act to decrease effective ice viscosity, resulting in transient ice sheet speedup and thinning.

Figure 1

(a) The nonlinear dependence of effective ice viscosity on ice temperature means that, for a given stress, strain is greater by an order of magnitude in ice at 0°C than in ice at −10°C. (b) Horizontal deformational ice velocity versus normalized ice thickness in 0 and −10°C isothermal ice under common stress. In the polythermal scenario, the deepest 15% of the −10°C ice column is warmed to 0°C.

In the same spirit as exploring the sea‐level rise associated with a potential marine collapse of the West Antarctic ice sheet [Bamber et al., 2009], we investigate the sea‐level rise associated with a potential thermal‐viscous collapse of the Greenland ice sheet over a multicenturial timescale. We define “thermal‐viscous collapse” as a transition from the polythermal ice sheet temperature distribution characteristic of the Holocene, to an isothermal ice sheet at the pressure melting point temperature throughout. Given their inherently hypothetical nature, end‐member assessments, such as marine instability induced collapse of the West Antarctic ice sheet or thermal‐viscous induced collapse of the Greenland ice sheet, can be difficult to interpret. We are therefore compelled to clearly acknowledge that we do not attempt to constrain the probability of thermal‐viscous collapse of the Greenland ice sheet, but rather merely demonstrate that initiating thermal‐viscous collapse appears plausible within four centuries, and assess the sea‐level rise contribution associated with this event. The exercise of considering end‐member cases such as this can provide a useful perspective for other analyses involving the thermal response of ice sheets.

Here, we explore the scenario of thermal‐viscous collapse of the Greenland ice sheet using a two‐pronged approach. First, we provide a conceptual model of three key physical processes required for thermal‐viscous collapse, and an assessment of their plausibility within the next few centuries (Section 2). Second, we employ numerical ice flow modeling, couched in a Monte Carlo approach, to provide a first‐order estimate of the potential sea‐level rise contribution associated with thermal‐viscous collapse over the five centuries following initiation of collapse (Section 3).

Conceptual Model

We consider three key processes required for a partial thermal‐viscous collapse: (1) sufficient energy available in projected ice sheet meltwater runoff (Section 2.1), (2) routing of meltwater to the bed of the ice sheet interior (Section 2.2), and (3) efficiency of energy transfer from meltwater to the ice via cryo‐hydrologic warming (Section 2.3; Figure 2). We define “partial thermal‐viscous collapse” as warming the deepest 15% of the ice sheet to the pressure melting point temperature in four centuries. Recalling that the majority of all shear occurs in the deepest 15% of an isothermal ice column, potentially warming this deep layer of the polythermal ice sheet to its pressure melting point temperature would achieve the majority of the enhanced ice deformation and flow associated with warming the entire ice sheet volume to pressure melting point temperature. We refer to warming the entire ice sheet volume to the pressure melting point temperature as “complete thermal‐viscous collapse.” As described in the following three subsections, we find that initiating a partial thermal‐viscous collapse, by warming the deepest 15% of the ice sheet to pressure melting point temperature, is plausible in four centuries.

Figure 2

Three key elements of a partial thermal‐viscous collapse: (1) Sufficient energy available in the projected Greenland meltwater runoff anomaly, (2) Routing of a fraction of meltwater to the bed of the ice sheet interior, and (3) Efficient energy transfer from meltwater to ice via cryo‐hydrologic warming. This cross‐sectional profile reflects mean observed Greenland ice surface and bedrock elevations between 74.1 and 76.4°N [Bamber et al., 2013]. Dashed lines illustrate stylized marine and land glacier termini.

Energy in Projected Meltwater Runoff

Initiating a complete thermal‐viscous collapse in four centuries, warming the entire ice sheet volume to pressure melting point temperature, would require an average energy delivery of 2.71 × 1020 J/a. This rate of energy delivery is equivalent to virtually all latent heat energy available in the projected year 2100 runoff anomaly (detailed below). Thus, initiating a complete thermal‐viscous collapse does not appear feasible in 400 years. It appears plausible, however, to meet the energetic requirements of a partial thermal‐viscous collapse in four centuries, whereby only the deepest 15% of the ice sheet is warmed to pressure melting point temperature. The regional climate model MAR projects a year 2100 meltwater runoff anomaly between 300 and 1400 Gt/a, depending on forcing scenario [Fettweis et al., 2013]. A meltwater runoff anomaly of 700 Gt/a, in the lower end of this range, is equivalent to a latent energy anomaly of 2.34 × 1020 J/a. The contemporary volume‐averaged temperature of the deepest 15% of the ice sheet (262.6 K) is substantially warmer than the volume‐averaged temperature of the entire ice sheet (252.0 K; Price et al., 2011]. Owing to melting point depression with depth, the volume‐averaged, pressure melting point temperature of the deepest 15% of the ice sheet is 271.4 K. Warming the deepest 15% of the ice sheet ∼8.8 K from contemporary temperatures to the pressure melting point requires 7.14 × 1021 J, or 1.79 × 1019 J/a for four centuries.

We now place the latent energy cost of a partial thermal‐viscous collapse in the context of the latent energy available in the projected meltwater runoff anomaly. As stated above, a year 2100 meltwater runoff anomaly of 700 Gt/a is equivalent to a latent energy anomaly of 2.34 × 1020 J/a, which is more than an order of magnitude greater than the energy delivery rate required for partial thermal‐viscous collapse over four centuries. While the majority of supraglacial runoff leaves the ice sheet, attempts to close the water budget of a well‐observed ice sheet watershed near Kangerlussuaq suggest that up to 50% of calculated supraglacial runoff does not leave the ice sheet during the summer melt season, but instead enters longer‐term englacial or subglacial storage systems [Rennermalm et al., 2013]. Glacio‐hydrological modeling near Ilulissat infers a mean glacier water residence time of ∼2 years between surface melt generation and discharge from the ice sheet [Colgan et al., 2011a]. The presence of refrozen ice, within relict crevasses and moulins in the ice sheet ablation zone, attests to the bulk refreezing of englacially and subglacially stored liquid water [e.g., McGrath et al., 2011]. Employing the range of year 2100 runoff anomalies simulated by MAR [Fettweis et al., 2013] and a range of englacial and subglacial meltwater retention fractions inferred by closure of the basin‐scale ice sheet water budget [Rennermalm et al., 2013], there is a wide range of plausible future meltwater anomalies and liquid water storage fractions that potentially satisfy the energy burden associated with warming basal the deepest 15% of the ice sheet to the pressure melting point temperature within four centuries (Figure 3). For example, englacial and subglacial retention of 8% of a 700 Gt/a supraglacial meltwater runoff anomaly represents 1.87 × 1019 J/a of latent energy, which is sufficient for warming the deepest 15% of the ice sheet to pressure melting point temperature in four centuries.

Figure 3

Cumulative latent energy anomaly (in Joules) within the ice sheet after 400 years, estimated by temporally integrating the range of year 2100 meltwater runoff anomalies simulated by MAR [Fettweis et al., 2013] and a range of fractional englacial and subglacial runoff retention informed by basin‐scale water budget closure [Rennermalm et al., 2013]. White contour denotes the 7.14 × 1021 J burden associated with warming the deepest 15% of the ice sheet to pressure melting point temperatures. Asterisk denotes the 700 Gt/a meltwater anomaly and 8% retention fraction discussed in the text.

While there is tremendous uncertainty associated with the role of annual englacial refreezing in basin‐scale water budgets [Flowers and Clarke, 2002; McGrath et al., 2011], invoking a small fraction of the catchment scale englacial and subglacial retention recently inferred in the ice sheet [e.g., Rennermalm et al., 2013] yields sufficient energy to initiate a partial thermal‐viscous collapse under year 2100 meltwater anomaly. As a consequence of both committed climate change and anticipated further climate forcing, however, the year 2100 meltwater anomaly may be a gross underestimate of the annual ice sheet meltwater anomaly characteristic of the subsequent three centuries [Collins et al., 2013]. In addition, the projected energy anomaly described above ignores the nontrivial sensible heat contained in meltwater runoff, which can reach temperatures of 0.1°C in the supraglacial environment and gain a further 0.2°C for each 100 m descent in elevation via the release of potential gravitational energy [Isenko et al., 2005]. Thus, it appears plausible that refreezing a small portion of the anticipated supraglacial meltwater runoff anomaly may be sufficient to overcome the energy burden associated with warming the deepest portion of the ice sheet, where the majority of shear occurs, to pressure melting point temperature in four centuries.

Routing Meltwater to the Interior Ice Sheet Bed

Observations and theory suggest that there is an active subglacial hydrological network in the ice sheet ablation area, in which liquid water persists throughout the winter [Catania and Neumann, 2010; Colgan et al., 2011a; Chandler et al., 2013]. Supraglacial streams over an extensive section of the ice sheet ablation and lower accumulation (“wet snow zone”) areas near Kangerlussuaq exclusively terminate in moulins, indicating that meltwater runoff is primarily routed through or under the ice, rather than across the surface, to the ice sheet margin [Smith et al., 2015]. In the mid (“percolation zone”) to upper (“dry snow zone”) accumulation areas, radar observations infer the widespread presence of subglacial water [Oswald and Gogineni, 2012] and refrozen basal ice masses [Bell et al., 2014], as well as englacial water within the perennial firn aquifer [Forster et al., 2014]. There is only a nascent understanding of the role of subglacial lakes in accumulation area hydrology. The handful of confirmed lakes are generally much smaller than their Antarctic counterparts (<10 km2), but have been postulated to be part of an open hydrological network, capable of moving water around the ice/bed interface [Palmer et al., 2013; Willis et al., 2015]. Notably, the Palmer et al. [2013] lakes were discovered near an ice flow divide in Northwest Greenland, where the ice sheet is commonly assumed to be frozen to the bed. Using first‐order gradients in hydraulic potential, Livingstone et al. [2013] postulate that hundreds of relatively small subglacial lakes may cover ∼1.2% of the ice sheet bed.

While Shreve [1972] is often cited as providing theoretical support for the equipotential gradients governing englacial and subglacial water flow closely following surface slope, in fact Shreve [1972] states that equipotential gradients are “∼11 times more sensitive to ice surface slopes than bedrock slopes,” or ultimately governed by ice thickness gradients. Acknowledging that equipotential theory assumes that subglacial water pressure is equivalent to ice overburden pressure, surface slope should only determine flow direction when bedrock slope is substantially less than the surface slope [Lewis and Smith, 2009]. In the virtually flat interior of the Greenland ice sheet, bedrock slope often exceeds 11 times surface slope (Figure 4). In the best‐surveyed regions of the ice sheet, with the least uncertainty in bed topography, equipotential theory suggests widespread areas where englacial and subglacial water flow is primarily governed by bedrock topography. Given that the interior bed of the ice sheet is depressed below sea level, meltwater introduced to the subglacial network sufficiently far inland should theoretically exhibit a tendency to flow to bedrock sink points in the ice sheet interior. Such reverse drainage and collection behavior has been modeled for smaller ice caps [Flowers et al., 2005]. Recent observations of discrete refrozen basal ice masses confirm that subglacial water can be routed inland toward relatively cold sink points in the ice sheet interior where it refreezes [Bell et al., 2014]. Although considerable uncertainty exists in the precise origin of this refrozen water (i.e., whether surface meltwater, basal meltwater, or ocean water), these observations suggest that nontrivial reverse drainage and interior refreezing do occur under contemporary climate.

Figure 4

(a) Bedrock elevation beneath the Greenland ice sheet (colorbar saturates at extreme values). (b) Regions of the Greenland ice sheet where bed slope is >11 times surface slope (red). (c) Uncertainty in bedrock elevation [Bamber et al., 2013]. Inset box highlights the relatively well‐surveyed Jakobshavn basin.

By increasing meltwater runoff, climate change not only enhances the piezometric potential for reverse drainage but also the delivery of meltwater from the ice sheet surface to its bed. Supraglacial lake drainage events, which can hydrofracture through thick (>1 km) cold ice and introduce large meltwater volumes to the subglacial network [Das et al., 2008], have been observed to occur more frequently and at higher elevations on the ice sheet during relatively intense melt years [Liang et al., 2012]. Similarly, subglacial lake drainage beneath ∼540 m of ice, within the accumulation area of a local ice cap, created a heavily crevassed supraglacial basin that subsequently provided efficient routing of meltwater to the subglacial environment [Willis et al., 2015]. While first‐order (strain rate threshold) models predict that crevasses are unlikely to form above 1600 m [Poinar et al., 2015], crevasse hydrofracture driven by supraglacial lake drainage has been observed to occur above 1800 m elevation under present‐day climate [Doyle et al., 2014]. The recent appearance of a summer velocity increase at 1840 m elevation, in the South Greenland accumulation area, has been interpreted as the recent upstream, or inland, propagation of an active subglacial hydrology network that enhances basal sliding on a seasonal basis [Doyle et al., 2014; van de Wal et al., 2015]. The absence of local crevasses, moulins or lake drainage events at this observation site (“S10”), suggests that subglacial water is being routed over larger distances than the longitudinal coupling length characteristic of such an inland location [Price et al., 2008; Colgan et al., 2012]. Combined with the present understanding of ice sheet englacial and subglacial hydrology and supraglacial lake hydrofracture mechanics, the inference of a transient expansion of active subglacial hydrology into the ice sheet interior suggests it is plausible that the subglacial network may migrate sufficiently far inland that bedrock topography induces subglacial flow toward interior sink points [Flowers et al., 2005; Bell et al., 2014], providing latent heat to the ice sheet interior.

Efficient Heat Transfer From Meltwater to Ice

Cryo‐hydrologic warming of glacier ice by the latent heat released by refreezing water has been both qualitatively [Bader and Small, 1955] and quantitatively [Jarvis and Clarke, 1974] documented. The term “cryo‐hydrologic warming” has been adopted by the IPCC to refer to material warming by latent heat release in glaciers [Vaughan et al., 2013]. In addition to advection, conduction, and production, cryo‐hydrologic warming may be conceptualized as a fourth term in the glacier thermal equation. While transient cryo‐hydrologic warming has been implicated as a potentially important mechanism of recent Greenland outlet glacier acceleration in some cases [van der Veen et al., 2011; Vaughan et al., 2013], thermo‐mechanical ice flow modeling with a four‐term heat equation that includes cryo‐hydrologic warming more faithfully reproduces historical West Greenland ice temperature and velocity profiles than the conventional three‐term heat equation [Phillips et al., 2013]. The contemporary velocity and temperature distributions of the ice sheet can therefore be interpreted to reflect features associated with latent heat transfer from water to ice prior to the onset of recent climate change [Harrington et al., 2015; Lüthi et al., 2015].

A process‐level understanding of cryo‐hydrologic warming remains somewhat elusive. Ultimately, it is the tremendous difference between the latent heat of fusion of water (3.34 × 105 J/kg) and the specific heat capacity of ice (2.01 × 103 J/kg/K), which means even a trace meltwater volume can act as a potent heat source in the ice sheet energy equation. One kilogram of refreezing meltwater releases sufficient energy to warm 100 kg of ice by ∼1.6 K, corresponding to a decrease in effective viscosity of ∼11% [Paterson, 1994]. The dual‐column parameterization of Phillips et al. [2010], whereby overlapping continua of water and ice transfer energy via a heat exchange term, offers one enthalpy‐based formulation for this process. This framework suggests that the efficiency of heat transfer from water to ice is far more sensitive to the characteristic spacing of elements of the glacier hydrologic network than the temperature gradient between water and ice [Phillips et al., 2010]. Given the ability for water‐filled crevasses and moulins to penetrate through sub‐freezing ice to reach the ice sheet bed [Boon and Sharp, 2003; van der Veen, 2007], characteristic crevasse and moulin spacing strongly modulates the efficiency of cryo‐hydrologic warming [Phillips et al., 2013].

Within the limited regions of West Greenland for which observations are available, both crevasse extent and moulin density have increased between 1985 and 2009 [Colgan et al., 2011b; Phillips et al., 2011]. The increase in moulin density likely results from an increase in specific meltwater runoff, whereby the ice sheet area associated with the threshold supraglacial runoff required for moulin establishment is increasing in a warming climate. Under the present climate, supraglacial streams already support moulin formation up to elevations of 1700 m in South Greenland, close to the elevation limit of moulins derived from supraglacial lake hydrofracture [Doyle et al., 2014; Smith et al., 2015]. Observed increases in crevasse extent likely stem from a combination of: (1) thinning and steepening of the ablation area, which increases tensile stress while decreasing lithostatic stress, (2) the reorientation of regional stress fields due to outlet glacier acceleration, and (3) more abundant surface meltwater available for hydrofracture [Lampkin et al., 2013]. Given that increased crevasse extent is both a response to changes in ice flow and a driver of cryo‐hydrologic warming, a potential crevasse‐warming‐velocity feedback has been postulated [Colgan et al., 2014; Figure 5]. Such a feedback has nontrivial implications. For instance, crevasse extent can increase by >10% on decadal timescales [Colgan et al., 2011b], water‐filled crevasses 20 m apart warm surrounding ice temperatures 3 times faster than when spaced 60 m apart [Phillips et al., 2010], and ice at 0°C is 10 times less viscous than ice at −10°C [Paterson, 1994].

Figure 5

Linkages between increased meltwater runoff, enhanced glaciological processes, and Greenland mass loss. A potentially nontrivial crevasse‐warming‐velocity feedback is highlighted [after Colgan et al., 2014].

For cryo‐hydrologic warming to become relatively efficient in the ice sheet interior, whereby ice and water temperatures reach a transient equilibrium on decadal timescales, a characteristic spacing of ∼100 m between elements of the glacier hydrologic network is required [Phillips et al., 2010]. In the absence of other processes, widely spaced supraglacial lake drainage events are unlikely to result in such a characteristic glacier hydrologic network length scale of ∼100 m. Efficient cryo‐hydrologic warming of the ice sheet interior therefore relies on the inland propagation of both crevasse fields and supraglacial stream‐derived moulins found in the contemporary ablation and lower accumulation (“wet snow zone”) areas to the ice sheet interior. Given that moulin formation is, and will likely continue to be, facilitated by an upward migration in equilibrium line altitude and consequent increase in bare and wet snow facie extents [Phillips et al., 2011; McGrath et al., 2013], and crevasse formation is, and will likely continue to be, facilitated by steepening of the ice sheet and increasing available meltwater for hydrofracture [van der Veen, 1998; Parizek and Alley, 2004], a continued inland migration of these supraglacial features accompanying the ongoing inland migration of snow facies does not appear implausible [Lampkin et al., 2013; McGrath et al., 2013]. Indeed, the potential deformational velocity gradients stemming from the periphery of the ice sheet warming to pressure melting point temperature before the interior of the ice sheet would likely challenge the continuum mechanics of ice, facilitating an inland propagation of fracturing and crevassing.

Numerical Model

We use a simple numerical model to provide a first‐order assessment of the potential changes in ice sheet form and flow associated with a complete thermal‐viscous collapse over five centuries. While we employ the end‐member scenario of a complete thermal‐viscous collapse, we note that even a partial collapse would realize the majority of the enhanced shear and volumetric decrease inferred by this simulation. To acknowledge substantial uncertainty in three key parameters (ice temperature, Wisconsin‐Holocene viscosity contrast, and surface boundary forcing), we employ a Monte Carlo statistical scheme. In this approach, an ensemble of ice sheets, each of which reflects the contemporary ice sheet form, flow and thermal structure under unique parameter and forcing perturbations, are subjected to a common ice temperature forcing, in which both ice temperatures and their warming rate are prescribed. In order to isolate the effects of thermal‐viscous collapse, we maintain surface forcing, basal sliding, and terminus position constant after model spin‐up.

Formulation

We employ a two‐dimensional (2D) ice flow model formulated according to the shallow ice approximation. We prescribe present‐day ice sheet temperatures from the three‐dimensional (3D) ice temperature field predicted by the Community Ice Sheet Model [CISM; Price et al., 2011]. To account for the influence of polythermal vertical ice temperature profiles on horizontal ice flux, horizontal velocity components are obtained by vertical integration of the first‐order momentum equation corresponding to the shallow ice approximation. Thus, although mass conservation is solved in 2D, we resolve the 3D velocity field of the ice sheet. The model is one‐way thermo‐mechanically coupled, whereby velocity computations do account for vertical temperature variations through the influence of temperature on the flow law parameter (effective ice viscosity). Given that solving ice temperature requires an iterative solution to the energy equation to be interlaced with an iterative solution to the momentum equation in each time‐step, one‐way thermal‐mechanical coupling has been previously adopted in 2D ice flow models to reduce computational burden [Marshall et al., 2005; Colgan et al., 2012]. The main advantage of one‐way coupling is that the computational efficiency associated with adopting one‐way thermal‐mechanical coupling permits us to perform multiple ice sheet scale simulations in Monte Carlo fashion, as well as have tight control over the ice temperature field, which we are forcing to evolve in a deliberate fashion that is not fully consistent with physical processes incorporated in conventional thermo‐mechanical models [Macpherson, 2013]. To acknowledge the uncertainty associated with the CISM 3D ice temperature field we employ, the initial ice temperature field [T(x,y,z)] is the first of three key uncertain parameters we perturb in each simulation. While the perturbation is uniform across all computational nodes within a simulation (x,y,z), the perturbation is random across simulations, normally distributed with a mean of 0 K and a standard deviation of 2 K, which is taken to represent uncertainty in the CISM‐derived ice temperature field [Price et al., 2011].

The vertically integrated formulation of the transient glacier continuity equation takes the form [Hooke, 2005]: where H is ice thickness, b is surface forcing, and and are depth‐averaged ice velocity in the x and y directions, respectively. Depth‐averaged ice velocities are calculated as where u and v denote basal sliding velocities in the x and y directions, ∂u/∂z and ∂v/∂z denote horizontal strain rates in the x and y directions, and z is ice elevation above bedrock. The limits of the vertical integration are bedrock (h) and ice surface (h) elevation. Numerical integration allows vertical profiles of both ice temperature and velocity to be resolved at each (x,y) node. Horizontal strain rates are calculated as [Hooke, 2005]: where ρ is bulk glacier density (taken as 910 kg/m3), g is gravitational acceleration (9.81 m/s2), α is absolute surface gradient, and A is effective ice viscosity. We assume the exponent (n) in the nonlinear flow law describing the relation between stress and strain to be 3 [Glen, 1958]. Effective ice viscosity, which is nonlinearly dependent on ice temperature (T), is calculated according to Paterson [1994] as where A is a reference effective viscosity (5.47 × 1010/Pa3/a when T ≥ 263 K and 1.14 × 10−5/Pa3/a when T < 263 K), Q is the creep activation energy of ice (139 kJ/mol when T ≥ 263 K and 60 kJ/mol when T < 263 K), and R is the ideal gas constant (8.314 J/mol/K). The dimensionless enhancement factor E accounts for the viscosity contrast between relatively less viscous Wisconsin ice (deposited prior to 11.7 kaBP) underlying relatively more viscous Holocene ice within the deepest layers of the ice sheet where the majority of shear occurs [Huybrechts, 1994]. This enhancement factor is the second of three key uncertain parameters involved in the model. In each simulation, we randomly prescribe E based on a normal distribution with a mean of 4 and a standard deviation of 1 [Budd et al., 2013; Lüthi et al., 2015]. E is taken as uniform across all computational nodes within a simulation (x,y,z). Although the transition depth between Wisconsin and Holocene ice is spatially heterogeneous across the ice sheet [Huybrechts, 1994], it is generally sufficiently deep that it influences the majority of shear [Reeh, 1985; Paterson, 1994]. For this reason, as well as avoiding the implementation of Lagrangian tracking for specific ice parcels whose relative depths change through time with local ice sheet thickness, we consider the full ice column “enhanced” by Wisconsin viscosity.

Boundary Conditions and Numerical Method

The surface boundary condition of the ice flow model is the assumption that the ice surface is a free surface, meaning it experiences zero stress [Hooke, 2005]. The surface forcing of the ice flow model is informed by the 1961–1990 mean surface mass balance simulated by the regional climate model MAR version 3.2 forced by ERA‐40 reanalysis data [Fettweis et al., 2013; Figure 6]. MAR surface mass balance, which is calibrated against in situ observations in both the ablation and accumulation areas of the ice sheet, accurately captures the spatial distribution of mass input to the ice sheet, and has been previously employed for surface forcing of 2D ice flow models [Quiquet et al., 2012; Goelzer et al., 2013; Colgan et al., 2015]. In order to achieve a fully transient spin‐up, with minimized trends in ice sheet volume, we decrease this surface mass balance field by a uniform 40% across the ice sheet. Without this 40% decrease in reference period surface mass balance, modeled ice sheet volume increases by approximately 2% per century during spin‐up. With this surface mass balance tuning, it instead decreases (corresponding to thinning) by 0.05% per century during spin‐up. Near‐term prognostic simulations often hold ice geometry constant during spin‐up, allowing the force balance to evolve within a prescribed geometry [e.g., Price et al., 2011]. Fully transient spin‐ups, whereby form evolves with flow, are generally limited to longer‐term diagnostic simulations [e.g., Solgaard et al., 2011]. Here, we impose observed ice sheet geometry as the target for a fully transient spin‐up for a near‐term prognostic simulation.

Figure 6

Mean reference period (1961–1990) surface mass balance over the Greenland ice sheet (in mWE/a), generated by the regional climate model MAR version 3.2 forced by ERA‐40 reanalysis data [Fettweis et al., 2013].

Rather than reflecting anomalous surface mass balance during the reference period, the substantial tuning of modeled surface mass balance compensates for inherent imperfections in the structure of our ice flow model, such as the assumption that n = 3 in interior regions of the ice sheet. In regions of relatively low deviatoric stress, such as near flow divides, ice flow is recognized to approach n = 1 [Pettit and Waddington, 2003]. Relative to n = 1 flow, n = 3 flow overestimates divide thickness and convexity. During spin‐up, the transient thickening in the ice sheet interior resulting from the n = 3 assumption is effectively compensated for, by decreasing surface mass balance. Similarly, as discussed below, the horizontal resolution of our model (Δx = Δy = 5 km) results in a slight underestimation of outlet glacier velocities, which also contributes to a positive bias in ice sheet mass balance during spin‐up. Given its inherent role in compensating such ice dynamic biases, we refer to b as “surface forcing” within the numerical framework, rather than strict “surface mass balance.”

The surface forcing parameter is the third of three key uncertainty parameters examined in this study. In each simulation, we randomly perturb the surface forcing (b). Similar to the ice temperature perturbation, this perturbation is uniform within a simulation, and random across simulations, in a normal distribution with a mean of 0% and a standard deviation of 25%. This means the 2σ uncertainty in surface forcing we assume (50%) is larger than the tuning we apply to MAR surface mass balance (40%). Simply put, we must acknowledge that the relative uncertainty in prescribed surface forcing is larger than the relative tuning we employ. We note that the relative uncertainty in untuned MAR surface mass balance is closer to 10% at ice sheet scale [Fettweis et al., 2013; Colgan et al., 2015]. While the relatively high uncertainty we assume in surface forcing contributes to a conservative ensemble spread in ice sheet volume, it also ensures that changes in ice sheet volume due to thermal‐viscous collapse must be larger than changes in ice sheet volume due to ±50% (or ±2σ) trends in surface forcing in order to be discernible.

The bottom boundary condition of the ice flow model consists of a prescribed basal sliding velocity. During spin‐up, at each computational node (x,y), we iteratively update prescribed basal sliding velocities (u and v) based on the difference between ice surface velocities observed by satellite interferometry ( and ; Joughin et al., 2010) and modeled ice surface velocities (u and v) according to where k denotes iteration. The quadratic summation explicitly acknowledges that absolute surface and basal velocities are comprised of both x‐ and y‐directional components. The dimensionless scaling parameter β forces the iterative update term in Equation (5) to decrease to 1% of the discrepancy between observed and modeled surface velocities by the end of spin‐up (Figure 7). We assume no forcing in the form of a basal mass balance, and no change in magnitude and spatial distribution of basal sliding after spin‐up to a transient equilibrium. This approach effectively attributes any discrepancy between modeled and observed surface velocities to basal sliding, regardless of whether basal ice is at the pressure melting point temperature. While undoubtedly an approximation, this approach does ensure that the transient velocity field very accurately reproduces the observed velocity field at the end of spin‐up.

Figure 7

Exponential growth of the scaling parameter (β) over the 200‐year transient spin‐up, in order to decay the iterative update of basal sliding velocity to extinction at the end of spin‐up (Equation (5)).

The lateral (downstream) boundary condition of the ice flow model is a prescribed zero thickness Dirichlet boundary condition (H = 0), at one grid spacing (5 km) beyond the contemporary ice sheet margin [e.g., Svendsen et al., 2014]. This boundary condition, which is imposed throughout spin‐up, forcing and spin‐down, effectively prescribes an upper limit to ice sheet area during our simulations [Saito and Abe‐Ouchi, 2010]. We note that recent decreases in surface mass balance have caused both land‐ and marine‐terminating portions of the ice sheet margin to retreat at a rate of >10 m/a [Kargel et al., 2012]. As an outcome of our surface forcing, we simulate analogous minor margin retreat during spin‐up. We therefore contend that 5 km is a reasonable upper limit for margin advance during transient forcing. The contemporary ice sheet area is defined by nodes at which observed ice thickness is greater than zero [Bamber et al., 2013]. As we must interpolate observed ice thickness values to a relatively coarse grid, which exceeds the width of most outlet glaciers, perimeter cells typically incorporate a mix of marine‐ and land‐terminating ice. Our prescription of observed margin and termini positions isolates the change in ice sheet volume due to thermal‐viscous collapse from potential changes in ice sheet volume associated with tidewater glacier mechanisms. Analogous fixed ice sheet boundary assumptions have been invoked to isolate transient changes in ice sheet volume due to changes in ice thickness stemming from specific ice dynamic processes [Reeh, 1985; Saito and Abe‐Ouchi, 2010].

The equations describing transient ice flow (Section 3.1) are discretized in space using first‐order finite volume methods with grid spacing Δx and Δy of 5 km. The semi‐discrete set of ordinary differential equations is then numerically integrated through 22 irregularly spaced vertical levels (Δz, increasing exponentially from 0.03 H at the bed to 0.10 H at the surface), to yield a suite of 2D (depth‐integrated) equations. These equations describing ice flow are solved using a semi‐implicit scheme and a Picard iteration with a time‐step (Δt) of 0.075 years. As part of the delicate dance between numerical stiffness and computational efficiency, we adopt a relatively relaxed Picard convergence criterion (δH) of 1 m between time steps. This yields numerically stable solutions with a tractable computational demand on 2.8 GHz cores (with 3 GB RAM per core). Under transient ice temperature forcing, the computational burden increases from ∼4 processor hours per century at the end of spin‐up to ∼100 processor hours per century four centuries later. This increase of nearly two orders of magnitude in computational demand is due to the greater number of iterations required to converge substantially higher ice fluxes after the onset of thermal‐viscous collapse (Figure 8). The super‐transient nature of these momentum equation terms under prescribed warming makes an interlaced iterative solution of the energy equation computationally intractable with the resources available to us. In total, the 50 simulations presented here reflect ∼17,000 processor hours. As the Janus Supercomputer uses ∼27 W per processor with a power usage effectiveness of 1.2, and assuming a coal‐based carbon burden of ∼1 kg/kWh associated with Colorado electricity [EIA, 2014], these processor hours are associated with 550 kg of CO2 emissions.

Figure 8

Computational burden of n = 50 simulations over the 200‐year spin‐up to transient equilibrium and the subsequent 500‐year combined transient forcing and spin‐down period.

Monte Carlo Simulations

We explore the potential changes in ice sheet form and flow associated with a transition from polythermal to temperate ice conditions with an ensemble of 50 simulations. Following spin‐up, our ensemble of ice sheets reflects the form, flow, and thermal stat of the contemporary ice sheet in transient equilibrium. In each simulation, we vary three key uncertain parameters, as described above: (1) present‐day (or initial) ice temperature (T), (2) enhancement factor (E) to the flow law parameter, and (3) surface forcing (b). In each simulation within our ensemble, we initialize the ice flow model with observed ice geometry [Bamber et al., 2013]. Our simulations reflect a 700‐year period divided into three stages: (1) the 200‐year transient spin‐up period during which ice temperatures are prescribed from the CISM, (2) a 400‐year transient forcing period during which ice temperatures are warmed to the pressure melting point at an a priori imposed rate of 0.05 K/a, and finally (3) a 100‐year transient spin‐down during which ice temperatures are maintained at the pressure melting point. The rate of warming we prescribe (0.05 K/a) stems from our desire to warm the volume‐averaged temperature of the ice sheet 20 K in 400 years, from the present 252.0 K [Price et al., 2011] to the volume‐averaged pressure melting point of 272.2 K. To place this rate in context, warming rates of ∼0.1 K/a have been observed in firn in Northwest Greenland at 1400–2500 m elevation between 1954 and 2013 [Polashenski et al., 2014], and warming rates of ∼0.5 K/a have been inferred for basal ice in West Greenland at equilibrium line altitude between 2001 and 2007 [Phillips et al., 2013]. The imposed warming reflects an energy burden of 1.09 × 1023 J to bring the ice sheet to pressure melting point temperature throughout. We again emphasize that this is an end‐member scenario, similar in spirit to exploring the potential sea‐level rise contribution associated with a marine instability collapse of the West Antarctic ice sheet [Bamber et al., 2009], and we do not attempt to constrain its probability of occurrence. The specific end‐member scenario we explore may be regarded as a first‐order upper limit on the Greenland sea‐level rise contribution due to cryo‐hydrologic warming.

Inferred Volume Evolution

Following the 200‐year transient spin‐up of 50 simulations, the ensemble mean ice volume is 2.94 × 106 km3, which is within 1% of the than observed ice volume [2.96 × 106 km; Bamber et al., 2013]. We quantify the variation in specific ice sheet characteristics across simulations by 1σ ensemble spread. For ice volume, this corresponds to a ±2.78 × 104 km3 ensemble spread. After spin‐up, the ensemble mean surface velocity is 43 ± 4 m/a, which is 12% lower than observed mean ice surface velocity [49 m/a; Joughin et al., 2010]. We speculate that the majority of this velocity underestimation stems from the fact that the 5 km surface velocity field used for model evaluation is derived from a statistical interpolation of a 500 m velocity product, which samples outlet glaciers with velocities orders of magnitude higher than surrounding ice, while the 5 km modeled velocity field only coarsely resolves outlet glaciers. By virtue of being forced by CISM ice temperatures and their associated uncertainty, the ensemble of simulations represents a best estimate of volume‐averaged ice temperature [252.0 K; Price et al., 2011], with an ensemble spread of ±2 K. Thus, while the velocities of some coarsely resolved outlet glaciers are underestimated in comparison to observations, the post spin‐up ensemble is broadly representative of the contemporary form, flow, and thermal state of the ice sheet (Figures 9 and 10). We interpret an apparent bias toward a steeper ice sheet periphery than observed as stemming from the limitations associated with the shallow ice approximation, as well as a negatively tuned surface forcing. The inclusion of longitudinal coupling stresses typically generates ice masses that are thinner than generated by the shallow ice approximation alone [Marshall et al., 2005].

Figure 9

Greenland ice sheet form, flow, and thermal state: ensemble mean ice surface elevation (a in m), ice surface velocity (b in m/a), and depth‐averaged ice temperature (c in °C) at the end of the 200‐year transient spin‐up.

Figure 10

Difference between ensemble mean ice surface elevation (a in m) and ice surface velocity (b in m/a) at the end of the 200‐year transient spin‐up relative to observed ice surface elevation [Bamber et al., 2013] and ice surface velocity [Joughin et al., 2010].

After the 200‐year transient spin‐up to equilibrium, ice temperatures are increased to the pressure melting point temperature during a 400‐year transient forcing period, and then maintained at the pressure melting point temperature for a 100‐year transient spin‐down. Over the 500‐year combined forcing and spin‐down periods, in response to the relaxation of effective ice viscosities, ensemble mean ice surface velocity increases by a factor of 2.9 (from 43 ± 4 to 126 ± 17 m/a) at 335 years post spin‐up (Figure 11; Movie S1, Supporting Information). While ice sheet volume is virtually constant during the 200‐year spin‐up, it decreases by 5 ± 2% by the end of the 500‐year combined forcing and spin‐down period. This corresponds to a cumulative sea‐level rise contribution of 33 ± 18 cm (Figure 12). Continued drift in mean velocity at 500 years post‐spin‐up indicates that the ice sheet has not reached a new equilibrium by that time. Given the aggressive end‐member ice temperature forcing scenario we adopt, we suggest that the full sea‐level rise contribution associated with a thermal‐viscous collapse of the ice sheet cannot be realized within 500 years.

Figure 11

Probability density time series of ensemble spread of n = 50 simulations in prescribed ice temperature (a), mean surface ice velocity (b), and ice volume (c), over the 200‐year spin‐up to transient equilibrium and the subsequent 500‐year combined transient forcing and spin‐down period.

Figure 12

Probability density time series of Greenland sea‐level rise contribution of n = 50 simulations over the 200‐year spin‐up to transient equilibrium and the subsequent 500‐year combined transient forcing and spin‐down period.

Summary Remarks

We have considered potential changes in Greenland ice sheet form and flow associated with increasing ice temperature, and relaxing effective ice viscosity, over a multicenturial timescale. We have defined “thermal‐viscous collapse” as the decrease in effective ice viscosity due to a transition from the polythermal ice temperature distribution characteristic of the Holocene to an isothermal ice sheet at the pressure melting point temperature throughout. Our conceptual model proposes three key elements required to initiate a partial thermal‐viscous collapse: (1) sufficient energy available in projected ice sheet meltwater runoff, (2) routing of meltwater to the bed of the ice sheet interior, and (3) efficient energy transfer from meltwater to the ice via cryo‐hydrologic warming. We suggest that initiating a partial thermal‐viscous collapse, whereby the deepest 15% of the ice sheet is warmed to the pressure melting point temperature, appears to be thermodynamically plausible within the next four centuries. While plausible, we do not attempt to constrain the probability of occurrence of this end‐member scenario, and merely seek to provide a useful first‐order upper limit calculation. As the majority of all deformational shear occurs within this deepest portion of the ice sheet, such a partial thermal‐viscous collapse would realize the majority of the enhanced shear and flow associated with a complete thermal‐viscous collapse.

A first‐order ice flow model was employed, couched in a Monte Carlo approach, to assess the potential sea‐level rise contribution associated with a complete thermal‐viscous collapse of the ice sheet over the next five centuries. Under this end‐member scenario, we find a decrease in ice sheet volume of 5 ± 2% by 500 years after the initiation of collapse. This is equivalent to a cumulative sea‐level rise contribution of 33 ± 18 cm from thermal‐viscous collapse alone. By maintaining surface forcing, basal sliding, and terminus position constant after simulation spin‐up, we have effectively ignored these additional mechanisms of ice sheet response. Over the next five centuries, these neglected mechanisms and their interconnected feedback responses will also influence the mass loss from the Greenland ice sheet. We note that the full sea‐level rise contribution associated with a complete thermal‐viscous collapse, and relaxation of effective ice viscosities, cannot be realized within five centuries. Given that the comparatively subtle rheological contrast between Wisconsin and Holocene ice results in a 15% change in mean ice sheet thickness over millennia [Reeh, 1985], the full volume decrease associated with thermal collapse is much likely larger than presented here, and would likely be realized over several millennia.

While the notion of a persistent upward and inward migration of the ice sheet runoff limit, ultimately initiating thermal‐viscous collapse, may appear radical, it is important to recall that the present rate of Arctic warming is virtually unprecedented [Miller et al., 2013]. Observations suggest that the West Greenland ablation area was ∼35 km wide in ca. 1990, grew to ∼80 km wide by ca. 2011, and is expected to be >120 km wide by 2020 [McGrath et al., 2013]. In the context of such highly transient recent change, as well as projected post‐2100 committed change, it seems difficult to expect a natural brake on the seemingly inexorable upward migration of snow zones. The Barnes ice cap, which is now entirely comprised of ablation area after its equilibrium line ascended above its highest elevations in the past 20 years, may serve as a micro‐analogue for the Greenland ice sheet [Gardner et al., 2012]. The notion of supraglacial stream networks reaching the ice cap summit would have presumably been greeted with skepticism decades ago [Baird et al., 1950]. Perhaps similarly, the notion of active supraglacial and subglacial hydrologic networks penetrating above 2500 m elevation in Greenland may seem less conjectural in subsequent centuries.

Movie S1

Click here for additional data file.