No general stability conditions for marine ice-sheet grounding lines in the presence of feedbacks

The “marine ice-sheet instability” hypothesis continues to be used to interpret the observed mass loss from the Antarctic and Greenland ice sheets. This hypothesis has been developed for conditions that do not account for feedbacks between ice sheets and environmental conditions. However, snow accumulation and the ice-sheet surface melting depend on the surface temperature, which is a strong function of elevation. Consequently, there is a feedback between precipitation, atmospheric surface temperature and ice-sheet surface elevation. Here, we investigate stability conditions of a marine-based ice sheet in the presence of such a feedback. Our results show that no general stability condition similar to one associated with the “marine ice-sheet instability” hypothesis can be determined. Stability of individual configurations can be established only on a case-by-case basis. These results apply to a wide range of feedbacks between marine ice sheets and atmosphere, ocean and lithosphere.

T he observed mass loss from the Antarctic and Greenland ice sheets associated with the retreat of their grounding lines [1][2][3] continues to be explained 2-4 by the "marine icesheet instability" hypothesis 5 , even though many recent studies have questioned its validity [6][7][8][9] . This hypothesis proposed by Weertman 5 , suggests that a marine ice sheet that resides on a bed below sea level and slopes towards its interior is "inherently unstable". He further suggested that "A stable ice sheet can occur if the bed slopes away from the center of the ice sheet." Weertman arrived to this conclusion by considering an unconfined marine ice sheet flowing into an unconfined ice shelf that gains its mass by accumulating snow at its surface with the accumulation rate constant in space and time. In later studies 10,11 , a similar ice-sheet configuration was considered and the same stability condition was derived by means of a linear stability analysis 12 .
Ice sheets interact with all other components of the Earth system-atmosphere, ocean, sea ice, lithosphere, land surface, and biosphere. These interactions give rise to a number of feedbacks 13 ; one of them is a feedback between the net accumulation/ablation rate and the ice-sheet height [14][15][16] . The net accumulation/ablation rate at a given location on the surface of the ice-sheet is determined by the difference between snow accumulation and surface melting (ablation). If the annual snow precipitation exceeds annual surface melt (if melting occurs at all) this location experiences net accumulation; if the reverse is true, it experiences net ablation. Both processes, snow precipitation, and surface melting, depend on in situ atmospheric conditions. Snow precipitation depends on the atmospheric moisture content and the surface ablation rate is determined by the surface energy balance. The latter involves a number of processes, e.g., absorption of short-and long-wave radiation, reflection due to albedo, wind, and many others. Both the snow precipitation rate and ablation rate strongly depend on the surface elevation 17 . Such a dependence produces a feedback between the surface elevation and the net accumulation/ablation rate 15 .
While several studies 14,15,18 have considered the effects of a feedback between the surface elevation and behavior of landbased ice sheets, all existing studies of stability of marine ice sheets 5,10-12,19-21 did not take into account feedbacks between the external conditions and the ice-sheet characteristics. It is, thus, not a priori clear whether the results of these studies remain valid in the presence of feedbacks.
Here we show that the previously derived stability conditions of unconfined marine ice sheets that disregard the feedback between net accumulation/ablation rate and the surface elevation do not generally hold if the feedback is taken into account. Furthermore, there is no general stability condition which is determined by properties of steady-state configurations at the grounding line. Stability of specific configurations can be determined only individually, on a case-by-case basis. These results apply to other feedbacks in interactions of marine ice sheets with other components of the climate system-ocean (the dependence of sub-ice-shelf melting on ice thickness), continental crust and lithosphere (changes in bed elevation due to erosion/deposition, glaciostatic adjustment and or/changes in the local sea level).

