The contrasting response of outlet glaciers to interior and ocean forcing

The dynamics of marine-terminating 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 surface-mass-balance forcing applied over the interior than to oceanic forcing applied at the grounding line. A recently developed reduced model represents outlet-glacier dynamics via two widely separated response timescales: a fast response associated with grounding-zone 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 surface-mass-balance 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.

1 Introduction

Marine-terminating 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 sea-level 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 marine-terminating 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 ice-sheet 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, catchment-specific factors continue to pose a challenge where observations are limited, as well as for ice-sheet-wide 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 ice-flow 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 marine-ice-sheet 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 long-term 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 Part 1 – dynamical responses to ocean and interior forcing

2.1 A simple outlet-glacier system

Before describing the dynamic models, it is useful to begin with the geometry and basic flux-balance 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 forward-sloping (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

$$\begin{array}{}\text{(1)}& Q_{\mathrm{g}}=\mathrm{\Omega }{h}_{\mathrm{g}}^{\mathit{\beta }},\end{array}$$

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

$$\begin{array}{}\text{(2)}& h_{\mathrm{g}}\left(L\right)=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}}b\left(L\right),\end{array}$$

where ρ_w and ρ_i are the densities of seawater and ice, respectively.

In this study, we consider the grounding line the boundary of the outlet-glacier 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 grounding-line 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 outlet-glacier system, we begin with a 1-D (flowline) version of the ice-sheet 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 ice-flux divergence. Conservation of mass requires that

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial h}{\partial t}}=S-{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{u}h\right)}{\partial x}},\end{array}$$

where S is the local surface mass balance and $\stackrel{\mathrm{‾}}{u}$ is depth-averaged 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,

$$\begin{array}{}\text{(4)}& {\mathit{\tau }}_{\mathrm{d}}={\mathit{\rho }}_{\mathrm{i}}gh{\displaystyle \frac{\partial s}{\partial x}},\end{array}$$

where ρ_i is ice density and g is acceleration due to gravity. Longitudinal stretching $\frac{\partial u}{\partial x}$ is represented by the shallow-shelf approximation, where stretching and basal drag τ_b balance driving stress:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{2}h{A}^{-\frac{\mathrm{1}}{n}}|{\displaystyle \frac{\partial u}{\partial x}}{|}^{\frac{\mathrm{1}}{n-\mathrm{1}}}{\displaystyle \frac{\partial u}{\partial x}}\right)-{\mathit{\tau }}_{\mathrm{b}}={\mathit{\tau }}_{\mathrm{d}}.\end{array}$$

We assume a typical power-law rheology, with coefficient A and exponent n (e.g., Glen, 1955). Internal shear deformation ($\frac{\partial u}{\partial z}$) is represented by the shallow-ice approximation. At a depth h−z, this is

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial u}{\partial z}}=\mathrm{2}A{\mathit{\tau }}_{\mathrm{e}}^{n-\mathrm{1}}{\mathit{\rho }}_{\mathrm{i}}g\left(h\text{-}z\right){\displaystyle \frac{\partial s}{\partial x}},\end{array}$$

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 Weertman-type power-law relationship:

$$\begin{array}{}\text{(7)}& u_{\mathrm{b}}={\left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}}{C}}\right)}^{\frac{\mathrm{1}}{m}},\end{array}$$

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 shallow-ice and shallow-shelf 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):

$$\begin{array}{}\text{(8)}& \mathrm{\Omega }={\left[A\left({\mathit{\rho }}_{\mathrm{i}}g{)}^{n+\mathrm{1}}{\left(\mathrm{\Theta }\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}}\right)\right)}^n\right({\mathrm{4}}^{n}C{)}^{-\mathrm{1}}\right]}^{\frac{\mathrm{1}}{m+\mathrm{1}}},\end{array}$$

$$\begin{array}{}\text{(9)}& \mathit{\beta }={\displaystyle \frac{m+n+\mathrm{3}}{m+\mathrm{1}}}.\end{array}$$

