the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The contrasting response of outlet glaciers to interior and ocean forcing
John Erich Christian
Alexander A. Robel
Cristian Proistosescu
Gerard Roe
Michelle Koutnik
Knut Christianson
The dynamics of marineterminating outlet glaciers are of fundamental interest in glaciology and affect mass loss from ice sheets in a warming climate. In this study, we analyze the response of outlet glaciers to different sources of climate forcing. We find that outlet glaciers have a characteristically different transient response to surfacemassbalance forcing applied over the interior than to oceanic forcing applied at the grounding line. A recently developed reduced model represents outletglacier dynamics via two widely separated response timescales: a fast response associated with groundingzone dynamics and a slow response of interior ice. The reduced model is shown to emulate the behavior of a more complex numerical model of ice flow. Together, these models demonstrate that ocean forcing first engages the fast, local response and then the slow adjustment of interior ice, whereas surfacemassbalance forcing is dominated by the slow interior adjustment. We also demonstrate the importance of the timescales of stochastic forcing for assessing the natural variability in outlet glaciers, highlighting that decadal persistence in ocean variability can affect the behavior of outlet glaciers on centennial and longer timescales. Finally, we show that these transient responses have important implications for attributing observed glacier changes to natural or anthropogenic influences; the future change already committed by past forcing; and the impact of past climate changes on the preindustrial glacier state, against which current and future anthropogenic influences are assessed.
Marineterminating outlet glaciers drain large portions of the Greenland and Antarctic ice sheets, conveying ice from interior catchments to the ocean. Their dynamic response to a changing climate is a critical component of projections of sealevel rise (e.g., Aschwanden et al., 2019), and dramatic increases in discharge over recent decades have been observed in Greenland (e.g., Howat et al., 2008; Moon and Joughin, 2008; Moon et al., 2012) and Antarctica (e.g., Pritchard et al., 2009; Miles et al., 2013; Cook et al., 2016). Ocean forcing is thought to play a major role in observed change (e.g., Nick et al., 2009; Joughin et al., 2012a; Straneo and Heimbach, 2013; Jenkins et al., 2016). Increased melt and runoff driven by atmospheric warming is also a major component of mass loss in Greenland (e.g., Van Den Broeke et al., 2009).
Yet despite the clear signals of mass loss from the Greenland and Antarctic ice sheets (Shepherd et al., 2018), observed changes in terminus positions, velocities, and ice thickness vary widely among marineterminating outlet glaciers, especially in Greenland (e.g., Moon et al., 2014; Csatho et al., 2014). This makes it difficult to establish consistent links between forcing and response, which are critical for projections of future glacier and icesheet change. Some reasons for the heterogeneity in glacier change have become clearer. For example, localized features in bedrock and fjord geometry can control differences in terminus retreat (Catania et al., 2018) and inland dynamic thinning (Felikson et al., 2017) for individual glaciers. Additionally, regional variations may be associated with differences in surface melt and subglacial hydrology (Moon et al., 2014) or with the ocean waters with which the glaciers are in contact (Straneo et al., 2012). However, catchmentspecific factors continue to pose a challenge where observations are limited, as well as for icesheetwide simulations. Further observations will continue to be critical, but this heterogeneity between outlet glaciers is also a strong motivation to investigate general principles at the same time, which can be expected to apply to a wide range of settings.
In this study, we focus on a principle common to all outlet glaciers: climate forcing predominantly comes either from the atmosphere, via changes to the surface mass balance in the interior catchment, or from the ocean, via changes to ice discharge at the marine margin. Our approach is to conduct idealized model experiments that isolate key physical principles of transient glacier responses to climate. We use a recently developed reduced model (Robel et al., 2018) and an established numerical iceflow model (Pollard and Deconto, 2012) as complementary tools. The reduced model yields simple physical and mathematical interpretations of the numerical model's response to forcing. Additionally, the reduced model can efficiently generate ensembles of responses for statistical analyses.
We will focus exclusively on stable geometries. The marineicesheet instability (e.g., Weertman, 1974; Schoof, 2007a) is of course a critical consideration for some catchments, especially in West Antarctica. However, understanding how a physical system approaches, or fluctuates around, a longterm equilibrium is a core analytical approach. This requires that we start with stable systems, but the resulting insights can often be applied to unstable regimes (e.g., Robel et al., 2019). The following section describes experiments that illustrate the fundamental system dynamics and establish the key differences between ocean and interior forcing. We then present three cases that illustrate the implications of these principles for interpreting observations and predicting the future behavior of outlet glaciers.
2.1 A simple outletglacier system
Before describing the dynamic models, it is useful to begin with the geometry and basic fluxbalance arguments of the system we are investigating. We consider an idealized stable outlet glacier, a schematic of which is presented in Fig. 1a (see also Schoof, 2007a; Robel et al., 2018). Ice enters as snow accumulation over the interior surface, flows from the interior towards the ocean on a forwardsloping (prograde) bed, and exits the system where it reaches flotation at the grounding line. Beyond this point, floating ice is assumed to calve into icebergs or melt due to contact with the ocean.
Ice flux across the grounding line (Q_{g}) has long been known to be a sensitive function of local ice thickness (e.g., Weertman, 1974; Thomas, 1979; Lingle, 1984; Schoof, 2007a). This function can be represented in general form as
where h_{g} is ice thickness at the grounding line and Ω and β depend on ice dynamics near the grounding line. A floating ice shelf, a tongue, or friction from valley sidewalls typically provides buttressing. This can be modeled explicitly or analytically represented in Ω and β. Typically, β≥1, and this nonlinearity has been shown to depend on assumptions about the basal rheology near the grounding line (Schoof, 2007a; Tsai et al., 2015) and the characteristics of a buttressing ice shelf (Haseloff and Sergienko, 2018). h_{g} is simply determined by a flotation criterion: let b(x) be the bed elevation relative to sea level (i.e., negative at the grounding line). Then, for a glacier of length L, h_{g} is
where ρ_{w} and ρ_{i} are the densities of seawater and ice, respectively.
In this study, we consider the grounding line the boundary of the outletglacier system. That is, we consider floating ice a part of the boundary condition for the system's output flux (Q_{g}). In equilibrium, Q_{g} equals the surface mass balance integrated over the entire catchment (neglecting basal and englacial melt). Climate changes can perturb this flux balance in two distinct ways: the atmosphere can affect mass balance over the interior surface or the ocean can modulate ice discharge at the grounding line (e.g., via a change in buttressing). In either case, Eqs. (1) and (2) reveal the basic system response to an imbalance: on a prograde bed, the terminus retreats into shallower water to decrease ice flux out of the system or advances into deeper water to increase flux out of the system. The groundingline migration needed to restore equilibrium depends strongly on bed slope and the degree of nonlinearity in Q_{g} (i.e., via β).
2.2 Flowline model
To simulate the dynamics of this outletglacier system, we begin with a 1D (flowline) version of the icesheet model developed by Pollard and Deconto (2012; hereafter PD12). The evolution of local ice thickness h at a grid point reflects the balance of mass exchange at the surface and horizontal iceflux divergence. Conservation of mass requires that
where S is the local surface mass balance and $\stackrel{\mathrm{\u203e}}{u}$ is depthaveraged horizontal ice velocity. The velocity profile has contributions from stretching, shear, and basal sliding, which are summarized as follows. In general, a sloping ice surface $\frac{\partial s}{\partial x}$ creates a driving stress,
where ρ_{i} is ice density and g is acceleration due to gravity. Longitudinal stretching $\frac{\partial u}{\partial x}$ is represented by the shallowshelf approximation, where stretching and basal drag τ_{b} balance driving stress:
We assume a typical powerlaw rheology, with coefficient A and exponent n (e.g., Glen, 1955). Internal shear deformation ($\frac{\partial u}{\partial z}$) is represented by the shallowice approximation. At a depth h−z, this is
where ${\mathit{\tau}}_{\mathrm{e}}^{\mathrm{2}}=\left(\frac{\mathrm{1}}{\mathrm{2}}\right){\mathit{\tau}}_{ij}{\mathit{\tau}}_{ij}$ is a scalar effective stress and τ_{ij} is the deviatoric stress tensor. Finally, sliding velocity is given by a Weertmantype powerlaw relationship:
where C is a friction coefficient and m is the sliding exponent. Here, m is the inverse of the flow exponent n, but note that this convention is flipped in some texts. The PD12 model uses a combination of the shallowice and shallowshelf approximations to solve for $\frac{\partial u}{\partial z}$ and $\frac{\partial u}{\partial x}$, respectively. In this study, however, we use parameters that yield flow dominated by longitudinal stretching. We refer the reader to PD12 for further description of this hybridization.
A crucial simplification in the PD12 model is that flux across the grounding line is parameterized in the form of Eq. (1). In PD12 and this study, Ω and β are taken from the analytical solution of Schoof (2007a):
Θ is a buttressing factor between 0 and 1, and all other parameters are defined as above. The PD12 model calculates groundingline flux based on a thickness h_{g} that is linearly interpolated from the height above flotation of the lastgrounded and firstfloating grid cells to the point where flotation is reached. The corresponding subgrid groundingline position is shown for all output from the PD12 model in this study. Although the PD12 model is typically run on coarse grids (𝒪∼1–10 km) for continentscale simulations over many millennia, we use a grid of 100 m to better resolve the details of groundingline variations. Finally, while the PD12 model does simulate a floating ice shelf, its dynamics are not important to our analyses as its buttressing effect is parameterized via Eq. (8).
Figure 1b shows the steadystate profile for an idealized outlet glacier simulated with the flowline model. The domain begins at an ice divide and thus has a zeroflux lateral boundary condition. The bed is constant in time and has an elevation of −100 m at the divide and constant prograde slope b_{x} of $\mathrm{2}\times {\mathrm{10}}^{\mathrm{3}}$. Surface mass balance S is 0.5 m yr^{−1} ice equivalent, assumed to be spatially uniform. Additional parameters are given in Table 1. The glacier has an equilibrium length (ice divide to grounding line) of ∼185 km, a maximum thickness of ∼1580 m at the divide, and a thickness of ∼526 m at the grounding line. These scales are comparable to many outletglacier catchments in Greenland, though direct comparisons are limited as the 1D flowline cannot capture the flow convergence of many catchment geometries, and we emphasize that we are not simulating a particular outlet glacier. We discuss two additional geometries for comparison in Sect. 3.1.
2.3 Response to forcing
We begin by comparing the flowline model's response to forcing from either surfacemassbalance changes in the interior or ocean forcing at the terminus. In all model experiments throughout this study, “interior forcing” and “ocean forcing” will refer to perturbations applied in the following manner. Interior surfacemassbalance anomalies are assumed to be spatially uniform. We represent ocean forcing very simply by perturbing the groundinglineflux coefficient, Ω. This broadly represents changes in the buttressing provided by an ice shelf or fjord walls, driven by anomalous melting or calving. In principle this could be targeted via Θ in Eq. (8), and other analytical formulations for Ω explicitly represent a buttressing ice shelf (Haseloff and Sergienko, 2018). Perturbing these parameters might be more realistic, but we focus on Ω so that we can very generally represent flux perturbations, which may in reality result from a host of ice–ocean interactions. Our primary interest is in how glacier dynamics respond to each forcing type, and representing ocean forcing in this way makes comparison straightforward. For example, a fractional change in S (surface mass balance) constitutes a flux anomaly with the same initial magnitude as the same fractional change in Ω.
Figure 1c shows the length response to step forcings at time t=0. The step is a 20 % decrease in S or a 20 % increase in Ω (i.e., more discharge for a given h_{g}) and, in both cases, forces the terminus to retreat into shallower water to reduce the output flux. Note that the responses to changes in S and Ω do not converge to the same final length. The equilibrium sensitivities differ because of different nonlinearities in the input and output fluxes (S⋅L and Q_{g}) as the system state evolves. However, our main focus is the marked difference between the initial responses to the step change, highlighted in the inset panel of Fig. 1c. Perturbing Ω results in a much faster initial retreat (blue curve) compared to perturbing S (orange curve), despite the larger equilibrium sensitivity to S. However, this faster retreat rate lasts only the first century or so and accommodates only ∼25 % of the total equilibrium response. The remaining retreat occurs at a rate similar to that driven by a change in S, which is an asymptotic approach to equilibrium over several millennia. Thus, the length response to Ω forcing has both a “fast” and “slow” component, whereas S forcing primarily drives a slow response.
As an alternative to an impulse forcing, we also consider the glacier's response to stochastic variability in either Ω or S. We apply this variability as interannual white noise, which by definition has equal power at all frequencies and no interannual persistence. We apply the exact same whitenoise time series as either Ω or S anomalies but with the opposite sign so that the corresponding icevolume anomalies match. The anomaly time series (hereafter Ω^{′} and S^{′}) are scaled to have standard deviations equal to 20 % of the mean values ($\stackrel{\mathrm{\u203e}}{\mathrm{\Omega}}$ and $\stackrel{\mathrm{\u203e}}{S}$). Figure 1d shows the resulting length responses. For both types of forcing, the glacier acts as a lowpass filter on the imposed climate anomalies, producing kilometerscale fluctuations with clear persistence. However, the length anomalies driven by Ω^{′} have much greater highfrequency content, and greater overall variance, than those driven by S^{′}. The highfrequency response is perhaps intuitive, in that forcing applied at the grounding line has an immediate effect on groundingline position. However, the response to Ω^{′} also contains millennialscale fluctuations onto which the faster variations are superimposed. These slow, wandering excursions are comparable to those driven by S^{′} and, like the multimillennial adjustment to step changes (Fig. 1c), suggest slow dynamics common to both responses.
Variability is intrinsic to climate, arising from the fundamentally chaotic nature of the natural system. The associated glacier fluctuations are a crucial part of characterizing glacier dynamics, and the implications have been extensively explored for mountain glaciers (e.g., Oerlemans, 2001; Roe, 2011; Roe et al., 2017) and more recently for marineterminating outlet glaciers (Robel et al., 2018) and ice streams (Mantelli et al., 2016). Figure 1 suggests a new and intriguing complication for outlet glaciers: the timescale and magnitude of glacier fluctuations depend on whether the climate variability comes in the form of surfacemassbalance or ocean anomalies. Additionally, the response to step changes (Fig. 1c) suggests that the type of forcing is also relevant for the response to nonstationary climate changes. In the next section, we investigate these contrasting responses using a recently developed reduced model.
2.4 The twostage model
Robel, Roe, and Haseloff (2018; hereafter RRH) developed a simple model for marineterminating outletglacier dynamics, derived from iceflux conservation and constrained by largescale glacier geometry. RRH showed that this model could accurately emulate a more complex numerical model on the multidecadal and longer timescales on which most glacier variance occurs. A full description and derivation can be found in RRH, but a summary of the model follows.
The RRH model reduces the outlet glacier to a system with two degrees of freedom: L, the length from ice divide to grounding line, and H, the characteristic interior ice thickness (Fig. 2a). The glacier can be conceptualized as a small groundingzone reservoir with a length l_{gz}≪L and thickness h_{g}, coupled to a large interior reservoir with thickness H and length L−l_{gz}. The geometry is further described by a bed topography with constant average slope b_{x}. Recall that the ice thickness at the grounding line, h_{g}, is a statedependent quantity, set entirely by L and b(x) (see Eq. 2).
The model dynamics are derived by balancing ice fluxes through these linked reservoirs (Fig. 2a). Ice thus enters the system via an accumulation flux S⋅L, flows via an interior flux Q to the grounding zone, and leaves the system as a flux across the grounding line Q_{g}. A steady state of H and L is one which balances these three fluxes. L ultimately controls the system's input and output fluxes by setting the total catchment area for accumulation and by controlling Q_{g} via h_{g} (Eqs. 1 and 2). Interior dynamic ice flux has the form
where α and γ can vary depending on the particular processes controlling interior flow. We use $\mathit{\alpha}=\mathrm{2}n+\mathrm{1}=\mathrm{7}$ and $\mathit{\gamma}=n=\mathrm{3}$, consistent with deformation dominated by longitudinal stretching (RRH).
Two coupled equations capture the transient adjustment of the two degrees of freedom, H and L, as they relax towards a steady state that balances all three fluxes:
Because achieving a steady state requires adjustment of both H and L, we refer to this model as the “twostage model”, following RRH. We discuss the operation of these stages further in Sect. 2.6. Note that Eq. (10) captures the nonlinear dependence of ice flux on driving stress (e.g., Eqs. 4–6) and ice thickness and that the first two terms of Eq. (11) resemble the continuity equation that governs ice thickness in the flowline model (Eq. 3). The twostage model makes significant simplifications but stems from the same physical principles as the flowline model and most other contemporary icesheet models (see RRH for further discussion).
RRH also linearized these equations for fluctuations H^{′} and L^{′} about a steady state, $\stackrel{\mathrm{\u203e}}{H}$ and $\stackrel{\mathrm{\u203e}}{L}$. The linear model yields analytical expressions for a number of important quantities, including two widely separated characteristic response timescales. With some simplifications detailed in RRH, the fast and slow response timescales (τ_{F} and τ_{S}, respectively) are
where ${s}_{\mathrm{T}}=\mathrm{1}+\frac{{\mathit{\rho}}_{\mathrm{w}}}{{\mathit{\rho}}_{\mathrm{i}}}\frac{\mathit{\beta}\stackrel{\mathrm{\u203e}}{{b}_{x}}\stackrel{\mathrm{\u203e}}{L}}{\stackrel{\mathrm{\u203e}}{{h}_{\mathrm{g}}}}$ is a parameter describing the stability of glacier response based on its geometry and groundingline dynamics (via β). Both response times contain a characteristic thickness divided by mass balance rate $\stackrel{\mathrm{\u203e}}{S}$, reminiscent of the canonical mountain glacier response time (Jóhannesson et al., 1989). τ_{F} is controlled by h_{g} and thus the volume of the system's small groundingzone reservoir, whereas τ_{S} relates to the volume of the interior reservoir via $\stackrel{\mathrm{\u203e}}{H}$. For scales typical of outlet glaciers, τ_{F} is on the order of decades to centuries, whereas τ_{S} is on the order of millennia. Thus, the twostage model approximates the outletglacier system as linked interior and groundingzone reservoirs with distinct timescales.
Here we generalize the RRH linearization to simultaneously include perturbations to interior surface mass balance (S^{′}) and groundingline flux (${Q}_{\mathrm{g}}^{\prime}$; proportional to Ω^{′}). The linearized equations take the form of a twodimensional dynamical system with two forcing terms:
A_{H}, A_{L}, B_{H}, and B_{L} are shorthand for the couplings between length, thickness, and flux changes, which are given in the Appendix. This linear system is readily discretized, and can be cast as a twodimensional autoregressive process (e.g., Box et al., 2008), which we also present in the Appendix.
2.5 Model comparison
RRH showed that the twostage model emulated the response of a flowline model forced with surfacemassbalance anomalies. Their flowline model (described in Schoof, 2007b; Robel et al., 2014) used a stressbased, as opposed to fluxbased, groundingline condition but was otherwise dynamically similar to the PD12 model. Here, we use the PD12 model and extend the comparison to groundinglineflux perturbations. Although the twostage model output includes H, we focus our comparison on L anomalies, which are more directly comparable between models. H is a more heuristic variable in the twostage model; however, it does govern the evolution of interior fluxes, which we discuss in the next section. Both models use the parameters given in Table 1, for which Eqs. (13) and (14) give response times of ∼76 and ∼2000 years, respectively.
Figure 2 shows output from both models in response to step and stochastic forcings. Figure 2b shows the groundingline retreat following a 20 % increase in Ω (blues) and 20 % decrease in S (orange and brown). For clarity, only the linearized twostage output is shown. Note that the nonlinear twostage model is constrained to have the same equilibrium response as the flowline model by Eq. (1). The linear and nonlinear responses match almost exactly for the stochastic fluctuations but disagree somewhat for the step changes. However, the disagreement is not severe, suggesting that we can reasonably use the linearized model and its response times as analytical tools. Importantly, the twostage model captures the faster initial response to forcing at the grounding line, with a slightly more pronounced slope break in the first few hundred years.
Figure 2c shows a 5 kyr sample of 100 kyr of length fluctuations driven by whitenoise forcing in S^{′} and Ω^{′}. Again, anomalies have no interannual persistence, are identical except for a sign reversal between Ω^{′} and S^{′}, and have a standard deviation of 20 % of $\stackrel{\mathrm{\u203e}}{\mathrm{\Omega}}$ or $\stackrel{\mathrm{\u203e}}{S}$. The twostage and flowline models agree closely for both types of forcing, although the twostage model slightly underestimates the magnitude of highfrequency fluctuations. Figure 2d shows the autocorrelation function of each length time series. This is a measure of the persistence of anomalies and reveals the memory within a dynamical system. For forcing from Ω^{′}, the autocorrelation drops off rapidly to approximately 0.4 after a couple of centuries, corresponding to the fast response time. However, it then requires several millennia to drop below ∼0.1, indicating the persistent influence of the slow timescale. For S^{′} forcing, the autocorrelation is dominated by the slow timescale. The power spectra of length fluctuations are shown in Fig. 2e. Because whitenoise forcing has a flat power spectrum, the shapes of the response spectra reveal the glacier's filtering properties. The glacier is a strong lowpass filter for both forcing types but damps high frequencies even more for S^{′} forcing. The twostage model underestimates the highfrequency response of the flowline model, but agreement is very good at the lower frequencies that contain the majority of the variance. Critically, both models capture the clear split between the spectra of Ω^{′} and S^{′}driven anomalies at frequencies of around 10^{−2} to 10^{−3} yr^{−1}, where the models agree quite well.
2.6 Interpretations
The twostage model's agreement with the flowline model suggests that it captures the essential dynamics responsible for the contrasting transient responses to interior and groundinglineflux anomalies. At this point, it is useful to discuss some of the interpretations enabled by this reduced model.
The linearized twostage model (Eq. 15) is a dynamical system with two eigenmodes, which, for the stable geometries considered here, are exponentially decaying functions with characteristic times τ_{F} and τ_{S} (the eigenvalues being ${\mathit{\tau}}_{\mathrm{F}}^{\mathrm{1}}$ and ${\mathit{\tau}}_{\mathrm{S}}^{\mathrm{1}}$). The responses of the system's two degrees of freedom (H^{′} and L^{′}) are linear combinations of these fast and slow eigenmodes. However, the two forcing types, Ω^{′} and S^{′}, project onto each mode with different relative weights. Forcing via S^{′} projects almost entirely onto the slow mode. For forcing in Ω^{′}, the fast mode makes a substantial contribution on multidecadal to centennial timescales. However, the slow mode still dominates the equilibrium response to an impulse, as well as the power spectrum of stochastic fluctuations.
These modes can also be conceptualized as a twostage lowpass filter on any forcing time series (e.g., Fig. 2e), and the type of forcing determines the order in which they are applied. Forcing over the interior (S^{′}) is first filtered by the slow mode, and the fast mode then operates on interior flux anomalies whose highfrequency content has already been strongly damped; there is little additional filtering to do. In contrast, if forcing is applied at the grounding line (${Q}_{\mathrm{g}}^{\prime}$ via Ω^{′}), the fast mode responds to unfiltered anomalies and the groundingline position exhibits a fast initial response before engaging the slow mode in the interior. Figures 1 and 2 suggest that the response to S^{′} could be reasonably approximated by the adjustment of a single multimillennial mode. However, we stress that the response to Ω^{′} is not simply the response of a single multidecadal mode, because the grounding zone is coupled to the slower but ultimately more sensitive interior reservoir.
While these mathematical interpretations may seem abstract, it is helpful to remember that the twostage model was derived from mass conservation, and that the linear response times (eigenmodes) reflect the largescale glacier geometry. Because a glacier is a Stokes (i.e., noninertial) fluid driven by potentialenergy gradients, the largescale glacier geometry must reflect the aggregate dynamics by which the system seeks flux balance. This relationship between geometry and dynamics is another way to interpret why Ω^{′} and S^{′} forcings project differently onto the fast and slow modes: these flux imbalances have different geometries, and thus the geometry and dynamics of the transient responses must also be distinct.
To help illustrate this, Fig. 3 tracks the evolution of fluxes following step changes in S and Ω. Five fluxes are shown for the flowline model: over the interior surface (S⋅L); across the grounding line (Q_{g}); and at conceptual flux gates located 5, 25, and 50 km upglacier of the grounding line (Fig. 3a). For the nonlinear twostage model, S⋅L, Q, and Q_{g} are shown (Fig. 3b). In all cases, flux anomalies are normalized by their final (equilibrium) change. A change in surface mass balance has an immediate, spatially distributed effect on interior thickness, which slowly alters driving stresses and thus fluxes throughout the interior. Anomalous ice flux arrives at the grounding zone, driving advance or retreat, which then modifies Q_{g} according to Eq. (1). Most of the disequilibrium during the transient response is between surfacemassbalance and interior fluxes, while Q_{g} can quickly keep pace with flux anomalies from the interior. In the twostage model, Q_{g} keeps up on the fast timescale. In contrast, a perturbation to Q_{g} (here, via Ω) is highly localized and first creates a large disequilibrium between Q_{g} and interior flux. The increase in Q_{g} causes the grounding zone to drain, drawing in ice from the interior, which transfers the disequilibrium from Q and Q_{g} to Q and S⋅L. Both models capture this transfer, and the flowline model flux gates show that it gradually propagates up from the grounding zone. The transient peak in fluxes is similar to a kinematic wave propagating from the terminus, which reaches well into the interior within a few multiples of τ_{F}. The ensuing drawdown of the interior reservoir (the slow mode in the twostage model) brings interior fluxes back down, and again, Q_{g} must follow via groundingline retreat.
Although the flowline model captures a more realistic and spatially distributed response, the twostage model contains enough geometric information to emulate the basic sequence of flux anomalies. Note that more discretized interior reservoirs could be added to the twostage model, which would eventually approach the form of the flowline model. However, a single reservoir could not capture the essential transient response. The distinct stages of the response to a perturbation in Q_{g}, borne out in both models, show that a twomode framework is useful for understanding the response to forcing applied at the grounding zone.
2.7 Persistence in variability
In the previous sections, we imposed stochastic variability with a flat power spectrum (i.e., white noise) because it allowed us to identify the influence of forcing type and model physics across all relevant frequencies. In reality, surfacemassbalance and ocean variability may exhibit different power spectra. Climate variables associated with the atmosphere (such as surface mass balance) often have little interannual memory (e.g., Medwedeff and Roe, 2017), but the ocean's thermal inertia can introduce persistence (e.g., Hasselmann, 1976). While the comparisons above show that interior and ocean forcing elicit different glacier responses due to the distinct geometries of the respective flux anomalies, we must also consider the possibility of climatic persistence in order to fully characterize an outlet glacier's response to natural climate variability. Roe and Baker (2016) showed that persistence in surfacemassbalance variability amplified the length fluctuations of mountain glaciers; we can expect similar principles to be especially relevant to outlet glaciers if they are sensitive to decadescale ocean variability.
We consider four types of synthetic forcing with different persistence characteristics: (1) white noise; (2) and (3) firstorder autoregressive (AR1) processes with persistence timescales τ_{AR1} of 4 and 20 yr, respectively; and (4) powerlaw persistence. These are all plausible models for the sort of persistence that may affect outletglacier forcing. For example, an AR1 process with a 4 yr memory is characteristic of seasurfacetemperature anomalies in the North Atlantic (e.g., Deser et al., 2003), while the 20 yr timescale would better correspond to subsurface anomalies (e.g., Gwyther et al., 2018; Jenkins et al., 2018). Additionally, the continuum of surfacetemperature variability from interannual to multimillennial timescales has been described with a powerlaw spectrum (e.g., Huybers and Curry, 2006).
We generate time series following Percival et al. (2001) and Roe and Baker (2016), first defining the desired spectral properties in frequency space and then applying the same set of random phases to each case so that the resulting anomalies are correlated in the time domain. The power spectrum of a discrete AR1 process as a function of frequency, f, is
where P_{0} scales the total variance and r is the autocorrelation at a lag of Δt (here, 1 yr). The memory timescale of the process is related to r by ${\mathit{\tau}}_{\mathrm{AR}\mathrm{1}}=\frac{\mathrm{\Delta}t}{\mathrm{1}r}$. Powerlaw noise has a spectrum defined by
where f_{0} is the highest sampled frequency. The power increases continually out to low frequencies, where ν is the slope of the spectrum in log–log space. We use ν=0.5, consistent with ν∼0.5 to 1 identified by Huybers and Curry (2006) in a large collection of paleoclimate records.
We normalize the forcing time series so that they all have the same variance. This ensures that, for a given choice of S^{′} or Ω^{′} forcing, differences in glacier response are due only to differences in persistence. Figure 4a shows 200 yr samples of the 10^{5} yr synthetic forcing time series, and Fig. 4c shows their power spectra. The AR1 processes (red, maroon) have clear persistence in the time domain. The lowfrequency content of powerlaw noise (blue) is hard to discern at short timescales, but the spectra show that it has similar power to the AR1 processes at millennial timescales (and, again, identical total variance). Powerlaw variability is thus difficult to constrain from centuryscale instrumental records alone (e.g., Percival et al., 2001). At a certain point, the spectrum of climate variability runs into timescales where paleoclimate reconstructions, rather than synthetic noise or instrumental observations, would be a more relevant choice for understanding natural glacier variability; we will return to this point in Sect. 3.3.
Figure 4b shows 10^{3} yr samples of the resulting length responses for the twostage model. The effects of persistence are dramatic: for both S^{′} or Ω^{′} forcing, 4 yr AR1 persistence more than doubles σ_{L} compared to length fluctuations driven by whitenoise forcing, and both 20 yr AR1 and powerlaw persistence cause a ∼5to6fold increase in σ_{L}. This is consistent with the approximate proportionality between ${\mathit{\sigma}}_{\mathrm{L}}^{\mathrm{2}}$ and τ_{AR1} predicted by theory (see, e.g., Roe and Baker, 2016; Robel et al., 2019).
The power spectra of the forcings and responses show why persistence has such a strong effect. The response to whitenoise forcing reveals the system's sensitivity across all timescales (Fig. 4d), because the input has equal power at all frequencies. The glacier response is set by its internal dynamics and the forcing type (S^{′} or Ω^{′}). Any forcing time series must be filtered by the same dynamics, and if the forcing has persistence, some of its spectral power is shifted towards lower frequencies where glacier sensitivity is higher. Thus, the shape of the response spectrum (Fig. 4e) reflects both the frequency content of the forcing (Fig. 4c) and the frequency dependence of the glacier dynamics or, in other words, the spectral shape of the glacier filter (Fig. 4d).
The practical takeaway is that persistence in climate forcing increases the total variance of length fluctuations. When combined with the finding that terminusflux anomalies excite the fast mode of response more than surfacemassbalance anomalies, the implication is that ocean variability – which tends to have persistence – may drive much larger terminus fluctuations than surfacemassbalance variability.
Figures 1–3 demonstrate that marineterminating outlet glaciers have different transient responses to interior and ocean forcing, because of how the fast and slow modes respond in each case. In the next sections, we examine the consequences for three key issues: the committed response to forcing, the attribution of an observed change, and a glacier's memory of past climate changes. The relative roles of ocean and interior forcing will, of course, vary widely among individual glaciers. Rather than conducting simulations of specific settings, we will use our simplified model framework to outline the general implications and assess how the combination of fast and slow dynamics applies to each question.
3.1 Committed change
Any system with a nonzero response time will lag applied forcing. “Committed change” refers to the total response such a system would need to undergo to attain equilibrium with the current level of forcing. In the context of a warming climate, committed change is an important lower bound on future change that is independent of future emission scenarios. It has long been recognized that surface temperatures lag CO_{2} forcing due to the ocean's thermal inertia (e.g., Hansen et al., 1985; Wigley and Schlesinger, 1985). In turn, other aspects of the climate system respond to warming with their own lags, including mountain glaciers (e.g., Christian et al., 2018) and global sea level (e.g., Meehl et al., 2005; Levermann et al., 2013). Industrialera climate forcing began in the 19th century (e.g., IPCC, 2013; Abram et al., 2016), and so the millennial response times (i.e., Eq. 14) of outlet glaciers imply a severe disequilibrium with current climate and a large committed change, even with no additional forcing.
To illustrate the current committed change in outlet glaciers, we use the twostage model and follow the framework of Christian et al. (2018), using idealized glacier geometries and forcings. Disequilibrium, and thus committed retreat, is defined in reference to an “instantaneous equilibrium” response, which is the state at which the glacier would be in equilibrium with climate at any given time. This is governed largely by the bed geometry and climate forcing. An idealized geometry (e.g., uniform bed slope and width) makes for a simple equilibrium response, but the physical principle of course also holds for more complex systems. For an idealized outlet glacier, the instantaneous equilibrium is the groundingline position (L) that yields ${Q}_{\mathrm{g}}=S\cdot L$.
We consider ramp forcing scenarios where S decreases, or Ω increases, by 30 % from 1880 to 2020 CE. The forcing is fixed at 2020, and the glacier is allowed to relax towards a new equilibrium. Figure 5 shows the responses of three glaciers with varied response times: (1) ${\mathit{\tau}}_{\mathrm{F}},{\mathit{\tau}}_{\mathrm{S}}\sim (\mathrm{76},\mathrm{2000})$ yr; (2) ${\mathit{\tau}}_{\mathrm{F}},{\mathit{\tau}}_{\mathrm{S}}\sim (\mathrm{56},\mathrm{1200})$ yr; and (3) ${\mathit{\tau}}_{\mathrm{F}},{\mathit{\tau}}_{\mathrm{S}}\sim (\mathrm{140},\mathrm{4600})$ yr. The parameters for each glacier are provided in Table 2. Dashed lines show the instantaneous equilibrium responses, which indicate the total committed response at any given time.
The most basic result is that, in all cases, the transient response as of 2020 is a small fraction (∼1 %–23 %) of the instantaneous equilibrium response. The disequilibrium for S forcing is particularly severe, as pointed out previously by RRH. Although changes in surface mass balance are already a major component of overall mass loss in some settings, basic glacier dynamics require that S anomalies also have an impact on the groundingline position, which is, as yet, essentially unrealized. Perturbations in Ω elicit a faster response over the industrial era, but Fig. 5a–c show that both types of forcing have a longterm commitment associated with the slow drawdown of interior ice.
The slow response (interior thinning → reduced interior fluxes → groundingline retreat) comprises the majority of the response to S forcing, but it is important to remember that it is also the second stage of the response to ocean forcing, and it dictates the total committed change. Following recent observed increases in ice discharge, Price et al. (2011) and Goldberg et al. (2015) examined the committed responses of several outlet glaciers in Greenland and Antarctica, respectively. These studies found substantial “dynamic” commitments associated with thinning and elevated fluxes, although they focused on nearterm sealevel rise and did not project past 2100. Our fluxbalance perspective is a complementary view of committed change and shows that dynamic thinning necessarily drives additional retreat on long timescales as a consequence of reduced interior fluxes. The slow equilibration of the interior would thus be a critical part of determining the total commitment due to ocean forcing, especially if the longterm retreat moves the terminus to more or lessstable slopes.
The slow mode also means that there can be large uncertainties in committed change if the magnitude of forcing is uncertain. Figure 5d shows the response to a trend in Ω, ranging from a 20 % to 40 % increase from the initial value. As of 2020, the differences in observed retreat are small but diverge as the system approaches equilibrium over subsequent millennia. In other words, the slow response limits how “wrong” shortterm simulations might be due to errors in the assumed forcing but makes no such constraints on the committed change. The same principle also applies to uncertainty in the outlet glacier's length sensitivity: note that glaciers 1 and 2 have nearly identical responses on centennial timescales but different equilibrium sensitivities. Similar arguments have been made for uncertainty in radiative forcing and equilibrium climate sensitivity (Armour and Roe, 2011).
3.2 The emergence and detectability of forced responses
We now turn to the topic of attributing outletglacier retreats to natural or anthropogenic forcing. The attribution of an observed change to a particular cause (i.e., an external forcing) can be a challenge because of factors specific to individual glaciers, such as complexities in bed geometry, regional climate, or the local collection of ice–ocean interactions. It can also be a challenge because of factors intrinsic to the transient ice dynamics, which affect the amount of the forced response that can be expressed over a given time. We focus here on this latter set, and in particular on the contrasting implications of ocean vs. interior forcing.
Attribution is often framed in terms of the signaltonoise ratio (SNR) of an observed trend. For a variable x, the SNR is often defined as $\frac{\mathrm{\Delta}x}{{\mathit{\sigma}}_{x}}$, where Δx is the change estimated by a linear fit and σ_{x} is the standard deviation of the residuals. For glacier retreat (or changes in any analogous system), this depends on the SNR of the relevant forcing variable (e.g., surface temperature, ocean heat content) and also on the memory timescale over which the glacier integrates anomalies in that variable. A glacier's memory damps its response to highfrequency noise, but the tradeoff is that it also delays its response to a persistent trend (see Fig. 5). Roe et al. (2017) showed that mountain glaciers exemplify this concept and, moreover, that their multidecadal response times are optimal for damping interannual climate variability while responding sensitively to centennial trends, thereby producing an amplified SNR.
For outlet glaciers, however, two types of forcing (ocean and interior) and much longer response times are at play. To understand how these factors interact on the timescales of historical anthropogenic forcing, we consider three idealized scenarios of stationary variability plus a trend (Fig. 6a). Cases 1 and 2 have variability with no interannual persistence (i.e., white noise) in S and Ω, respectively. Case 3 has multidecadal variability in Ω, generated by an AR1 process with a memory timescale of 20 yr (see Sect. 2.7). To compare glacier responses, we standardize the forcing scenarios by their SNR: in each case, a linear trend begins in 1880 CE (negative in S for case 1, positive in Ω for 2 and 3). We stipulate that the trend reaches a 20 % departure from the mean by 2020 CE and continues thereafter at the same rate. The stationary variability again has σ=20 % of the initial mean. Thus, the imposed trends in climate forcing all have a SNR that reaches 1 by 2020.
Figure 6b shows the resulting groundingline anomalies of the test glacier described in Sect. 2, generated with the twostage model. In all cases, the simulation was initialized in the year 0 CE. For reference, responses without variability are also shown (thin lines), as well as responses with variability but no trend. Shading indicates the 1σ and 2σ bounds of groundingline fluctuations, determined from 10^{7} yr simulations of stationary variability.
In case 1, noisy surface mass balance drives small length fluctuations, but the very slow response to the trend in S means that the forced response is slow to emerge from the noise; it remains within 2σ_{L} until the late 21st century. In contrast, the faster response to a trend in Ω exceeds 5σ_{L} by 2020, for case 2 (interannual variability). The multidecadal fast response time thus acts to amplify the climate trend's SNR roughly by a factor of 5, consistent with the findings of Roe et al. (2017) for mountain glaciers. However, because persistence amplifies natural glacier variability (Sect. 2.7), the multidecadal variability in case 3 inflates the envelope of natural glacier fluctuations such that the forced response is again hard to detect on centennial timescales. The forced response is roughly 1σ_{L} by 2020, meaning that glacier dynamics no longer improve upon the SNR of the forcing. Thus, a multidecadal fast timescale does not necessarily amplify the SNR of a climate trend if a glacier is subjected to multidecadal ocean variability.
These idealized cases are useful for understanding factors that affect the detectability of a length trend, but they also demonstrate that the SNR may be a problematic metric in a practical sense. Because most of the natural glacier variability is expressed at very low frequencies, σ_{L} (i.e., the noise) is undersampled by direct observations, which cover only several decades for most outlet glaciers. Additionally, slow natural excursions could dominate the baseline against which a trend is measured. For example, in case 3, the initial length anomaly due to random variability is greater in magnitude in 1880 than the forced component of the retreat from 1880 to 2020 (i.e., the difference between orange and blue lines).
An alternative approach is to evaluate an observed ΔL against the distribution of natural trends that occur over the same interval of time as the observation. This incorporates information about the rate and persistence of natural length fluctuations. The simplicity of the twostage model allows us to generate large synthetic ensembles of random fluctuations for a range of glacier parameters and forcing scenarios. We develop distributions of naturally forced trends as follows: if Δt_{obs} is the length of a hypothetical observational period, we draw the last Δt_{obs} years from model runs of 10^{4} years, each of which are forced with stationary noise and no trend. For each realization, ΔL_{null} is the trend in groundingline position over the observational period. We find ΔL_{null} for 10^{4} independent simulations for each type of variability, yielding distributions of ΔL_{null} for a given Δt_{obs}. The resulting distributions are shown in Fig. 6c, for Δt_{obs} of 50, 100, 250, and 500 years.
For variability in S, the distribution of ΔL_{null} is narrow for short Δt_{obs} due to the strong damping of high frequencies but widens as Δt_{obs} increases. In other words, large fluctuations are possible, but ice dynamics constrain them to be slow. As expected, Ω variability increases the likelihood of observing larger trends on centennial timescales, and the effects of persistence are once again dramatic: a 1 km retreat in 50 yr is very rare with no persistence (99th percentile) but fairly commonplace with multidecadal variability (70th percentile).
As Δt_{obs} increases out to 500 years, the distributions of ΔL_{null} do not widen as much for variability in Ω as for variability in S. This suggests that centennial timescales (i.e., a few multiples of τ_{F}) are nearly optimal for sampling large, persistent trends driven by stochastic Ω variability alone (note that in the limit of Δt_{obs}≫τ_{S} , the distributions would narrow back to zero). This raises the importance of the fast mode for attributing changes in the time frame of one to two centuries of anthropogenic climate forcing.
The widely separated dynamical timescales characteristic of outlet glaciers pose a unique challenge for understanding their modern changes. The slow mode means that the overall variance (i.e., ${\mathit{\sigma}}_{\mathrm{L}}^{\mathrm{2}}$) must be considered, as stochastic variability over the last millennia may be important for the preindustrial state. Yet if the glacier is sensitive to variability at the terminus, the fast mode enables shorterterm fluctuations that may obscure the early response to anthropogenic forcing, especially if the variability has significant persistence. The sensitivity of groundingline flux to ocean variability and large calving events, as well as the associated timescales, will thus be critical for attribution studies. Finally, it should be borne in mind that these detection challenges arise from slow dynamics and not from low sensitivity. As is clear from Fig. 5, only a small fraction of the committed response is available for attribution studies today.
3.3 Inherited conditions
We have thus far considered outletglacier fluctuations due to stationary interannual to multidecadal climate variability. However, long response times also imply some memory of climate variations throughout the Holocene. Icesheet models are often spun up using paleoclimate proxy data precisely for this reason (e.g., Bindschadler et al., 2013). However, these strategies do not always reproduce the same modern ice extent (e.g., Goelzer et al., 2018), and the simulated icesheet history can depend strongly on the choice of proxy and its implementation in the model (e.g., Nielsen et al., 2018; Buizert et al., 2018). Here, we compare outletglacier responses to ocean versus interior forcing over the Holocene as an additional factor for the modern state. We focus on climate anomalies that are welldocumented in the Northern Hemisphere, but similar considerations would apply to Antarctic outlet glaciers.
First, we consider a climate scenario with idealized representations of three events: a deglacial warming at 11 kyr, a Little Ice Age (LIA) cool period from 1450 to 1850 CE, and an anthropogenic warming trend from 1880 to 2100 CE. The deglacial and LIA transitions are smoothed by error functions of 500 and 50 yr widths, respectively. The magnitudes of deglacial, LIA, and anthropogenic events have a ratio of $\mathrm{10}:\mathrm{1}:\mathrm{4}$, and the intervening Holocene climate is assumed constant. We scale forcings linearly to the climate anomalies, where warming corresponds to negative S^{′} or positive Ω^{′}. Obviously, different combinations of these forcings could be expected in reality. Rather than choosing a particular combination, we examine each in isolation and normalize the glacier responses. These experiments thus serve as limiting cases to illustrate the relative influences of ocean and interior forcing over different timescales. This scenario is designed to explore two practical points: (1) the glacier's memory of large, longago events compared to smaller, more recent events and (2) the glacier's relative memory of past ocean vs. interior forcing.
We use the twostage model, linearized with respect to the Holocene climate with parameters for glacier 1 (see Table 2; ${\mathit{\tau}}_{\mathrm{F}},{\mathit{\tau}}_{\mathrm{S}}\sim \mathrm{76}$, 2000 yrs). The advantage of using the linear model here is that it has uniform sensitivity to fractional perturbations in S and Ω. This allows us to more clearly distinguish the signatures of fast and slow dynamics; the tradeoff is that it ignores nonlinearities that are surely a factor for very large climate changes. Accurate simulations over such transitions would depend not only on nonlinearities in ice dynamics but also on spatial information (e.g., bed topography) that is simplified in reduced models. Nevertheless, the linear model is a straightforward tool for demonstrating consequences of having both fast and slow dynamics. Focus should be directed towards the responses to the idealized LIA, for which the linearization is more valid. Including the deglaciation signal primarily serves to account for residual, albeit faint, disequilibrium implied by millennial response times.
Figure 7a shows the idealized climate (top) and twostage model responses. Anomalies are shown relative to the midHolocene and normalized to the large deglacial transition. Figure 7b shows 1000 to 2000 CE in more detail, including a scenario with no anthropogenic warming (dashed). For reference, the length at which the glacier would be in equilibrium with the LIA climate is also plotted (gray). For both forcings, the groundingline retreat due to the deglacial signal is nearly complete by the onset of the LIA. However, the LIA advance is $\sim \mathrm{2}\times $ greater for forcing in Ω, because it can engage the fast mode to a greater degree. Yet, even forcing in Ω yields a muted response: the transient response only reaches ∼35 % equilibration before the period ends. The slow mode, which takes up the majority of the response for both forcing types, barely feels our 400 yr LIA before it is reversed. This is worth bearing in mind whenever the duration of glacier excursions and climate anomalies are less than τ_{S}, because the system never achieves equilibrium. This would be an issue particularly if such events are used to tune glacier sensitivity in models.
The idealized LIA is a much more discrete “event” than is supported by paleoclimate records and ignores other variations in the Holocene. Thus, we also consider the response to a more realistic forcing time series. Again, this is not a reconstruction of actual terminus changes, but it is useful to see how glacier memory integrates the continuum of variations found in paleoclimate records. We use a time series of temperatures for Disko Bay, Greenland, from the regional reconstruction of Buizert et al. (2018) (Fig. 7c–d). The forcings are scaled as above and normalized to 1880 CE.
The glacier responses in Fig. 7d show the same essential behavior as the more idealized case in Fig. 7b. Even with a more subtle LIA, the terminus response is much more pronounced if it is driven by Ω anomalies that engage the fast response. It is worth noting, though, that most of the glacier disequilibrium in the 1800s is due to cooling over the previous millennia, to which the slow mode responds with a pronounced lag. Thus, the preindustrial state of an outlet glacier may depend significantly on both LIA and longerterm ocean cooling, via its two distinct response timescales.
Marineterminating outlet glaciers are sensitive to two fundamental types of forcing: changes in surface mass balance, which are distributed over a large interior catchment, and changes in their flux at the grounding line, which are typically driven by the ocean. We have used two models of different complexity to explore and contrast the dynamic responses to these two categories of forcing. Our key findings are as follows:

Ocean forcing (via the groundinglineflux coefficient, Ω) and interior surfacemassbalance forcing elicit fundamentally different transient groundingline responses (Figs. 1 and 2). The response to ocean forcing is characterized by a fast initial groundingline migration, followed by a second, much slower stage of adjustment. In contrast, the response to interior forcing is dominated by slow groundingline migration, with very strong damping at short timescales.

The twostage model (Robel et al., 2018) captures the evolution of flux anomalies that gives rise to these two contrasting dynamic responses (Fig. 3). This lends further confidence to the twomode interpretation of outletglacier dynamics: a fast mode associated with adjustment of the grounding zone and a slow mode associated with the interior reservoir. We have shown here that ocean forcing projects onto both modes, while interior forcing projects almost entirely onto the slow mode.

Persistence in stochastic climate forcing amplifies natural groundingline fluctuations (Fig. 4). Increased persistence means more of the variance in the forcing occurs at low frequencies, where the glacier is ultimately more sensitive. In particular, multidecadal ocean variability can drive large fluctuations on multidecadal to millennial timescales. Understanding the magnitude and persistence of ocean forcing is thus critical for attributing observed terminus changes and detecting the response to anthropogenic forcing (Fig. 6).

Despite contrasting responses on short timescales, slow (i.e., millennial) dynamics dominate the full system response for both types of forcing. This implies a large committed response (Fig. 5), as well as a memory of past climate fluctuations (Fig. 7), for both types of forcing. The slow response of interior ice and the attendant consequences have been explored in a number of previous studies (e.g., Nye, 1960; Levermann et al., 2013; Robel et al., 2018); the results presented here demonstrate how a slow response is fundamental to ocean forcing as well.
Given the rapid observed changes linked to ocean forcing in the past few decades, it is useful to discuss points (3) and (4) further. Increases in discharge from Greenland outlet glaciers have been linked to regional climate variability and warming of the subpolar North Atlantic (Straneo and Heimbach, 2013). Records for a smaller selection of glaciers extend up to a century (Andresen et al., 2012; Bjørk et al., 2012) and further indicate sensitivity to regional oceanic forcing on relatively short timescales. Evidence for the response to decadal variability is strong in Antarctica, too (e.g., Jenkins et al., 2016, 2018). In sum, the capacity for outletglacier grounding lines to react quickly to changes at the terminus is quite clear from the observational record.
Our results show that climate and glacier fluctuations over the historical record also have a longertimescale context. The widely separated glacier response timescales mean that fast fluctuations are essentially superimposed on millennial fluctuations. This is clear even in the flowline model, which does not make an explicit approximation of only two timescales (e.g., Fig. 1d). Slow fluctuations are difficult to resolve in observational records, but the practical point is that short periods of terminus “stabilization”, or fluctuations about a shortterm mean state, do not necessarily preclude a large disequilibrium between fluxes near the grounding line and fluxes from the interior. This disequilibrium should be reflected in icethickness and velocity profiles and has been noted for some systems whose retreat or discharge rates have leveled off, including Jakobshavn Isbræ (e.g., Joughin et al., 2012b; Khazendar et al., 2019) and Pine Island Glacier (e.g., Christianson et al., 2016). Careful modeling might integrate such observations with climate reconstructions to inform attribution studies or model initialization. Wherever the forcing history is uncertain, multiple realizations of past variability should be considered, and we again emphasize that the persistence of this variability is a critical parameter.
A number of simplifying assumptions throughout our study warrant some discussion. First of all, we have assumed a constant, prograde bed slope. In reality, variations in bed topography can have a strong effect on retreat rates and sensitivity, both at the terminus (Catania et al., 2018) and in the interior (Felikson et al., 2017). The effect of these variations on the fast and slow modes remains to be investigated. Felikson et al. (2017) showed that bedrock knickpoints may serve as barriers to rapid inland thinning following forcing at the terminus. Depending on where the knickpoint is, this could mean that the dimensions of the effective interior reservoir are different for ocean vs. interior forcing, potentially changing the relevant response timescales.
A related issue is the instability associated with retrograde slopes (e.g., Weertman, 1974; Schoof, 2007a). In an unstable configuration, the linearized timescales diverge, but RRH showed that fast and slow tendencies still govern the rates of unstable retreat. More recently, Robel et al. (2019) showed that instability magnifies variability in transient retreat rate due to internal climate variability. Together, these points suggest that instability might also magnify the difference between retreats initiated by interior and groundinglineflux anomalies.
Another simplification is that we have focused on spatially uniform interior forcings, whereas in reality, surfacemassbalance anomalies are likely to be greater near the marine margin. We conducted several experiments with the flowline model in which S anomalies were concentrated towards the grounding line and were also amplified to produce the same volume anomalies. This enhanced the highfrequency length response, although there was still a clear difference from the transient response to Ω forcing. Compared to a uniform S anomaly, a step change in S concentrated on the lower half of the glacier (with doubled magnitude) nearly doubles the groundingline response after 100 yr, while the 100 yr response following a step in Ω is 4 times greater. For whitenoise forcing, σ_{L} increases ∼12 % if S anomalies are concentrated in the same way but increases >50 % if applied as Ω anomalies. Thus, the temporal distinction between interior and ocean forcing will depend on the spatial pattern of mass balance, but it would likely take extremely large and concentrated anomalies to match the fast response of Ω forcing. We expect that the distinction would also depend on the horizontal catchment geometry: convergence would amplify the effects of anomalies from the deep interior even if they are smaller in magnitude. Further experiments with 2D geometries would help to clarify these potentially competing effects. We have also neglected orographic and surface–elevation feedbacks, which would depend on the spatial pattern of thinning and may thus affect fast and slow responses differently. These would provide another interesting avenue for future analyses.
Finally, we chose to impose ocean forcing via the groundinglineflux coefficient, which adjusts the parameterized relationship between ice thickness and groundingline flux. This allowed us to compare flux perturbations in a very general way, but it would be more realistic to directly force a dynamic ice shelf or to perturb calving processes. The analytical flux conditions of Haseloff and Sergienko (2018) could be an intermediate route and would be feasible for models that parameterize groundingline flux. RRH showed that these conditions can introduce nonlinear sensitivity to perturbations in a buttressing ice shelf and thus another key difference between glacier responses to interior and ocean forcing. For example, this nonlinearity would skew the distribution of fluctuations driven by natural ocean variability, changing the thresholds for trend detection (see Fig. 6).
Previous studies have employed a variety of forcing strategies for glaciers without buttressing ice shelves. Some have directly perturbed stresses at the grounding line (e.g., Nick et al., 2009; Price et al., 2011), while others have implemented forcing via frontal melt rates (e.g., Morlighem et al., 2016; Aschwanden et al., 2019). Another approach is to impose terminus retreat based on observations and then allow upglacier ice thickness and flux to evolve. This approach has been used to isolate the effects of bed topography on inland thinning (Felikson, 2018).
Regardless of how forcing is implemented, our analyses show that the coupling between groundingzone and interior dynamics is a key part of an outlet glacier's response to ocean forcing. In particular, increased discharge is what eventually precipitates the slow (but large) second stage of groundingline retreat. Our idealized ocean forcing has limitations, of course, but the mechanism for a second stage of retreat is physically robust: elevated interior fluxes must eventually fall as the interior drains, creating a tendency towards further groundingline retreat (Fig. 3). The actual retreat could be modulated depending on the bed topography through which the grounding line migrates, and this would be yet another factor to investigate with more realistic geometries.
At the icesheet scale, ocean forcing must often be simplified considerably, and the optimal strategy remains to be determined. The Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6; Nowicki et al., 2016) defines several experimental protocols for use across various icesheet models. The ISMIP6 oceanforcing parameterizations for Greenland were recently established, based on an empirical study of a large collection of terminus observations and regional climate data (Slater et al., 2019, 2020). They proposed two strategies: one in which length change is imposed directly as a function of climate forcing (“retreat implementation”) and one in which submarine melt rates are prescribed (“melt implementation”). Both leave the evolution of ice flux up to the icesheet model, but the retreat implementation does not allow interior ice dynamics to feed back into the terminus position. The results presented here suggest that, because the timescale of interior dynamics is much longer than the observational records used to calibrate the imposed retreats, the retreat implementation would project less longterm terminus retreat compared with the melt implementation.
The impact that the choice of parameterization would have on projections of nearterm sealevel rise is not immediately clear, however. Slater et al. (2020) note that these parameterizations are intended primarily for 21stcentury projections. In a scenario with severe warming, the fast response to progressive forcing might plausibly outweigh the contribution of the slow mode over this period. On the other hand, in a scenario of climate stabilization, the two approaches could yield more disparate results. Under the retreat implementation, the terminus stabilizes if the climate stabilizes (Slater et al., 2019), whereas under the melt implementation, it could in principle continue retreating due to disequilibrium between interior and groundingzone fluxes (much like in Fig. 5). In such a case, the choice of parameterization would affect the difference between responses to lowemission versus highemission scenarios (e.g., representative concentration pathway 2.6 vs. 8.5). The principles from our idealized experiments suggest that a fixed versus free terminus would yield a different partition between fast and slow glacier responses, but further experimentation is needed to investigate this. Although much attention is directed towards sealevel rise by 2100, longertimescale comparisons within the ISMIP6 framework, with both idealized and comprehensive models, may also prove illuminating. Understanding the timescales of glacier dynamics is a longstanding pursuit in glaciology (e.g., Nye, 1960; Jóhannesson et al., 1989; Oerlemans, 2001; Roe and Baker, 2014; Robel et al., 2018). Here, we have explored for marineterminating outlet glaciers how these timescales manifest under different types of forcing. Interpreting observations and making useful projections will always depend partly on characteristics that are unique to each outlet glacier and on processes that remain to be investigated further. However, the basic constraints of flux conservation and largescale geometry will always apply and, as we have shown here, are enough on their own to yield consequential differences in glacier response. Obviously, there are many intermediate levels of complexity between the models we have used here and those used for detailed simulations and projections. Alongside the need for more comprehensive glacier and icesheet projections in a warming climate, there also remains the need to improve understanding by analyzing model experiments positioned throughout this hierarchy of model complexity (e.g., Held, 2005). Reduced models and experiments with idealized glacier geometries, in which the fundamentals are laid bare, can thus serve as a useful foundation for the needed analyses.
Here we present the linearized twostage model in more detail, and derive a discrete autoregressive form. Again, a full model derivation can be found in RRH; we start here from the linearized equations with perturbations in Q_{g} and S:
A_{H}, A_{L}, B_{H}, and B_{L} are the linear couplings between length and thickness anomalies:
where $\mathit{\lambda}=\frac{{\mathit{\rho}}_{w}}{{\mathit{\rho}}_{i}}$. For numerical implementation, Eqs. (A1) and (A2) can be discretized using a backward Euler method, where Δt is the time step between glacier states [i] and [i−1]:
Solving for H_{[i]} and L_{[i]} gives
which are still semiimplicit in terms of the current state [i]. However, we can substitute Eq. (A10) for the ${L}_{\left[i\right]}^{\prime}$ term in Eq. (A9) and vice versa for H_{[i]}. Solving again for H_{[i]} and L_{[i]}, we arrive at explicit expressions that depend on only the past state and current forcing. This takes the form of a twodimensional autoregressive process (e.g., Box et al., 2008) with two forcing terms. For compactness, we redefine several parameter combinations at this point. Let $\mathit{\u03f5}=(\mathrm{1}{B}_{L}\mathrm{\Delta}t{)}^{\mathrm{1}}$, $\mathit{\eta}=(\mathrm{1}{A}_{H}\mathrm{\Delta}t{)}^{\mathrm{1}}$, and $\mathit{\kappa}=(\mathrm{1}\mathit{\eta}\mathit{\u03f5}{A}_{L}{B}_{H}\mathrm{\Delta}{t}^{\mathrm{2}}{)}^{\mathrm{1}}$. Then,
where the autoregressive coefficients are contained in the matrix C and vectors d and e:
This form lends itself to a straightforward solution algorithm, given steadystate glacier parameters and forcing time series as inputs. Note that for the ocean forcing used in this study, ${Q}_{\mathrm{g}}^{\prime}$ is proportional to Ω^{′}, but other relationships could easily be implemented.
Code used for the analyses in this study is available as a public repository at https://github.com/johnerich/outletglacforcing (last access: 12 July 2020; Christian, 2020). Output from the flowline model is available from the corresponding author upon request.
All authors contributed to the study design. JEC carried out the analyses and wrote the manuscript, with input from all authors.
The authors declare that they have no conflict of interest.
We are very grateful to David Pollard (Penn State University) for model code and a generous introduction to using the flowline model. We also thank Ginny Catania, Fiamma Straneo, and Donald Slater for insightful conversations on outletglacier dynamics. Finally, we thank Martin Lüthi and the anonymous reviewer for comments that helped clarify the manuscript, and we thank the editor Olivier Gagliardini. John Erich Christian was supported by the NSF Graduate Research Fellowship Program.
This research has been supported by the National Science Foundation, Division of Graduate Education (grant no. DGE1256082) and National Science Foundation, Office of Polar Programs (grant no. OPP1542756).
This paper was edited by Olivier Gagliardini and reviewed by Martin Lüthi and one anonymous referee.
Abram, N. J., McGregor, H. V., Tierney, J. E., Evans, M. N., McKay, N. P., Kaufman, D. S., Thirumalai, K., Martrat, B., Goosse, H., Phipps, S. J., Steig, E. J., Kilbourne, K. H., Saenger, C. P., Zinke, J., Leduc, G., Addison, J. A., Mortyn, P. G., Seidenkrantz, M. S., Sicre, M. A., Selvaraj, K., Filipsson, H. L., Neukom, R., Gergis, J., Curran, M. A., and Von Gunten, L.: Early onset of industrialera warming across the oceans and continents, Nature, 536, 411–418, https://doi.org/10.1038/nature19082, 2016. a
Andresen, C. S., Straneo, F., Ribergaard, M. H., Bjørk, A. A., Andersen, T. J., Kuijpers, A., NørgaardPedersen, N., Kjær, K. H., Schjøth, F., Weckström, K., and Ahlstrøm, A. P.: Rapid response of Helheim Glacier in Greenland to climate variability over the past century, Nat. Geosci., 5, 37–41, https://doi.org/10.1038/ngeo1349, 2012. a
Armour, K. C. and Roe, G. H.: Climate commitment in an uncertain world, Geophys. Res. Lett., 38, F02030, https://doi.org/10.1029/2010GL045850, 2011. a
Aschwanden, A., Fahnestock, M. A., Truffer, M., Brinkerhoff, D. J., Hock, R., Khroulev, C., Mottram, R., and Abbas Khan, S.: Contribution of the Greenland Ice Sheet to sea level over the next millennium, Sci. Adv., 5, eaav9396, https://doi.org/10.1126/sciadv.aav9396, 2019. a, b
Bindschadler, R. A., Nowicki, S., AbeOUCHI, A., Aschwanden, A., Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Khroulev, C., Levermann, A., Lipscomb, W. H., Martin, M. A., Morlighem, M., Parizek, B. R., Pollard, D., Price, S. F., Ren, D., Saito, F., Sato, T., Seddik, H., Seroussi, H., Takahashi, K., Walker, R., and Wang, W. L.: Icesheet model sensitivities to environmental forcing and their use in projecting future sea level (the SeaRISE project), J. Glaciol., 59, 195–224, https://doi.org/10.3189/2013JoG12J125, 2013. a
Bjørk, A. A., Kjær, K. H., Korsgaard, N. J., Khan, S. A., Kjeldsen, K. K., Andresen, C. S., Box, J. E., Larsen, N. K., and Funder, S.: An aerial view of 80 years of climaterelated glacier fluctuations in southeast Greenland, Nat. Geosci., 5, 427–432, https://doi.org/10.1038/ngeo1481, 2012. a
Box, G. E. P., Jenkins, G. M., and Reinsel, G. C.: Time Series Analysis: Forecasting and Control: Fourth Edition, John Wiley, https://doi.org/10.1002/9781118619193, 2008. a, b
Buizert, C., Keisling, B. A., Box, J. E., He, F., Carlson, A. E., Sinclair, G., and DeConto, R. M.: GreenlandWide Seasonal Temperatures During the Last Deglaciation, Geophys. Res. Lett., 45, 1905–1914, https://doi.org/10.1002/2017GL075601, 2018. a, b, c
Catania, G. A., Stearns, L. A., Sutherland, D. A., Fried, M. J., Bartholomaus, T. C., Morlighem, M., Shroyer, E., and Nash, J.: Geometric Controls on Tidewater Glacier Retreat in Central Western Greenland, J. Geophys. Res.Earth Surf., 123, 2024–2038, https://doi.org/10.1029/2017JF004499, 2018. a, b
Christian, J. E.: outletglacforcing, GitHub, available at: https://github.com/johnerich/outletglacforcing, last access: 12 July 2020. a
Christian, J. E., Koutnik, M., and Roe, G.: Committed retreat: Controls on glacier disequilibrium in a warming climate, J. Glaciol., 64, 675–688, https://doi.org/10.1017/jog.2018.57, 2018. a, b
Christianson, K., Bushuk, M., Dutrieux, P., Parizek, B. R., Joughin, I. R., Alley, R. B., Shean, D. E., Abrahamsen, E. P., Anandakrishnan, S., Heywood, K. J., Kim, T. W., Lee, S. H., Nicholls, K., Stanton, T., Truffer, M., Webber, B. G., Jenkins, A., Jacobs, S., Bindschadler, R., and Holland, D. M.: Sensitivity of Pine Island Glacier to observed ocean forcing, Geophys. Res. Lett., 43, 10817–10825, https://doi.org/10.1002/2016GL070500, 2016. a
Cook, A. J., Holland, P. R., Meredith, M. P., Murray, T., Luckman, A., and Vaughan, D. G.: Ocean forcing of glacier retreat in the western Antarctic Peninsula, Science, 353, 283–286, https://doi.org/10.1126/science.aae0017, 2016. a
Csatho, B. M., Schenka, A. F., Van Der Veen, C. J., Babonis, G., Duncan, K., Rezvanbehbahani, S., Van Den Broeke, M. R., Simonsen, S. B., Nagarajan, S., and Van Angelen, J. H.: Laser altimetry reveals complex pattern of Greenland Ice Sheet dynamics, P. Nat. Acad. Sci. USA, 111, 18478–18483, https://doi.org/10.1073/pnas.1411680112, 2014. a
Deser, C., Alexander, M. A., and Timlin, M. S.: Understanding the Persistence of Sea Surface Temperature Anomalies in Midlatitudes, J. Climate, 16, 57–72, https://doi.org/10.1175/15200442(2003)016<0057:UTPOSS>2.0.CO;2, 2003. a
Felikson, D.: Geometric controls on the inland extent of dynamic thinning for Greenland Ice Sheet outlet glaciers, Phd thesis, University of Texas, Austin, available at: http://hdl.handle.net/2152/72430 (last access: 12 July 2020), 2018. a
Felikson, D., Bartholomaus, T. C., Catania, G. A., Korsgaard, N. J., Kjær, K. H., Morlighem, M., Noël, B., Van Den Broeke, M., Stearns, L. A., Shroyer, E. L., Sutherland, D. A., and Nash, J. D.: Inland thinning on the Greenland ice sheet controlled by outlet glacier geometry, Nat. Geosci., 10, 366–369, https://doi.org/10.1038/ngeo2934, 2017. a, b, c
Glen, J. W.: The creep of polycrystalline ice, P. Roy. Soc. London A, 228, 519–538, https://doi.org/10.1098/rspa.1955.0066, 1955. a
Goelzer, H., Nowicki, S., Edwards, T., Beckley, M., AbeOuchi, A., Aschwanden, A., Calov, R., Gagliardini, O., GilletChaulet, F., Golledge, N. R., Gregory, J., Greve, R., Humbert, A., Huybrechts, P., Kennedy, J. H., Larour, E., Lipscomb, W. H., Le clec'h, S., Lee, V., Morlighem, M., Pattyn, F., Payne, A. J., Rodehacke, C., Rückamp, M., Saito, F., Schlegel, N., Seroussi, H., Shepherd, A., Sun, S., van de Wal, R., and Ziemen, F. A.: Design and results of the ice sheet model initialisation experiments initMIPGreenland: an ISMIP6 intercomparison, The Cryosphere, 12, 1433–1460, https://doi.org/10.5194/tc1214332018, 2018. a
Goldberg, D. N., Heimbach, P., Joughin, I., and Smith, B.: Committed retreat of Smith, Pope, and Kohler Glaciers over the next 30 years inferred by transient model calibration, The Cryosphere, 9, 2429–2446, https://doi.org/10.5194/tc924292015, 2015. a
Gwyther, D. E., O'Kane, T. J., GaltonFenzi, B. K., Monselesan, D. P., and Greenbaum, J. S.: Intrinsic processes drive variability in basal melting of the Totten Glacier Ice Shelf, Nat. Commun., 9, 3141, https://doi.org/10.1038/s41467018056182, 2018. a
Hansen, J., Russell, G., Lacis, A., Fung, I., Rind, D., and Stone, P.: Climate response times: Dependence on climate sensitivity and ocean mixing, Science, 229, 857–859, https://doi.org/10.1126/science.229.4716.857, 1985. a
Haseloff, M. and Sergienko, O. V.: The effect of buttressing on grounding line dynamics, J. Glaciol., 64, 417–431, https://doi.org/10.1017/jog.2018.30, 2018. a, b, c
Hasselmann, K.: Stochastic climate models Part I. Theory, Tellus, 28, 473–485, https://doi.org/10.3402/tellusa.v28i6.11316, 1976. a
Held, I. M.: The gap between simulation and understanding in climate modeling, B. Am. Meteorol. Soc., 86, 1609–1614, https://doi.org/10.1175/BAMS86111609, 2005. a
Howat, I. M., Joughin, I., Fahnestock, M., Smith, B. E., and Scambos, T. A.: Synchronous retreat and acceleration of southeast Greenland outlet glaciers 200006: Ice dynamics and coupling to climate, J. Glaciol., 54, 646–660, https://doi.org/10.3189/002214308786570908, 2008. a
Huybers, P. and Curry, W.: Links between annual, Milankovitch and continuum temperature variability, Nature, 441, 329–332, https://doi.org/10.1038/nature04745, 2006. a, b
IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Tech. rep., Cambridge, https://doi.org/10.1017/CBO9781107415324, 2013. a
Jenkins, A., Dutrieux, P., Jacobs, S., Steig, E. J., Gudmundsson, G. H., Smith, J., and Heywood, K. J.: Decadal ocean forcing and Antarctic Ice Sheet response: lessons from the Amundsen Sea, Oceanography, 29, 106–117, 2016. a, b
Jenkins, A., Shoosmith, D., Dutrieux, P., Jacobs, S., Kim, T. W., Lee, S. H., Ha, H. K., and Stammerjohn, S.: West Antarctic Ice Sheet retreat in the Amundsen Sea driven by decadal oceanic variability, Nat. Geosci., 11, 733–738, https://doi.org/10.1038/s4156101802074, 2018. a, b
Jóhannesson, T., Raymond, C., and Waddington, E.: TimeScale for Adjustment of Glaciers to Changes in Mass Balance, J. Glaciol., 35, 355–369, https://doi.org/10.3189/s002214300000928x, 1989. a, b
Joughin, I., Alley, R. B., and Holland, D. M.: Icesheet response to oceanic forcing, Science, 338, 1172–1176, https://doi.org/10.1126/science.1226481, 2012a. a
Joughin, I., Smith, B. E., Howat, I. M., Floricioiu, D., Alley, R. B., Truffer, M., and Fahnestock, M.: Seasonal to decadal scale variations in the surface velocity of Jakobshavn Isbrae, Greenland: Observation and modelbased analysis, J. Geophys. Res.Earth Surf., 117, F02030, https://doi.org/10.1029/2011JF002110, 2012b. a
Khazendar, A., Fenty, I. G., Carroll, D., Gardner, A., Lee, C. M., Fukumori, I., Wang, O., Zhang, H., Seroussi, H., Moller, D., Noël, B. P., van den Broeke, M. R., Dinardo, S., and Willis, J.: Interruption of two decades of Jakobshavn Isbrae acceleration and thinning as regional ocean cools, Nat. Geosci., 12, 277–283, https://doi.org/10.1038/s4156101903293, 2019. a
Levermann, A., Clark, P. U., Marzeion, B., Milne, G. A., Pollard, D., Radic, V., and Robinson, A.: The multimillennial sealevel commitment of global warming, P. Natl. Acad. Sci. USA, 110, 13745–13750, https://doi.org/10.1073/pnas.1219414110, 2013. a, b
Lingle, C. S.: A numerical model of interactions between a polar ice stream and the ocean: Application to ice stream E, West Antarctica, J. Geophys. Res., 89, 3523–3549, https://doi.org/10.1029/JC089iC03p03523, 1984. a
Mantelli, E., Bertagni, M. B., and Ridolfi, L.: Stochastic ice stream dynamics, P. Natl. Acad. Sci. USA, 113, E4594–E4600, https://doi.org/10.1073/pnas.1600362113, 2016. a
Medwedeff, W. G. and Roe, G. H.: Trends and variability in the global dataset of glacier mass balance, Clim. Dynam., 48, 3085–3097, https://doi.org/10.1007/s003820163253x, 2017. a
Meehl, G. A., Washington, W. M., Collins, W. D., Arblaster, J. M., Hu, A., Buja, L. E., Strand, W. G., and Teng, H.: How much more global warming and sea level rise?, Science, 307, 1769–1772, https://doi.org/10.1126/science.1106663, 2005. a
Miles, B. W., Stokes, C. R., Vieli, A., and Cox, N. J.: Rapid, climatedriven changes in outlet glaciers on the Pacific coast of East Antarctica, Nature, 500, 563–566, https://doi.org/10.1038/nature12382, 2013. a
Moon, T. and Joughin, I.: Changes in ice front position on Greenland's outlet glaciers from 1992 to 2007, J. Geophys. Res., 113, F02022, https://doi.org/10.1029/2007JF000927, 2008. a
Moon, T., Joughin, I., Smith, B., and Howat, I.: 21stcentury evolution of Greenland outlet glacier velocities, Science, 336, 576–578, https://doi.org/10.1126/science.1219985, 2012. a
Moon, T., Joughin, I., Smith, B., Van Den Broeke, M. R., Van De Berg, W. J., Noël, B., and Usher, M.: Distinct patterns of seasonal Greenland glacier velocity, Geophys. Res. Lett., 41, 7209–7216, https://doi.org/10.1002/2014GL061836, 2014. a, b
Morlighem, M., Bondzio, J., Seroussi, H., Rignot, E., Larour, E., Humbert, A., and Rebuffi, S.: Modeling of Store Gletscher's calving dynamics, West Greenland, in response to ocean thermal forcing, Geophys. Res. Lett., 43, 2659–2666, https://doi.org/10.1002/2016GL067695, 2016. a
Nick, F. M., Vieli, A., Howat, I. M., and Joughin, I.: Largescale changes in Greenland outlet glacier dynamics triggered at the terminus, Nat. Geosci., 2, 110–114, https://doi.org/10.1038/ngeo394, 2009. a, b
Nielsen, L. T., Adalgeirsdóttir, G., Gkinis, V., Nuterman, R., and Hvidberg, C. S.: The effect of a Holocene climatic optimum on the evolution of the Greenland ice sheet during the last 10 kyr, J. Glaciol., 64, 477–488, https://doi.org/10.1017/jog.2018.40, 2018. a
Nowicki, S. M. J., Payne, A., Larour, E., Seroussi, H., Goelzer, H., Lipscomb, W., Gregory, J., AbeOuchi, A., and Shepherd, A.: Ice Sheet Model Intercomparison Project (ISMIP6) contribution to CMIP6, Geosci. Model Dev., 9, 4521–4545, https://doi.org/10.5194/gmd945212016, 2016. a
Nye, J. F.: The response of glaciers and icesheets to seasonal and climatic changes, P. Roy. Soc. London A, 256, 559–584, https://doi.org/10.1098/rspa.1960.0127, 1960. a, b
Oerlemans, J. J.: Glaciers and climate change, A.A. Balkema Publishers, ISBN 10: 9026518137, 2001. a, b
Percival, D. B., Overland, J. E., and Mofjeld, H. O.: Interpretation of North Pacific Variability as a Shortand LongMemory Process, J. Climate, 14, 4545–4559, https://doi.org/10.1175/15200442(2001)014<4545:IONPVA>2.0.CO;2, 2001. a, b
Pollard, D. and DeConto, R. M.: Description of a hybrid ice sheetshelf model, and application to Antarctica, Geosci. Model Dev., 5, 1273–1295, https://doi.org/10.5194/gmd512732012, 2012. a, b
Price, S. F., Payne, A. J., Howat, I. M., and Smith, B. E.: Committed sealevel rise for the next century from Greenland ice sheet dynamics during the past decade, P. Natl. Acad. Sci. USA, 108, 8978–8983, https://doi.org/10.1073/pnas.1017313108, 2011. a, b
Pritchard, H. D., Arthern, R. J., Vaughan, D. G., and Edwards, L. A.: Extensive dynamic thinning on the margins of the Greenland and Antarctic ice sheets, Nature, 461, 971–975, https://doi.org/10.1038/nature08471, 2009. a
Robel, A. A., Schoof, C., and Tziperman, E.: Rapid grounding line migration induced by internal ice stream variability, J. Geophys. Res.Earth Surf., 119, 2430–2447, https://doi.org/10.1002/2014JF003251, 2014. a
Robel, A. A., Roe, G. H., and Haseloff, M.: Response of MarineTerminating Glaciers to Forcing: Time Scales, Sensitivities, Instabilities, and Stochastic Dynamics, J. Geophys. Res.Earth Surf., 123, 2205–2227, https://doi.org/10.1029/2018JF004709, 2018. a, b, c, d, e, f, g, h
Robel, A. A., Seroussi, H., and Roe, G. H.: Marine ice sheet instability amplifies and skews uncertainty in projections of future sealevel rise, P. Natl. Acad. Sci. USA, 116, 14887–14892, https://doi.org/10.1073/pnas.1904822116, 2019. a, b, c
Roe, G. H.: What do glaciers tell us about climate variability and climate change?, J. Glaciol., 57, 567–578, https://doi.org/10.3189/002214311796905640, 2011. a
Roe, G. H. and Baker, M. B.: Glacier response to climate perturbations: An accurate linear geometric model, J. Glaciol., 60, 670–684, https://doi.org/10.3189/2014JoG14J016, 2014. a
Roe, G. H. and Baker, M. B.: The response of glaciers to climatic persistence, J. Glaciol., 62, 440–450, https://doi.org/10.1017/jog.2016.4, 2016. a, b, c
Roe, G. H., Baker, M. B., and Herla, F.: Centennial glacier retreat as categorical evidence of regional climate change, Nat. Geosci., 10, 95–99, https://doi.org/10.1038/ngeo2863, 2017. a, b, c
Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res.Earth Surf., 112, https://doi.org/10.1029/2006JF000664, 2007a. a, b, c, d, e, f, g
Schoof, C.: Marine icesheet dynamics. Part 1. The case of rapid sliding, J. Fluid Mech., 573, 27–55, https://doi.org/10.1017/S0022112006003570, 2007b. a
Shepherd, A., Ivins, E., Rignot, E., Smith, B., Van Den Broeke, M., Velicogna, I., Whitehouse, P., Briggs, K., Joughin, I., Krinner, G., Nowicki, S., Payne, T., Scambos, T., Schlegel, N., Geruo, A., Agosta, C., Ahlstrøm, A., Babonis, G., Barletta, V., Blazquez, A., Bonin, J., Csatho, B., Cullather, R., Felikson, D., Fettweis, X., Forsberg, R., Gallee, H., Gardner, A., Gilbert, L., Groh, A., Gunter, B., Hanna, E., Harig, C., Helm, V., Horvath, A., Horwath, M., Khan, S., Kjeldsen, K. K., Konrad, H., Langen, P., Lecavalier, B., Loomis, B., Luthcke, S., McMillan, M., Melini, D., Mernild, S., Mohajerani, Y., Moore, P., Mouginot, J., Moyano, G., Muir, A., Nagler, T., Nield, G., Nilsson, J., Noel, B., Otosaka, I., Pattle, M. E., Peltier, W. R., Pie, N., Rietbroek, R., Rott, H., SandbergSørensen, L., Sasgen, I., Save, H., Scheuchl, B., Schrama, E., Schröder, L., Seo, K. W., Simonsen, S., Slater, T., Spada, G., Sutterley, T., Talpe, M., Tarasov, L., Van De Berg, W. J., Van Der Wal, W., Van Wessem, M., Vishwakarma, B. D., Wiese, D., and Wouters, B.: Mass balance of the Antarctic Ice Sheet from 1992 to 2017, Nature, 558, 219–222, https://doi.org/10.1038/s415860180179y, 2018. a
Slater, D. A., Straneo, F., Felikson, D., Little, C. M., Goelzer, H., Fettweis, X., and Holte, J.: Estimating Greenland tidewater glacier retreat driven by submarine melting, The Cryosphere, 13, 2489–2509, https://doi.org/10.5194/tc1324892019, 2019. a, b
Slater, D. A., Felikson, D., Straneo, F., Goelzer, H., Little, C. M., Morlighem, M., Fettweis, X., and Nowicki, S.: Twentyfirst century ocean forcing of the Greenland ice sheet for modelling of sea level contribution, The Cryosphere, 14, 985–1008, https://doi.org/10.5194/tc149852020, 2020. a, b
Straneo, F. and Heimbach, P.: North Atlantic warming and the retreat of Greenland's outlet glaciers, 504, 36–43, https://doi.org/10.1038/nature12854, 2013. a, b
Straneo, F., Sutherland, D. A., Holland, D., Gladish, C., Hamilton, G. S., Johnson, H. L., Rignot, E., Xu, Y., and Koppes, M.: Characteristics of ocean waters reaching greenland's glaciers, Ann. Glaciol., 53, 202–210, https://doi.org/10.3189/2012AoG60A059, 2012. a
Thomas, R. H.: The Dynamics of Marine Ice Sheets, J. Glaciol., 24, 167–177, https://doi.org/10.3189/s0022143000014726, 1979. a
Tsai, V. C., Stewart, A. L., and Thompson, A. F.: Marine icesheet profiles and stability under Coulomb basal conditions, J. Glaciol., 61, 205–215, https://doi.org/10.3189/2015JoG14J221, 2015. a
Van Den Broeke, M., Bamber, J., Ettema, J., Rignot, E., Schrama, E., Van Berg, W. J. D., Van Meijgaard, E., Velicogna, I., and Wouters, B.: Partitioning recent Greenland mass loss, Science, 326, 984–986, https://doi.org/10.1126/science.1178176, 2009. a
Weertman, J.: Stability of the Junction of an Ice Sheet and an Ice Shelf, J. Glaciol., 13, 3–11, https://doi.org/10.3189/s0022143000023327, 1974. a, b, c
Wigley, T. M. L. and Schlesinger, M. E.: Wigley_Schlesinger_1985, Nature, 315, 649–652, https://doi.org/10.1038/315649a0, 1985. a