Results
Accumulation/ablation rate dependence on surface temperature. Precipitation is determined by the atmospheric moisture content. A good proxy for the moisture content is the saturation vapor pressure, governed by the Clausius-Clapeyron equation, which in turn, depends on the surface temperature 17 . Ablation depends on the (near) surface temperature in two ways-the latter determines whether surface melting occurs (if the surface temperature is above the freezing point), and it determines the magnitude of the surface melting, the ablation rate. Although there are many complex processes that govern interactions between the ice-sheet surface and atmosphere, the atmospheric surface temperature can be used as a sole proxy parameter to estimate both precipitation and surface melting, and consequently the net accumulation/ablation rate.
Atmospheric temperature is a strong function of elevation: it decreases according to the lapse rate as elevation increases. As a result, the net accumulation/ablation rate implicitly depends on the ice-sheet surface elevation via the surface temperature.
Recent developments in the regional atmospheric modeling for polar regions have allowed for more accurate representation of the near-surface atmospheric processes 22,23 . In the absence of the direct continent-wide long-term observations, this has made results of these models' simulations a useful tool for analysis of the polar atmospheric conditions. For example, results of simulations of the regional atmospheric model, MAR, have been used to construct a parameterization of changes in the net accumulation/ablation due to changes in the surface elevation of Greenland Ice Sheet based on a probabilistic approach 24 .
In this study, the results of MAR simulations, for both Antarctica 25 and Greenland 26 are used to construct an empirical relationship between the net accumulation/ablation rate and the surface temperature. The simulated accumulation/ablation rates, _ a, as a function of the annual mean surface temperatures for the RCP 8.  26 ). a Net accumulation/ablation rate _ a(m/yr) as function of the surface temperature T ( ∘ C); green dots indicate fit described by expression (1). b Surface temperature T ( ∘ C) versus ice sheet surface elevation S (m); green lines are expression (2) with Γ = 9.8 ∘ C/km, T sl = −15 ∘ C for Antarctica, and T sl = −5 ∘ C for Greenland, respectively. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29892-3 mean surface temperature T is described by Eq. (1) (see the "Methods" section).
Results of MAR simulations (Fig. 1b) show a linear decrease of the annual mean surface temperature with elevation supporting the relationship (2). Also, these results indicate that simulated temperatures around Greenland (yellow dots) are consistently warmer at all elevations than around Antarctica (blue dots), and substantially warmer at sea level T s (S = 0). Consequently, considering the accumulation/ablation rate of Antarctica and Greenland together, and simulated for the high-end emission scenario, allows construction of an empirical relationship for a wider range of surface temperatures and surface elevations that can simultaneously capture the net accumulation and ablation regimes.
Relating the net accumulation/ablation to the surface temperature has other advantages, in addition to the physical basis of the chosen relationship. As (Fig. 1b) indicates, depending on the climate conditions temperature at the same surface elevations can be different; relating the net accumulation/ablation rate to the surface temperature and not surface elevation removes this ambiguity. A simple dependence on the annual mean surface temperature at sea level T sl makes the empirical relationship (1) together with (2) an ideal tool for theoretical and conceptual studies that aim to establish fundamental aspects of interactions between ice sheets and atmosphere, and also for paleo-climate studies as it eliminates the need to compute precipitation and the surface melting separately. "Methods" section and Supplementary  Fig. 3 illustrate how sea-level temperature T sl controls the behavior of _ aðT S ðSÞÞ, and demonstrate that the empirical relationships (1) and (2) can capture regimes characteristic of MAR simulations for Antarctica (blue dots) and Greenland (yellow dots in Fig. 1a).
The effect of a feedback on stability of marine ice sheets. A steady-state configuration of an unconfined marine ice sheet can be described by a solution of the problem (3) described in the "Methods" section, which is the same as a model used in ref. 10 . An example of a steady-state marine ice sheet, for which accumulation/ablation rate varies with the surface elevation according to relationships (1) and (2), is shown in Fig. 2a. Its grounding line x g is located on the up-sloping (or "retrograde") bed. According to both a stability condition derived by Schoof 12 and a more complex one that accounts for the bed curvature and the gradient of the accumulation/ablation rate _ a x derived by Sergienko and Wingham 20 , this grounding line position is stable.
These stability analyses, however, have been performed under the assumption that _ a varies only along the length of the ice sheet, as illustrated in Fig. 2b, and does not depend on any characteristics of the ice sheet, e.g., its surface elevation. To determine whether this steady-state configuration is stable in the presence of a feedback between the accumulation/ablation rate and the surface elevation, Eqs. (1) and (2), a numerical timedependent model simulation, in which the grounding line position is perturbed from its steady-state position, is performed. The numerical model solves the time-variant problem described by Eqs. (3a), (3c)-(3e). The mass-balance evolves in time according to Eq. (4), which accounts for the feedback between the net accumulation/ablation rate and the surface elevation i.e., _ a ¼ _ aðT S ðSÞÞ. (Details of the time-dependent simulations are described in Supplementary Methods 1). Snapshots of the timeevolving ice-sheet surface and the corresponding spatial distributions of _ aðT S ðSÞÞ are shown in Supplementary Fig. 6. Evolution of the simulated grounding line position is shown in Fig. 3a. The grounding line advances beyond its steady-state position (marked by the black dashed line) indicating that the steady-state position is unstable. This contradicts stability conditions 12,20 that do not account for feedbacks.
However, if _ a is treated as a function of the position, and not of the surface elevation, i.e., _ a ¼ _ aðxÞ as shown in Fig. 2b, the simulated advance of the grounding line from a similarly perturbed position terminates at its steady-state location, indicating that in this case the steady-state configuration is stable (Fig. 3b), which is in agreement with traditional analytical stability conditions 12,20 that do not account for feedbacks. This experiment illustrates that it is the presence of the feedback that changes stability of the grounding line from stable if _ a ¼ _ aðxÞ to unstable if _ a ¼ _ aðT S ðSÞÞ. To shed light on such a different behavior of the grounding lines of ice sheets, which accumulation/ablation rates vary either with the surface elevation or only with the distance along the ice sheet, we analyze a simplified problem (5) used to establish a more complex stability condition 20 for the case of the net accumulation/ablation rate being a function of the position, _ a ¼ _ aðxÞ. A linear stability analysis for the case of the feedback between the net accumulation/ablation rate and the surface elevation _ a ¼ _ aðT S ðSÞÞ follows the same standard approach 20 , and considers small, time-dependent perturbations around the steadystate solutions, e.g., x g ¼x g þx g ðtÞ, wherex g is a steady-state position of the grounding line andx g ðtÞ is a small time-variant perturbation. The evolution of the grounding line in the corresponding linearized perturbation problem (Supplementary Methods 2) is described byx g ðtÞ ¼x g ðt ¼ 0Þe Λt , where Λ is an eigenvalue. A stability condition for the case of _ a ¼ _ aðxÞ relies on properties of the eigenfunctions of the perturbation problem 27 and is determined by the sign of the largest eigenvalue-either the fastest growing, in the case of an unstable grounding line position, or the slowest decaying in the case of a stable grounding line position 20 . In the case of _ a ¼ _ aðT S ðSÞÞ no inferences about either the sign or magnitude of eigenvalues can be made; and the perturbation eigenvalue problem has to be solved explicitly to determine the largest eigenvalue. This implies that there is no general stability condition that depends on properties of the steady-state ice-sheet configuration in the presence of the feedback between the net accumulation/ablation rate and the surface elevation. As Supplementary Eq. (14) indicates, the absence of the general stability condition in the presence of the feedbacks is insensitive to specifics of the functional form of the empirical relationship described by Eqs. (1) and (2).
Results of a numerical solution of the linear-perturbation eigenvalue problem are illustrated in Fig. 4a that shows the first ten eigenvalues. The first eigenvalue is positive indicating that the small perturbations from the steady-state configuration grow with time, i.e., the steady-state configuration and its grounding line position is unstable. This is in agreement with the results of the temporal evolution of the grounding line perturbed from its steady-state location (Fig. 3a). Eigenvalues for the case of _ aðxÞ are shown in Fig. 4b. The largest eigenvalue is negative, indicating that the small perturbations from the steady-state solution decay in time, and the steady-state solution is stable. This is again in agreement with the results of the temporal evolution of the grounding line perturbed from its steady-state location (Fig. 3b).