Θ is a buttressing factor between 0 and 1, and all other parameters are defined as above. The PD12 model calculates grounding-line flux based on a thickness h_g that is linearly interpolated from the height above flotation of the last-grounded and first-floating grid cells to the point where flotation is reached. The corresponding subgrid grounding-line 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 continent-scale simulations over many millennia, we use a grid of 100 m to better resolve the details of grounding-line 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 steady-state profile for an idealized outlet glacier simulated with the flowline model. The domain begins at an ice divide and thus has a zero-flux 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⁻¹ 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 outlet-glacier catchments in Greenland, though direct comparisons are limited as the 1-D 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 surface-mass-balance 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 surface-mass-balance anomalies are assumed to be spatially uniform. We represent ocean forcing very simply by perturbing the grounding-line-flux 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 white-noise time series as either Ω or S anomalies but with the opposite sign so that the corresponding ice-volume anomalies match. The anomaly time series (hereafter Ω^′ and S^′) are scaled to have standard deviations equal to 20 % of the mean values ($\stackrel{\mathrm{‾}}{\mathrm{\Omega }}$ and $\stackrel{\mathrm{‾}}{S}$). Figure 1d shows the resulting length responses. For both types of forcing, the glacier acts as a low-pass filter on the imposed climate anomalies, producing kilometer-scale fluctuations with clear persistence. However, the length anomalies driven by Ω^′ have much greater high-frequency content, and greater overall variance, than those driven by S^′. The high-frequency response is perhaps intuitive, in that forcing applied at the grounding line has an immediate effect on grounding-line position. However, the response to Ω^′ also contains millennial-scale 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 marine-terminating 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 surface-mass-balance 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.

Table 1 Parameters used for flowline and two-stage models. Values for A and C follow previous idealized case studies (Schoof, 2007a; Robel et al., 2018).

Figure 1 Model schematic and responses to ocean forcing (via Ω) and interior forcing (via S). (a) An idealized marine-terminating glacier system. (b) Initial steady-state profile simulated by the flowline model. (c) Length responses to 20 % step increase in Ω (blue) and 20 % step decrease in S (orange) at t=0. Inset panel zooms into the first 0.5 kyr of response. Note the difference in initial retreat rates depending on the type of forcing. (d) Length responses to white-noise interannual variability in Ω or S. Right panel shows a 0.5 kyr segment in more detail. Anomalies (black) have a standard deviation of 20 % of the mean value and have opposite signs for Ω and S so that length anomalies are correlated for comparison. The response to ocean forcing (blue) contains much more high-frequency content than the response to interior forcing (orange).

2.4 The two-stage model

Robel, Roe, and Haseloff (2018; hereafter RRH) developed a simple model for marine-terminating outlet-glacier dynamics, derived from ice-flux conservation and constrained by large-scale 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 grounding-zone 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 state-dependent 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

$$\begin{array}{}\text{(10)}& Q={\left({\displaystyle \frac{{\mathit{\rho }}_{i}g}{C}}\right)}^n{\displaystyle \frac{H^{\mathit{\alpha }}}{L^{\mathit{\gamma }}}},\end{array}$$

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:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}=S-{\displaystyle \frac{Q}{L}}-{\displaystyle \frac{H}{h_{\mathrm{g}}L}}\left(Q-Q_{\mathrm{g}}\right),\end{array}$$

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{d}L}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{h_{\mathrm{g}}}}\left(Q-Q_{\mathrm{g}}\right).\end{array}$$

Because achieving a steady state requires adjustment of both H and L, we refer to this model as the “two-stage 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 two-stage model makes significant simplifications but stems from the same physical principles as the flowline model and most other contemporary ice-sheet models (see RRH for further discussion).

RRH also linearized these equations for fluctuations H^′ and L^′ about a steady state, $\stackrel{\mathrm{‾}}{H}$ and $\stackrel{\mathrm{‾}}{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

$$\begin{array}{}\text{(13)}& {\mathit{\tau }}_{\mathrm{F}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}{\stackrel{\mathrm{‾}}{S}}}{\left(\mathit{\alpha }+\mathit{\gamma }+\mathrm{1}-s_{\mathrm{T}}\right)}^{-\mathrm{1}},\end{array}$$

