A simple equation for the melt elevation feedback of ice sheets

In recent decades, the Greenland Ice Sheet has been losing mass and has thereby contributed to global sealevel rise. The rate of ice loss is highly relevant for coastal protection worldwide. The ice loss is likely to increase under future warming. Beyond a critical temperature threshold, a meltdown of the Greenland Ice Sheet is induced by the self-enforcing feedback between its lowering surface elevation and its increasing surface mass loss: the more ice that is lost, the lower the ice surface and the warmer the surface air temperature, which fosters further melting and ice loss. The computation of this rate so far relies on complex numerical models which are the appropriate tools for capturing the complexity of the problem. By contrast we aim here at gaining a conceptual understanding by deriving a purposefully simple equation for the self-enforcing feedback which is then used to estimate the melt time for different levels of warming using three observable characteristics of the ice sheet itself and its surroundings. The analysis is purely conceptual in nature. It is missing important processes like ice dynamics for it to be useful for applications to sea-level rise on centennial timescales, but if the volume loss is dominated by the feedback, the resulting logarithmic equation unifies existing numerical simulations and shows that the melt time depends strongly on the level of warming with a critical slowdown near the threshold: the median time to lose 10 % of the present-day ice volume varies between about 3500 years for a temperature level of 0.5 C above the threshold and 500 years for 5 C. Unless future observations show a significantly higher melting sensitivity than currently observed, a complete meltdown is unlikely within the next 2000 years without significant ice-dynamical contributions.

Numerical simulations suggest that a decline of the Greenland Ice Sheet is inevitable once its surface temperature permanently exceeds a certain threshold (Charbit et al., 2008;Greve, 2000;Huybrechts and De Wolde, 1999;Huybrechts et al., 2011;Ridley et al., 2005Ridley et al., , 2010Robinson et al., 2012;Solgaard and Langen, 2012). If and when this temperature threshold is passed depends critically on past and future greenhouse gas emissions Goelzer et al., 2013;Gregory et al., 2004a;Rae et al., 2012). Even if emissions were reduced to zero, temperatures would not drop significantly for thousands of years because of the long lifetime of anthropogenic CO 2 in the atmosphere and reduced oceanic heat uptake if oceanic convection is extenuated (Allen et al., 2009;Solomon et al., 2009;Zickfeld et al., 2013). This implies a possible commitment of a melt down of the Greenland Ice Sheet in the near future, which would eventually raise global sea-level by more than 7 m (Gregory et al., 2004a). Whether this occurs on a multi-centennial or rather a multi-millennial timescale is of relevance for coastal planning.
In this article we first recap the Vialov profile and add a simple representation of the melt elevation feedback towards a governing equation for a steady-state ice sheet in 1 dimension, then we derive the critical warming threshold for the existence of an ice sheet in this simple model (Sect. 2). In Sect. 3 we derive a simple time evolution equation for the decay of the ice sheet after surface temperatures have exceeded the threshold. Finally we use observational estimates of the three characteristics that enter the model to estimate the decay time of the ice sheet under melting above the threshold (Sect. 4). Here solid ice discharge is neglected as well as any other ice sheet dynamics (Andresen et al., 2012;Howat and Eddy, 2012;Moon et al., 2012;Nick et al., 2009;Price et al., 2011;Straneo et al., 2011;Walsh et al., 2012). The framework that we introduce here can be used to include new physical processes that might be discovered in the future, e.g. potential changes in surface albedo through melting (Box et al., 2012) or aerosol-induced surface melt or the lack thereof (Polashenski et al., 2015).

Governing equation for shallow-ice steady states under melt elevation feedback
A non-linear threshold behaviour is generally associated with a fundamental self-enforcing feedback and thereby an associated system memory . For the Greenland Ice Sheet, such a feedback is given by the interaction between surface elevation and surface melting (Weertman, 1961). For illustration, we include this feedback in a well-established highly idealized ice profile of an ice sheet in 1 dimension, the so-called Vialov profile (Vialov, 1958). We introduce the melt elevation feedback in the simplest possible way by assuming that the surface melt rate depends linearly on the surface temperature and that the temperature decreases linearly with the height of the ice surface following a constant atmospheric lapse rate.