Discussion
In order to investigate how feedbacks may affect stability conditions of marine ice sheets, we have constructed an empirical relationship between the surface temperature, which varies linearly with elevation according to the lapse rate, and the net accumulation/ablation rate. The constructed formulation uses the fundamental relationships between the atmospheric surface temperature and precipitation (Clausius-Clapeyron equation) and the surface temperature and ice-sheet surface melting (energy balance). Expressed in terms of temperature at the sea level, T sl this formulation can be an efficient, physically based parameterization of the net accumulation/ablation for paleo-climate simulations.
Numerical and analytical analyses of the marine ice-sheet stability performed with this empirical relationship produce results markedly different from those of previous analyses that did not account for any feedbacks. Stability conditions derived for the cases with no feedbacks between the ice-sheet characteristics (e.g., the ice-sheet surface elevation) and other parameters (e.g., the net accumulation/ablation rate) do not hold in general if such feedbacks are present. Moreover, in the demonstration of this considered here, there are no general stability conditions that depend on characteristics of the steady-state configurations at the grounding line. In the presence of feedbacks, stability of a particular ice-sheet steady-state configuration can be determined only by either numerically simulating evolution of this configuration perturbed from its steady state, or by solving a perturbation eigenvalue problem to determine the largest eigenvalue, which sign indicate stability (all eigenvalues are negative) or instability (at least one eigenvalue is positive).
A physical interpretation of these results is the following: If the net accumulation/ablation rate depends on the surface elevation then the specific details of such a dependence (e.g., how _ a varies with S), together with the characteristics of the steady-state configuration, determine whether this steady-state configuration is stable or unstable. Changes in net accumulation/ablation rate  caused by changes in surface elevation due to small perturbations from steady-state configurations could suppress or amplify these perturbations, and consequently stabilize or destabilize the grounding line.
The presented analysis can be applied to other processes that include feedbacks. If sub-ice-shelf melting occurs at the grounding line, the melting rate is determined by the pressure melting point, which depends on the ice thickness. This introduces a feedback between ice thickness and sub-ice-shelf melting at the grounding line. If the bed elevation varies due to, for instance, glaciostatic adjustment and erosion/deposition processes, or due to changes in the local sea level, the surface elevation and temperature vary as a result of these processes. In all these circumstances the presented results remain valid.
A number of processes that have feedbacks with the ice-sheet characteristics are in play on the Antarctic and Greenland ice sheets 13 . Consequently, their stability conditions are not determined by the bed slopes at their grounding lines alone as suggested by the traditional "marine ice-sheet instability hypothesis" 5 . Instead, stability conditions depend on the interplay between buttressing (if the ice shelves are confined), basal conditions (bed elevation, slope, curvature and sliding properties) and feedbacks between the ice sheet configuration and atmosphere, ocean, and lithosphere interactions.