$$\begin{array}{}\text{(14)}& {\mathit{\tau }}_{\mathrm{S}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{H}}{\stackrel{\mathrm{‾}}{S}{s}_{\mathrm{T}}\mathit{\alpha }}}\left(\mathit{\alpha }+\mathit{\gamma }+\mathrm{1}-s_{\mathrm{T}}\right),\end{array}$$

where $s_{\mathrm{T}}=\mathrm{1}+\frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}\frac{\mathit{\beta }\stackrel{\mathrm{‾}}{b_x}\stackrel{\mathrm{‾}}{L}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}$ is a parameter describing the stability of glacier response based on its geometry and grounding-line dynamics (via β). Both response times contain a characteristic thickness divided by mass balance rate $\stackrel{\mathrm{‾}}{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 grounding-zone reservoir, whereas τ_S relates to the volume of the interior reservoir via $\stackrel{\mathrm{‾}}{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 two-stage model approximates the outlet-glacier system as linked interior- and grounding-zone reservoirs with distinct timescales.

Here we generalize the RRH linearization to simultaneously include perturbations to interior surface mass balance (S^′) and grounding-line flux ($Q_{\mathrm{g}}^{\prime }$; proportional to Ω^′). The linearized equations take the form of a two-dimensional dynamical system with two forcing terms:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left[\begin{array}{c}{H}^{\prime }\\ L^{\prime }\end{array}\right]& =\left[\begin{array}{cc}{A}_H& A_L\\ B_H& B_L\end{array}\right]\left[\begin{array}{c}{H}^{\prime }\\ L^{\prime }\end{array}\right]+\left[\begin{array}{c}{\stackrel{\mathrm{‾}}{L}}^{-\mathrm{1}}\left(\frac{\stackrel{\mathrm{‾}}{H}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}-\mathrm{1}\right)\\ {\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}^{-\mathrm{1}}\end{array}\right]{Q_{\mathrm{g}}}^{\prime }\\ & +\left[\begin{array}{c}\mathrm{1}\\ \mathrm{0}\end{array}\right]S^{\prime }.\end{array}\end{array}$$

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 two-dimensional autoregressive process (e.g., Box et al., 2008), which we also present in the Appendix.

2.5 Model comparison

RRH showed that the two-stage model emulated the response of a flowline model forced with surface-mass-balance anomalies. Their flowline model (described in Schoof, 2007b; Robel et al., 2014) used a stress-based, as opposed to flux-based, grounding-line condition but was otherwise dynamically similar to the PD12 model. Here, we use the PD12 model and extend the comparison to grounding-line-flux perturbations. Although the two-stage 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 two-stage 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 grounding-line retreat following a 20 % increase in Ω (blues) and 20 % decrease in S (orange and brown). For clarity, only the linearized two-stage output is shown. Note that the nonlinear two-stage 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 two-stage 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 white-noise 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{‾}}{\mathrm{\Omega }}$ or $\stackrel{\mathrm{‾}}{S}$. The two-stage and flowline models agree closely for both types of forcing, although the two-stage model slightly underestimates the magnitude of high-frequency 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 white-noise forcing has a flat power spectrum, the shapes of the response spectra reveal the glacier's filtering properties. The glacier is a strong low-pass filter for both forcing types but damps high frequencies even more for S^′ forcing. The two-stage model underestimates the high-frequency 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⁻² to 10⁻³ yr⁻¹, where the models agree quite well.

Figure 2 Comparing the transient response of flowline and two-stage models to ocean forcing (via Ω) and interior forcing (via S). (a) Schematic of the two-stage model. (b) Flowline and linearized two-stage model length responses to step changes in Ω (blues) and S (orange and brown). The thinner lines show the two-stage output. (c) Length anomalies in response to white-noise forcing in Ω (blues) and S (orange and brown). As in Fig. 1, anomalies have a standard deviation of 20 % of the mean value and have opposite signs. Of the 100 kyr model runs, 5 kyr is shown. (d) Autocorrelation function for the stochastic length fluctuations shown in (c). (e) Power spectral density for the stochastic length fluctuations, which is estimated throughout this study via Welch's method with windows of 1∕16 of the time series length and 50 % overlap. Note the clear split between fluctuations driven by Ω^′ (blues) and S^′ (orange and brown) at millennial and shorter timescales. This split is robust across models. The flowline model has more high-frequency power for both types of forcing, but note that the spectral power is orders-of-magnitude lower at such high frequencies.

2.6 Interpretations

The two-stage model's agreement with the flowline model suggests that it captures the essential dynamics responsible for the contrasting transient responses to interior and grounding-line-flux anomalies. At this point, it is useful to discuss some of the interpretations enabled by this reduced model.

The linearized two-stage 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 two-stage low-pass 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 high-frequency 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 grounding-line 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 two-stage model was derived from mass conservation, and that the linear response times (eigenmodes) reflect the large-scale glacier geometry. Because a glacier is a Stokes (i.e., noninertial) fluid driven by potential-energy gradients, the large-scale 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 up-glacier of the grounding line (Fig. 3a). For the nonlinear two-stage 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 surface-mass-balance and interior fluxes, while Q_g can quickly keep pace with flux anomalies from the interior. In the two-stage 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 two-stage model) brings interior fluxes back down, and again, Q_g must follow via grounding-line retreat.