Governing equation
We consider a highly simplified flow line model for an isothermal ice sheet grounded on a flat and rigid bed. The solution of the shallow-ice approximation in 1 dimension for the ice sheet elevation under these simplifying assumptions is the Vialov profile: where h m is the maximum surface elevation and n is Glen's flow law exponent (Glen, 1955). x denotes the horizontal position and L the horizontal limit of the ice sheet. The inherent assumption of isothermal ice is a strong simplification, which needs to be kept in mind when interpreting the results.
The aim of this derivation is purposefully not a comprehensive representation of the ice flow but to derive a measure of the average height of the ice sheet and its dependence on changes in the surface mass balance. The surface mass balance is considered to be spatially and temporally constant at a value, a, which will later be considered to be dependent on the surface elevation and thereby temporally variable. The overall horizontal extension of the ice sheet is set to L, and it is thereby assumed that any ice flow across this point is calved off into icebergs. This situation represents a confined ice-bearing bedrock topography as in most of Greenland's interior (Howat et al., 2014). The mean surface elevation can then be computed to be It is proportional to the maximum surface elevation h m with a proportionality factor which only depends on the flow law exponent. The maximum surface elevation is determined by the surface mass balance a and the ice softness A: with ρ being the ice density and g the gravity constant. We normalize all three quantities by defining h ≡ ω · h m /h 0 , a ≡ a/a 0 and A ≡ A/A 0 , where a 0 is the accumulation rate on the ground, i.e. in the absence of an ice sheet, and A 0 = a 0 / (ρg) n ( · L) (n+1) with = h 0 /L being the typical height-to-width ratio. h 0 is the equilibrium line altitude of the considered ice sheet in the initial equilibrium situation. Values for a 0 , h 0 and L are later chosen to resemble the conditions of the Greenland Ice Sheet. The non-dimensional surface elevation, h, of the ice sheet can then be expressed as For the Vialov profile, m = 2(n + 1) where the Glen flow law exponent is commonly chosen to be around n = 3, which yields m = 8. We introduce the melt elevation feedback in its simplest form through a dependency of the surface melt rate on the surface elevation: with the atmospheric lapse rate > 0. γ denotes the melting sensitivity of the ice surface, i.e. the increase in surface melt rate per degree of warming, which is regularly measured and comprises a large number of physical processes (e.g. Box, 2013). For simplicity we rescale the surface mass balance by the constant ice softness parameter, A, to obtain h = (a 0 + γ ·h) 1/m . The steady-state solution for the surface elevation of the ice sheet is thus governed by the following equation: which has two positive solutions for h as long as the surface mass balance on the ground is negative, i.e. a 0 < 0. Note that the surface mass balance can be positive even if a 0 < 0. If the ice sheet is in an unstable configuration, a slight perturbation will either cause it to converge into the stable state with a positive surface mass balance or to melt completely. Our simple approach qualitatively captures the basic hysteresis behaviour of the Greenland Ice Sheet caused by the melt elevation feedback (Fig. 1, in which we have assumed the surface mass balance to depend linearly on temperature): For a given surface temperature, a stable state of the ice sheet (red line) annihilates an external perturbation in surface elevation by changes in surface mass balance (grey arrows). The unstable solution branch defines the basin of attraction for the stable state. A surface elevation that is lower than the unstable solution branch cannot be sustained. In that case the melting reduces the surface elevation to practically zero even without further external perturbation (grey arrows). Beyond a certain surface temperature threshold (vertical dotted line), no ice sheet can be sustained.

Critical surface mass balance in steady state
As illustrated in Fig. 1, there is a critical temperature above which the ice sheet is not sustainable. Let us denote the corresponding surface elevation by h c . The critical point (T c , h c ) has to fulfill two conditions, i.e. being a solution of the governing Eq. (7) and minimum of the function which we can determine by setting the derivative of F to zero. Consequently, Inserting this into the governing equation yields the critical surface mass balance at the ground: For illustrative purposes we have assumed a 0 to decline linearly with the surrounding temperature and plotted the solution of Eq. (7) against that temperature with an arbitrary offset in Fig. 1.