Methods
Empirical relationship between _ a and S. The fitting (green) line in Fig. 1a has the following expression where a 1 , a 2 , T 0 , and σ are fitting parameters and T S is the temperature at the surface elevation S, which is where T sl is the temperature at sea level and Γ is the lapse rate. In this study, we use the dry adiabatic lapse rate Γ = 9.8 ∘ C/km. In general, surface temperature at the sea level depends on insolation, which is function of latitude. In this study, we disregard this dependence and treat sea-level surface temperature T sl as constant.
The magnitude of the sea-level temperature T sl controls whether an ice sheet can experience net ablation or it experiences only net accumulation. The steadystate configuration shown in Fig. 2a has been computed with T sl = −4 ∘ C; the corresponding net ablation/accumulation rate shown in Fig. 2b as a function of distance, in Supplementary Fig. 4a as a function of surface elevation and in Supplementary Fig. 4b as a function of surface temperature experiences net ablation (negative values of _ a) at the low elevations close to the grounding line. These results indicate that this relationship captures MAR simulations for Greenland in warmer climate of the RCP 8.5 scenario (yellow dots in Fig. 1a). For colder sea-level temperature T sl = −10 ∘ C, a steady-state configuration computed with Eqs. (1) and (2) does not experience net ablation (Supplementary Fig. 3). As Supplementary Fig. 3d illustrates, for this value of T sl the constructed empirical relationship captures the behavior of _ a shown in Fig. 1a obtained in MAR simulations for Antarctica (blue dots).
Marine ice-sheet model. We use the model of ref. 10 . As in many previous studies of unconfined marine ice sheets flowing into unconfined ice shelves, the momentum balance of the ice shelf can be integrated analytically and the problem is reduced to one of the marine ice-sheet grounded parts only 28 . In a steady-state it is where u(x) is the depth-averaged ice velocity, h(x) ice thickness, B(x) is bed elevation (negative below sea level and positive above sea level), S = h + B is the surface elevation, A −1/n is the ice stiffness parameter (assumed to be constant), n is the exponent of Glen's flow law, g is the acceleration due to gravity, τ b ¼ C u j j mÀ1 u is basal shear with constant parameters C and m = 1/n, x d and x g are positions of the ice divide (taken here to be x d = 0) and the grounding line, respectively, and _ aðT S Þ is described by (1), with T S described by (2); g 0 is the reduced gravity defined as g 0 ¼ δg, where δ ¼ ρ w Àρ ρ w is the buoyancy parameter. In the case of temporal evolution of the ice-sheet configuration, the mass balance (3b) takes the form Simplified model. An approximation of (3) with (4) instead of (3b) written in terms of the ice thickness h and the ice flux q = uh is To determine stability conditions, a linear stability analysis of (5) is performed by considering small perturbations from the steady-state solutions. The details of the linear stability analysis in the presence of a feedback between the net accumulation/ ablation rate and the ice-sheet surface, i.e., _ a ¼ _ aðT S ðSÞÞ are described in Supplementary Methods 2; the details of the linear stability analysis in the absence of such a feedback i.e., _ a ¼ _ aðxÞ are described in ref. 20 .

Data availability
Outputs of MAR simulations for Antarctica 25