Although the flowline model captures a more realistic and spatially distributed response, the two-stage model contains enough geometric information to emulate the basic sequence of flux anomalies. Note that more discretized interior reservoirs could be added to the two-stage 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 two-mode framework is useful for understanding the response to forcing applied at the grounding zone.

Figure 3 Flux changes following interior vs. ocean forcing. (a) Flowline model flux anomalies to perturbations in S and Ω. The schematic (top) shows the flux gates considered. Below, normalized flux responses to the perturbations are shown for each flux gate. (b) As for (a) but for the two-stage model, which has three flux variables. The two-stage model's interior flux variable (Q) captures the basic tendencies of interior fluxes in the flowline model. In both cases, the response to a perturbation in Ω clearly shows the two-timescale nature of outlet-glacier dynamics.

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, surface-mass-balance 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 surface-mass-balance 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 decade-scale ocean variability.

We consider four types of synthetic forcing with different persistence characteristics: (1) white noise; (2) and (3) first-order autoregressive (AR-1) processes with persistence timescales τ_AR1 of 4 and 20 yr, respectively; and (4) power-law persistence. These are all plausible models for the sort of persistence that may affect outlet-glacier forcing. For example, an AR-1 process with a 4 yr memory is characteristic of sea-surface-temperature 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 surface-temperature variability from interannual to multimillennial timescales has been described with a power-law 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 AR-1 process as a function of frequency, f, is

$$\begin{array}{}\text{(16)}& P_r=P_{\mathrm{0}}{\left(\mathrm{1}+r^{\mathrm{2}}-\mathrm{2}r\cos \left(\mathrm{2}\mathit{\pi }\mathrm{\Delta }tf\right)\right)}^{-\mathrm{1}},\end{array}$$

where P₀ scales the total variance and r is the auto-correlation 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}$. Power-law noise has a spectrum defined by

$$\begin{array}{}\text{(17)}& P_{\mathit{\nu }}=P_{\mathrm{0}}{\left({\displaystyle \frac{f_{\mathrm{0}}}{f}}\right)}^{\mathit{\nu }},\end{array}$$