A simple temporal equation for the melt elevation feedback
Once the critical surface mass balance and surface elevation threshold (as derived in the previous Sect. 2) is transgressed, a meltdown of the ice sheet is inevitable in our conceptual model. Let us define the time τ α as the time it takes to melt a fraction α of the initial ice volume and the threshold temperature T c as the temperature above the pre-industrial level at which the surface mass balance becomes negative. Robinson et al. (2012) find a range of 0.8-3.2 • C for the threshold warming beyond which no ice sheet can be sustained on Greenland. Their best estimate for the threshold is 1.6 • C above pre-industrial level. The study uses a regional climate model of intermediate complexity (Robinson et al., 2010) coupled to the SICOPOLIS (SImulation COde for POLythermal Ice Sheets) ice sheet model (Greve, 1997). Using a different model combination, Ridley et al. (2010) find that in their model the ice sheet cannot be sustained for a warming of 2 • C. They combine the HadCM3 atmosphere-ocean general circulation model (Gordon et al., 2000) with an atmospheric resolution of 2.5 • ×3.75 • (Pope et al., 2000) to an ice sheet model of 20 km horizontal resolution (Huybrechts and De Wolde, 1999).
Some studies assume that the threshold is associated with a mean negative surface mass balance (Gregory et al., 2004b;Ridley et al., 2005;Toniazzo et al., 2004). In Fig. 2 we use 1.6 • C as a threshold value for both models because this value is given by Robinson et al. (2012)  percentage ice volume change, a constant horizontal ice surface area was assumed which renders the analysis conceptual in nature. Thus the quantitative interpretation of the melt times are subject to this additional simplification. For a fixed anomalous melt rate a 0 = −γ · T in response to an anomalous temperature increase T = T − T c above this threshold temperature, T c , the decay time without any feedbacks would be Since the surface temperature increases with decreasing elevation, this zero-order estimate for the decay time is higher than the actual value. As a first-order correction to the situation of fixed melting, let us assume that the anomalous surface mass balance behaves as where τ γ = 1/(γ · ). From the relation dh/dt = a, we then obtain if h ≡ h 0 − h is defined as the reduction in height. For a time-dependent melting induced by surface warming a 0 = −γ · T , the general solution of Eq. (13) is This equation corresponds to a linear response theory with the melting −γ · T as forcing and an exponential response function Linear response theory states that the convolution of Eq. (14) yields the linear response of the system (Good et al., 2011;Winkelmann and Levermann, 2013). Note that linear response theory is generally used as an approximation of a non-linear system to relatively weak forcing. In these circumstances the response function has to decline with time because it represents the history of the system's response to past perturbation. For example, if the response function was a declining exponential R(t ) = e −t , this would mean that the effect of forcing that occurred in the past, i.e. prior to the time t that is considered, becomes exponentially less relevant for the current system response. Here, however, the response function is increasing with time, which means that the past deviation from the steady state is amplified as expected near an unstable fixed point. The exponent 1/τ γ can be considered the Lyaponov exponent of the system. Given the boundary condition h(t = 0) = 0 for a constant temperature increase T , Eq. (14) becomes The decay time for a relative volume reduction of α is then given by where log denotes the natural logarithm. Equation (17) is denoted the decay time equation hereafter.