where f₀ 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⁵ yr synthetic forcing time series, and Fig. 4c shows their power spectra. The AR-1 processes (red, maroon) have clear persistence in the time domain. The low-frequency content of power-law noise (blue) is hard to discern at short timescales, but the spectra show that it has similar power to the AR-1 processes at millennial timescales (and, again, identical total variance). Power-law variability is thus difficult to constrain from century-scale 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³ yr samples of the resulting length responses for the two-stage model. The effects of persistence are dramatic: for both S^′ or Ω^′ forcing, 4 yr AR-1 persistence more than doubles σ_L compared to length fluctuations driven by white-noise forcing, and both 20 yr AR-1 and power-law persistence cause a ∼5-to-6-fold 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 white-noise 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 terminus-flux anomalies excite the fast mode of response more than surface-mass-balance anomalies, the implication is that ocean variability – which tends to have persistence – may drive much larger terminus fluctuations than surface-mass-balance variability.

Figures 1–3 demonstrate that marine-terminating 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.

Figure 4 The effect of climatic persistence on glacier fluctuations. (a) Segments of synthetic anomalies with different levels of persistence: white noise (no persistence; black), red noise with τ_AR1=4 yr (red) and τ_AR1=20 yr (maroon), and power-law noise with ν=0.5 (blue). Time series shown are dimensionless and offset for clarity, with vertical ticks of 1σ. (b) The glacier-length variations generated by the two-stage model, when the synthetic anomalies are applied as Ω^′ (dark lines) and S^′ (pastel lines). Time series are again offset, with vertical ticks of 1 km. (c) Power spectra of the synthetic climate forcings. For reference, vertical lines correspond to $\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{F}}}$ and $\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{S}}}$. (d) Spectra of the glacier's response to white noise, which illustrates the glacier's sensitivity as a function of frequency and thus its temporal filtering properties. This depends on whether forcing is applied as Ω^′ (black) or S^′ (gray). (e) Spectra of glacier response to Ω^′ (dark lines) or S^′ (pastel lines). Responses depend on the spectra of the forcing and the glacier's filtering properties. Forcing with persistence has more power at low frequencies where glacier sensitivity is highest, increasing the overall glacier variance.

3 Part 2 – implications

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₂ 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). Industrial-era 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 two-stage 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 grounding-line 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 grounding-line 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 long-term commitment associated with the slow drawdown of interior ice.

The slow response (interior thinning → reduced interior fluxes → grounding-line 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 near-term sea-level rise and did not project past 2100. Our flux-balance 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 long-term retreat moves the terminus to more- or less-stable 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” short-term 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).

Table 2 Parameters varied between three idealized glaciers (top) and the resulting steady-state values (bottom). The steady state is calculated by the nonlinear two-stage model, and the linearized response times are given by Eqs. (13) and (14).

Figure 5 Glacier responses to idealized climate forcing over the industrial era. (a) Responses for glacier 1 (see Table 2). A 30 % change in S (brown) or Ω (blue) is realized as a linear trend from 1880 to 2020 CE. Dotted lines show the grounding-line position that would balance accumulation and grounding-line fluxes. (b) As for (a) but for glacier 2, which has faster response times and smaller equilibrium sensitivity. (c) As for (a) and (b), for glacier 3, which has slower response times and greater equilibrium sensitivity. (d) The range of responses of glacier 1 due to an uncertain forcing, idealized as a trend in Ω reaching a 20 %–40 % change by 2020. The slow transient response means a small range in absolute length change in the modern era, but uncertainty in the forcing projects directly onto the committed equilibrium response. Note the break in the time-axis scale at 2100.

3.2 The emergence and detectability of forced responses

We now turn to the topic of attributing outlet-glacier 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 signal-to-noise 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 high-frequency 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 AR-1 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 grounding-line anomalies of the test glacier described in Sect. 2, generated with the two-stage 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 grounding-line fluctuations, determined from 10⁷ 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 two-stage 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⁴ years, each of which are forced with stationary noise and no trend. For each realization, ΔL_null is the trend in grounding-line position over the observational period. We find ΔL_null for 10⁴ 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 shorter-term fluctuations that may obscure the early response to anthropogenic forcing, especially if the variability has significant persistence. The sensitivity of grounding-line 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.

Figure 6 Detecting the response to a climate trend in the presence of natural variability. Three types of natural variability are considered in each panel. The top row corresponds to interannual variability in S, the middle row to interannual variability in Ω, and the bottom row to multidecadal variability (τ_AR1=20 yr) in Ω. (a) The idealized climate trends (plus variability) beginning in 1880. In all cases, the linear trend reaches a SNR of 1 by 2020. (b) Grounding-line responses to each idealized trend. Shaded regions are the 1σ_L and 2σ_L bounds for each type of noise. The orange lines indicate when the trend was applied. Thinner lines show the grounding-line response without variability. (c) Probability density functions for grounding-line trends driven by each type of natural variability but no external forcing over time periods from 50 to 500 years. Note the different length scales in each case. Trends on the order of kilometers per century are extremely unlikely to occur due to variability in S alone (top) but commonplace if the glacier is sensitive to multidecadal variability in Ω (bottom).

3.3 Inherited conditions

We have thus far considered outlet-glacier fluctuations due to stationary interannual to multidecadal climate variability. However, long response times also imply some memory of climate variations throughout the Holocene. Ice-sheet 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 ice-sheet 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 outlet-glacier responses to ocean versus interior forcing over the Holocene as an additional factor for the modern state. We focus on climate anomalies that are well-documented 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, long-ago events compared to smaller, more recent events and (2) the glacier's relative memory of past ocean vs. interior forcing.

We use the two-stage 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 two-stage model responses. Anomalies are shown relative to the mid-Holocene 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 grounding-line 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 longer-term ocean cooling, via its two distinct response timescales.

Figure 7 Response to long-term climate variations. (a) An idealized climate scenario representing the last deglaciation, Little Ice Age (LIA), and anthropogenic forcing (black line). Normalized length responses of the linear model are shown for forcing applied via S (brown) and Ω (blue). The length at which the glacier would remain in equilibrium with climate is shown for reference (gray). (b) As for (a) but zoomed into the last 1000 yr. Ω forcing produces a much larger response on multicentury timescales, but the glacier does not reach equilibrium with the LIA in either case. (c) A more realistic forcing time series, based on reconstructed atmospheric temperatures for Disko Bay (black line; Buizert et al., 2018). Normalized length responses are again shown for both forcing types, which yield similar responses on multimillennial timescales. (d) As for (c) but zoomed into the last 1000 yr. Again, the LIA is expressed to a greater degree through ocean forcing. Note also the long-term advance driven by gradual cooling over the last several millennia.

4 Summary and discussion

Marine-terminating 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:

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

2. The two-stage 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 two-mode interpretation of outlet-glacier 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.

3. Persistence in stochastic climate forcing amplifies natural grounding-line 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).

4. 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 outlet-glacier 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 longer-timescale 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 short-term 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 ice-thickness 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 grounding-line-flux anomalies.

Another simplification is that we have focused on spatially uniform interior forcings, whereas in reality, surface-mass-balance 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 high-frequency 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 grounding-line response after 100 yr, while the 100 yr response following a step in Ω is 4 times greater. For white-noise 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 2-D 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 grounding-line-flux coefficient, which adjusts the parameterized relationship between ice thickness and grounding-line 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 grounding-line 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 up-glacier 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 grounding-zone 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 grounding-line 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 grounding-line 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 ice-sheet 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 ice-sheet models. The ISMIP6 ocean-forcing 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 ice-sheet 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 long-term terminus retreat compared with the melt implementation.

The impact that the choice of parameterization would have on projections of near-term sea-level rise is not immediately clear, however. Slater et al. (2020) note that these parameterizations are intended primarily for 21st-century 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 grounding-zone fluxes (much like in Fig. 5). In such a case, the choice of parameterization would affect the difference between responses to low-emission versus high-emission 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 sea-level rise by 2100, longer-timescale comparisons within the ISMIP6 framework, with both idealized and comprehensive models, may also prove illuminating. Understanding the timescales of glacier dynamics is a long-standing 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 marine-terminating 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 large-scale 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 ice-sheet 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.