Estimating the melt time of the Greenland Ice Sheet from observables
In this simplified approach, the collapse time is thus a function of three observable quantities: the equilibrium line altitude, h 0 , the atmospheric lapse rate, , and the melting sensitivity to temperature, γ . The average equilibrium line altitude of the Greenland Ice Sheet is at about 1150 m (Box and Steffen, 2001). The observed range for the atmospheric lapse rate is estimated to be between 5 ± 2 • C km −1 (Fausto et al., 2009;Gardner and Sharp, 2009) and current estimates for the melting sensitivity scatter around 4.4±2 cm year −1 • C −1 (Box, 2013). In order to obtain an estimate of the decay time and the uncertainty around this estimate we use Eq. (17) and choose the lapse rate and melting sensitivity uniformly randomly from these observed intervals ( the melt time of the Greenland Ice Sheet, assuming uniform probability distributions for both and γ within the above intervals. Figure 2 shows the histograms of the time until 10 % of its present-day ice volume (corresponding to 0.7 m global sea-level rise) is melted for different warming scenarios. The melt time strongly depends on the level of warming beyond the temperature threshold: the median estimate varies from more than 2000 years for a warming of +1 • C to less than 500 years for a warming of +5 • C. Existing numerical simulations of a decay of the Greenland Ice Sheet (Ridley et al., 2010;Robinson et al., 2012) differ in their trajectories for the total ice volume, but exhibit a characteristic functional form when the relative ice volume is expressed as a function of the temperature anomaly above the critical temperature threshold (Fig. 2). This characteristic relation is captured by our first-order equation for the decay time, embedding the results from process-based models into a simple analytical framework. This approach provides a good approximation if, on the one hand, the volume loss is large enough for the melt elevation feedback to become relevant and, on the other hand, the melting dominates the ice loss in contrast to the dynamic ice discharge.
Since the simple equation provided here does not account for any dynamic discharge or even ice motion, the results from Eq. (17) strongly deviate from numerical simulations when the ice has time to adjust dynamically to the volume loss. This can be seen for a stronger ice loss of 50 % of the initial volume where the functional dependence between the decay time and the temperature anomaly clearly follows a different functional form than predicted by Eq. (17) (Fig. 3).
Since the melt time is a monotonically decreasing function of both the lapse rate and the melting sensitivity, the upper and lower limits of the estimates can be directly computed from the observed uncertainty interval of these quantities. However, the functional form of Eq. (17) introduces a Shown are the median (black line) and the likely (18-83 % percentiles, dark blue shading) and very likely (5-95 % percentiles, light blue shading) ranges for the time to melt 50 % of the presentday ice volume, estimated via the equation for the decay time τ α . The red crosses indicate the results from process-based model simulations by Robinson et al. (2012). specific structure into the histogram of the melt time which is highly skewed towards the low end (Table 1 and Fig. 4). For increasing warming levels the histogram shifts towards lower decay times. At the same time the histogram narrows and higher decay times become less frequent within the chosen parameter range (see description above).

Discussion and conclusion
Our estimate for the decay time captures the characteristic slowdown near the critical threshold as can be seen from the divergence of the decay time, τ α , in the limit of vanishwww.the-cryosphere.net/10/1799/2016/ The Cryosphere, 10, 1799-1807, 2016 ing warming above the threshold (Eq. 17). The simple equation of the decay time quantitatively reproduces the range given by simulations with process-based models. The relative speed-up of ice loss due to the melt elevation feedback (Fig. 5)  mospheric lapse rate = 5 • C km −1 and melting sensitivity γ = 4.4 cm year −1 • C −1 . The feedback becomes more dominant near the threshold compared to larger temperature increases for which the external climatic forcing is more relevant. The simple equation provided here is clearly limited in its applicability. The role of the ice material properties is comprised into one parameter, the melting sensitivity of the ice to a temperature increase at the surface. This sensitivity will in general vary not only with time but also spatially and due to the melting itself. Similarly, the feedback role of the surrounding climate is represented by only one parameter, the atmospheric lapse rate which will again vary spatially but also with time as the ice surface declines.
Ice dynamics are deliberately excluded in our simple conceptual approach in order to separate and characterize the melt elevation feedback. In reality, ice dynamics of course play an important role in the ice sheet mass balance: radar (ERS-2) and laser (ICESat) altimetry observations show that mass changes in Greenland were dominated by changes in the surface mass balance (SMB) between 1995 and 2001, and both SMB and dynamics contributed equally to mass loss from the Greenland Ice Sheet between 2001 and 2009 (Hurkmans et al., 2014). Fürst et al. (2015) estimate that 40 % of the recent loss (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) is due to an increase in ice dynamic discharge, 60 % due to changes in the surface mass balance. Their results suggest that the future volume loss from the Greenland Ice Sheet might be predominantly The Cryosphere, 10, 1799-1807, 2016 www.the-cryosphere.net/10/1799/2016/ caused by surface melting and dynamic discharge is limited by margin thinning and retreat. Some studies suggest (Graversen et al., 2010;Price et al., 2011) that the dynamic discharge from Greenland is strongly limited by the ice sheet's bottom topography, for which estimates yield an upper bound of approximately 5-13 cm during the next century. Over a period during which the ice loss is dominated by the feedback and the ice-dynamic effect is limited, our approach provides a quantitative estimate of the melt time based on observable quantities. Equation (17) can thus be used if new observations suggest an altered melting sensitivity or changes in the atmospheric response to Greenland ice loss.
For a temperature increase of 5 • C, which could be reached within this century (IPCC, 2013), the median rate of sea-level contribution is about 1.4 mm year −1 , which is about 4 times that of its current contribution of about 0.4 mm year −1 (Rignot et al., 2011;Shepherd et al., 2012). Even for extremely high temperatures, however, the Greenland Ice Sheet cannot melt infinitely fast -our results show that a complete disintegration within the next 2 millennia is highly unlikely unless ice dynamics effects become dominant or the melting sensitivity is significantly higher than currently observed. For a global mean temperature increase below 2 • C, as agreed upon during the 2015 Paris UNFCCC climate summit, the threshold temperature would only be exceeded mildly and the decay time of the Greenland Ice Sheet would be multimillennial.