Appendix A: Linearized two-stage model

Here we present the linearized two-stage 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:

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial H^{\prime }}{\partial t}}=A_H{H}^{\prime }+A_L{L}^{\prime }+{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{L}}}({\displaystyle \frac{\stackrel{\mathrm{‾}}{H}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}}-\mathrm{1})Q_{\mathrm{g}}^{\prime }+S^{\prime },\text{(A2)}& {\displaystyle \frac{\partial L^{\prime }}{\partial t}}=B_H{H}^{\prime }+B_L{L}^{\prime }-{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}}{Q}_{\mathrm{g}}^{\prime }.\end{array}$$

A_H, A_L, B_H, and B_L are the linear couplings between length and thickness anomalies:

$$\begin{array}{}\text{(A3)}& {\displaystyle }{A}_H(\stackrel{\mathrm{‾}}{H},\stackrel{\mathrm{‾}}{L})=& {\displaystyle }-\stackrel{\mathrm{‾}}{Q_{\text{g}}}\mathit{\alpha }{\stackrel{\mathrm{‾}}{h_{\text{g}}}}^{-\mathrm{1}}{\stackrel{\mathrm{‾}}{L}}^{-\mathrm{1}},\text{(A4)}& {\displaystyle }{A}_L(\stackrel{\mathrm{‾}}{H},\stackrel{\mathrm{‾}}{L})=& {\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\text{g}}}}{{\stackrel{\mathrm{‾}}{L}}^{\mathrm{2}}}}[\mathrm{1}+\mathit{\gamma }{\displaystyle \frac{\stackrel{\mathrm{‾}}{H}}{\stackrel{\mathrm{‾}}{h_{\text{g}}}}}+\mathit{\beta }\mathit{\lambda }{b}_x{\displaystyle \frac{\stackrel{\mathrm{‾}}{L}}{\stackrel{\mathrm{‾}}{h_{\text{g}}}}}(\mathrm{1}-{\displaystyle \frac{\stackrel{\mathrm{‾}}{H}}{\stackrel{\mathrm{‾}}{h_{\text{g}}}}}\left)\right],\text{(A5)}& {\displaystyle }{B}_H(\stackrel{\mathrm{‾}}{H},\stackrel{\mathrm{‾}}{L})=& {\displaystyle }\stackrel{\mathrm{‾}}{Q_{\text{g}}}\mathit{\alpha }{\stackrel{\mathrm{‾}}{H}}^{-\mathrm{1}}{\stackrel{\mathrm{‾}}{h_{\text{g}}}}^{-\mathrm{1}},\text{(A6)}& {\displaystyle }{B}_L(\stackrel{\mathrm{‾}}{H},\stackrel{\mathrm{‾}}{L})=& {\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\text{g}}}}{\stackrel{\mathrm{‾}}{h_{\text{g}}}}}(\mathit{\beta }\mathit{\lambda }{b}_x{\stackrel{\mathrm{‾}}{h_{\text{g}}}}^{-\mathrm{1}}-{\displaystyle \frac{\mathit{\gamma }}{\stackrel{\mathrm{‾}}{L}}}),\end{array}$$

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]:

$$\begin{array}{}\text{(A7)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{H_{\left[i\right]}^{\prime }-H_{[i-\mathrm{1}]}^{\prime }}{\mathrm{\Delta }t}}& =A_H{H}_{\left[i\right]}^{\prime }+A_L{L}_{\left[i\right]}^{\prime }\\ & +{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{L}}}\left({\displaystyle \frac{\stackrel{\mathrm{‾}}{H}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}}-\mathrm{1}\right)Q_{g\left[i\right]}^{\prime }+S_{\left[i\right]}^{\prime },\end{array}\text{(A8)}& {\displaystyle \frac{L_{\left[i\right]}^{\prime }-L_{[i-\mathrm{1}]}^{\prime }}{\mathrm{\Delta }t}}=B_H{H}_{\left[i\right]}^{\prime }+B_L{L}_{\left[i\right]}^{\prime }-{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}}{Q}_{g\left[i\right]}^{\prime }.\end{array}$$

Solving for H_[i] and L_[i] gives

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{H}_{\left[i\right]}^{\prime }& ={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-A_H\mathrm{\Delta }t}}\left[A_L\mathrm{\Delta }t{L}_{\left[i\right]}^{\prime }+H_{[i-\mathrm{1}]}^{\prime }\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{L}}}({\displaystyle \frac{\stackrel{\mathrm{‾}}{H}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}}-\mathrm{1})\mathrm{\Delta }t{Q}_{g\left[i\right]}^{\prime }+\mathrm{\Delta }t{S}_{\left[i\right]}^{\prime }\right],\end{array}\end{array}$$

$$\begin{array}{}\text{(A10)}& L_{\left[i\right]}^{\prime }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-B_L\mathrm{\Delta }t}}[B_H\mathrm{\Delta }t{H}_{\left[i\right]}^{\prime }+L_{[i-\mathrm{1}]}^{\prime }+{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}}\mathrm{\Delta }t{Q}_{g\left[i\right]}^{\prime }],\end{array}$$

which are still semi-implicit 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 two-dimensional autoregressive process (e.g., Box et al., 2008) with two forcing terms. For compactness, we redefine several parameter combinations at this point. Let $\mathit{ϵ}=(\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{ϵ}{A}_L{B}_H\mathrm{\Delta }{t}^{\mathrm{2}}{)}^{-\mathrm{1}}$. Then,

$$\begin{array}{}\text{(A11)}& \left[\begin{array}{c}{H}_{\left[i\right]}^{\prime }\\ L_{\left[i\right]}^{\prime }\end{array}\right]=\mathbf{\text{C}}\left[\begin{array}{c}{H}_{[i-\mathrm{1}]}^{\prime }\\ L_{[i-\mathrm{1}]}^{\prime }\end{array}\right]+\mathit{d}{Q_{\mathrm{g}}}_{\left[i\right]}^{\prime }+\mathit{e}{S}_{\left[i\right]}^{\prime },\end{array}$$

where the autoregressive coefficients are contained in the matrix C and vectors d and e:

$$\begin{array}{}\text{(A12)}& {\displaystyle }& {\displaystyle }\mathbf{\text{C}}=\mathit{\kappa }\left[\begin{array}{cc}\mathit{\eta }& \mathit{ϵ}\mathit{\eta }{A}_L\mathrm{\Delta }t\\ \mathit{ϵ}\mathit{\eta }{B}_H\mathrm{\Delta }t& \mathit{ϵ}\end{array}\right],\text{(A13)}& {\displaystyle }& {\displaystyle }\mathit{d}=\mathit{\kappa }\mathrm{\Delta }t\left[\begin{array}{c}\mathit{\eta }\left({\stackrel{\mathrm{‾}}{L}}^{-\mathrm{1}}\right(\frac{\stackrel{\mathrm{‾}}{H}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}-\mathrm{1})-\mathit{ϵ}{A}_L{h_{\mathrm{g}}}^{-\mathrm{1}}\mathrm{\Delta }t)\\ \mathit{ϵ}\left(B_H{\stackrel{\mathrm{‾}}{L}}^{-\mathrm{1}}\right(\frac{\stackrel{\mathrm{‾}}{H}}{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}-\mathrm{1})\mathrm{\Delta }t-{\stackrel{\mathrm{‾}}{h_{\mathrm{g}}}}^{-\mathrm{1}})\end{array}\right],\text{(A14)}& {\displaystyle }& {\displaystyle }\mathit{e}=\mathit{\kappa }\mathit{\eta }\mathrm{\Delta }t\left[\begin{array}{c}\mathrm{1}\\ \mathit{ϵ}{B}_H\mathrm{\Delta }t\end{array}\right].\end{array}$$

This form lends itself to a straightforward solution algorithm, given steady-state 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.