the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Biases in ice sheet models from missing noiseinduced drift
Vincent Verjans
Aminat A. Ambelorun
Most climatic and glaciological processes exhibit internal variability, which is omitted from many ice sheet model simulations. Prior studies have found that climatic variability can change ice sheet sensitivity to the longterm mean and trend in climate forcing. In this study, we use an ensemble of simulations with a stochastic largescale ice sheet model to demonstrate that variability in frontal ablation of marineterminating glaciers changes the mean state of the Greenland Ice Sheet through noiseinduced drift. Conversely, stochastic variability in surface mass balance does not appear to cause noiseinduced drift in these ensembles. We describe three potential causes for noiseinduced drift identified in prior statistical physics literature: noiseinduced bifurcations, multiplicative noise, and nonlinearities in noisy processes. Idealized simulations and Reynolds decomposition theory show that for marine ice sheets in particular, noiseinduced bifurcations and nonlinearities in variable ice sheet processes are likely the cause of the noiseinduced drift. We argue that the omnipresence of variability in climate and ice sheet systems means that the state of realworld ice sheets includes this tendency to drift. Thus, the lack of representation of such noiseinduced drift in spinup and transient ice sheet simulations is a potentially ubiquitous source of bias in ice sheet models.
 Article
(935 KB)  Fulltext XML
 BibTeX
 EndNote
The Earth system exhibits internal variability in many processes on a wide range of timescales. As one component of the Earth system, ice sheets are subject to variability in climatic processes, including snowfall, atmospheric temperatures, and ocean currents. Ice sheets also exhibit internal variability of their own in processes related to hydrology, ice fracture and ice flow. In general, numerical ice sheet modeling studies focus on the ice sheet response to changes in the mean forcing, often without including internal variability in climate or glaciological systems (e.g., Golledge et al., 2015; DeConto et al., 2021). The central assumption of such studies is that the longterm state of glaciers and ice sheets is set only by the multidecadal mean and trend in climate forcing. This assumption is based on the long equilibrium timescale of glaciers and ice sheets (Nye, 1960; Oerlemans and Van Der Veen, 1984; Robel et al., 2018). However, critically, this long response timescale does not imply that glaciers and ice sheets are insensitive to shorttimescale climatic fluctuations (Roe and O’Neal, 2009). Several recent studies, most using idealized glacier and ice sheet models, have demonstrated that this assumption may not hold in many circumstances known to exist in the real world. In landbased ice sheets with stochastic variability in surface temperature (Mikkelsen et al., 2018; Lauritzen et al., 2023) or marinebased ice sheets with periodic variability in ice viscosity (Hindmarsh and Le Meur, 2001), stochastic variability in ice shelf length (Robel et al., 2018), or seasonal variability in the calving front (Felikson et al., 2022), the inclusion of variability causes drift of the ice sheet state. This phenomenon of “noiseinduced drift” is well known in the statistical physics community, where many useful mathematical tools have been developed to understand the cause of this phenomenon (e.g., Kloeden and Platen, 1995; Horsthemke and Lefever, 1984).
In this study, we show that noiseinduced drift in response to stochastic frontal ablation is expected to occur in real marine ice sheets and numerical modeling of marine ice sheets. This is demonstrated with ensemble simulations of the Greenland Ice Sheet, resembling modern conditions with realistic stochastic variability in frontal ablation. Ensembles with stochastic forcing in surface mass balance do not exhibit the same noiseinduced drift, though other studies using stochastic surface temperature forcing in parameterized surface mass balance schemes do exhibit such drift (Mikkelsen et al., 2018). We describe the three different potential mechanisms of noiseinduced drift in generic stochastic systems and identify which of these mechanisms are likely to cause noiseinduced drift in real ice sheets. We close by arguing that modern ice sheet models omitting variability in climate and glaciological processes could produce biased estimates of the ice sheet mean state and the ice sheet response to climate change. We provide two potential solutions for this problem in the initialization and forcing of ice sheet models.
The central goal of this study is to demonstrate that the response of ice sheets to longterm (decadal–millennial) climatic forcing depends on the inclusion and magnitude of variability in climate and glaciological processes. To achieve this goal, we run four ensembles of Greenland Ice Sheet simulations using the Stochastic IceSheet and SeaLevel System Model (StISSM; Verjans et al., 2022). The core of this model is ISSM, which solves for the ice thickness and velocity on a finiteelement mesh refined in locations of interest (Larour et al., 2012). In this study, we use the shallow shelf approximation (SSA; MacAyeal, 1989) and refine the mesh at 11 large marineterminating glacier catchments, where the ice sheet margin evolves dynamically. All simulations are initialized at a deterministic steady state. This configuration is meant to resemble the modern state of the Greenland Ice Sheet but deviates somewhat from the real ice sheet which is not at a steady state (Otosaka et al., 2022). This initial deterministic steady state comes from a long spinup run over 31 000 years with temporally constant forcing in surface mass balance (SMB) and ablation at calving fronts (described in more detail in Verjans et al., 2022). SMB at model mesh elements is set according to an elevationdependent profile, which is fit separately in 19 catchments encompassing the entire ice sheet (Zwally et al., 2012), to resemble mean 1961–1990 SMB simulated in RACMO2 (Ultee et al., 2024; Ettema et al., 2009). Each marineterminating catchment has a prescribed rate of ocean melt at calving fronts based on thermal forcing from Wood et al. (2021). In the spinup, calving rates at each catchment are calibrated to produce a steadystate ice sheet configuration resembling the presentday ice sheet. We apply the Budd sliding law (Budd et al., 1979):
where τ_{b} is the basal friction, u_{b} is the basal sliding speed, and C^{2} is a spacevarying coefficient. Effective pressure, N, is set to maintain local hydrostatic equilibrium with the ocean throughout the icecovered model domain (Tsai et al., 2015):
where ρ_{i} and ρ_{w} are the densities of ice and water, respectively; g is the acceleration due to gravity; h is the ice thickness; and b is the bed elevation. Initialized from this steady state, a deterministic control run with temporally constant forcings exhibits an increase in ice mass of only 0.07 % in 2000 years. The spatial pattern of ice thickness change in this deterministic control run (not plotted) shows weak thickness changes which are uniformly distributed over catchments, indicating no significant changes to glacier termini.
We run ensembles of 10member simulations each, applying stochastic variability separately in SMB and calving rate, and we quantify the role of each forcing in setting the ice sheet state. Realistic stochastic parameterizations for SMB and ocean thermal forcing (which determines frontal melt) were described in previous studies (Ultee et al., 2024; Verjans et al., 2023). These studies found that variability in both SMB and ocean thermal forcing around Greenland is best described by autoregressive moving average models of low order. In this study, for ease of interpretability, we conservatively apply simple white noise to different forcing variables with a mean that remains constant in time and equal to deterministic steadystate values. White noise is characterized by independent random perturbations drawn from a Gaussian distribution and with no autocorrelation in time. For both stochastic ensembles, the standard deviation of the stochastic variable in each catchment is set to onethird of the mean in that catchment. This amplitude of variability is chosen for simplicity but is similar to variability from observations and highfidelity models of SMB and ocean forcing. In particular, Fettweis et al. (2020) find that averaged across 13 different SMB models calibrated against observations, Greenlandwide SMB has a temporal standard deviation, which is approximately 40 % of the mean. Hanna et al. (2011) develop observationbased reanalyses of Greenland SMB over the 20th century, which also indicate a temporal standard deviation which is approximately 25 %–35 % of the mean (depending on the calibration dataset used). Verjans et al. (2023) find that interannual variability in thermal forcing (which drives frontal ablation at glaciers) in the Estimating the Circulation and Climate of the Ocean (ECCO; Nguyen et al., 2012) Arctic reanalysis product typically ranges between 10 % and 60 % around Greenland. As a point of comparison, we also run a fourth ensemble with the standard deviation of the stochastic calving rate equal to a conservatively low onesixth of the mean calving rate.
In implementing white noise forcing in SMB and frontal ablation rate, we introduce symmetric variability directly in terms of the mass conservation equations for the ice sheet. This simplifies the task of identifying potential causes of resulting noiseinduced drift since the only dynamics to consider are those related to ice sheet flow. However, it may be that in reality, symmetric variability occurs in variables more removed from ice sheet dynamics such as atmospheric or ocean temperatures. Then, asymmetries or nonlinearities in the dependence of mass fluxes on these variables can be an additional source of noiseinduced drift, as previously discussed by Mikkelsen et al. (2018) and Lauritzen et al. (2023). Our goal in this study is to identify mechanisms of noiseinduced drift that are inherent in the fundamental dynamics of ice sheet flow. Such mechanisms would be common to all ice sheet models and not dependent on the modelspecific parameterizations of mass fluxes as a function of climate forcing.
Ensemble simulations are run for 2000 years in order to observe the ice sheet evolution towards a new state. However, we note that an ice sheet the size of Greenland likely requires more than 10 000 years to reach a new steady state in response to an icesheetwide change in forcing due to longterm dynamic adjustment extending through the interior. Such long simulations are computationally challenging to perform for the entire Greenland Ice Sheet on a wellresolved mesh. The design of this ensemble was initially inspired by the larger Greenland Ice Sheet ensemble used to benchmark StISSM in Verjans et al. (2022), which showed that just 10 ensemble members are sufficient to constrain the ensemblemean ice sheet mass to less than 0.1 % of the converged values (albeit under different stochastic forcing). We also note here that in this depthaveraged model, the dynamic influence of calving and ocean melt at glacier termini is identical. We have chosen to implement stochastic calving in this study, but the results would be identical if stochastic frontal melt were implemented instead.
Figure 1 shows the evolution of Greenland Ice Sheet mass over time from these ensemble simulations (colored lines and shading) in comparison to the deterministic control simulation (black line). The most striking result is that stochastic variability in calving at marineterminating glaciers causes substantial drift in the ensemblemean ice sheet mass (yellow and blue lines). This drift is apparent in all ensemble members and exceeds the spread of intraensemble variability after the first few years of the simulation (i.e., all ensemble members drift almost immediately). In the first 100 years of the simulation ensemble, the drift amounts to approximately 1 cm of global sea level equivalent, which is 5 %–10 % of the median projected Greenland contribution to sea level rise by 2100 in ISMIP6 (Goelzer et al., 2020). At the end of the 2000year simulation ensemble with highest variability amplitude (yellow line), the drift is larger than 1.5 % of total initial ice mass or about 12 cm of sea level equivalent. Based on these two ensembles, we conclude that the rate of drift increases with the amplitude of the variability in calving rate. As a point of comparison, the dashed line shows a single simulation, without stochastic variability but with a 270 % increase in the mean calving rate at all 11 marineterminating glaciers for which we simulate terminus migration. The spatial pattern of ice thickness change in this simulation (not plotted) is very similar to the stochastic calving ensemble with highest variability amplitude, indicating that the noiseinduced drift in the stochastic ensemble occurs due to increased mass loss at the terminus. This indicates that model drift due to a realistic level of noise in just the annual calving rate is equivalent to ice loss from a substantial increase in calving rate without noise. Calibrating a deterministic model to match the observed ice sheet state, which is subject to variability from climatic and glaciological processes, would require tuning parameters to very different values. We discuss the resulting biases in Sect. 4.
Variability in SMB (green line) does not drive discernible drift in the ice sheet volume, in contrast to the study of Lauritzen et al. (2023), which found strong noiseinduced drift in an ensemble of Greenland Ice Sheet simulations in response to temperature variability applied through a positivedegreeday model. We do not use such a model to parameterize SMB. Instead, we specify stochastic variability directly in SMB on a catchmentbycatchment basis.
While these stochastic ensembles exhibit less than 2 % changes in their total Greenland Ice Sheet mass after 2000 years, the local change in ice thickness at some of the largest marineterminating glaciers in Greenland is a substantial fraction of their initial ice thickness (Fig. 2c). At some glaciers, there is thinning in some ensemble members and thickening at others. At other glaciers, all ensemble members show thinning. To show the expression of this noiseinduced drift at different glaciers, we further plot profiles of ice thickness for all ensemble members at Sermeq Kujalleq (also called Jakobshavn Isbræ) in Fig. 2a–b and Petermann Glacier in Fig. 2d–e. Under a sufficiently large amplitude of variability in calving rate, retreat of the terminus of Sermeq Kujalleq occurs episodically with timing that is variable across ensemble members (Fig. 2a). At Petermann Glacier, retreat of the terminus is monotonic and nearly uniform across ensemble members during the early parts of simulations (Fig. 2d–e). The different expressions of this drift indicate that there is likely to be more than one mechanism responsible for producing the drift, as explored in the next section.
Many systems, including the climate system (Penland, 2003), exhibit noiseinduced drift, wherein inclusion of noise causes a change in the mean system state. To explain the potential causes of noiseinduced drift, we start from a generic stochastic differential equation:
where x is the model state, t is time, f(x) is a function describing the deterministic model dynamics, g(x) is a function describing how the amplitude of the noise forcing the system may depend on model state, η(t) is a random noise term drawn from some distribution (typically Gaussian), and β is an exponent. For the sake of simplicity, we treat Eq. (3) in a scalar form, but it can be generalized to a vectorvalued case without loss of generality. In the case where $f\left(x\right)=\mathit{\alpha}x$, g(x)=1, β=1, and η(t) is a random variable drawn from a Gaussian distribution, this is the Langevin equation describing Brownian motion of a particle without drift. However, in many more complex systems, real physical processes described by the components of this equation lead to noiseinduced drift. For a more technical review of noiseinduced drift, the interested reader is referred to Horsthemke and Lefever (1984).
Here, we describe three causes of noiseinduced drift that are potentially relevant to ice sheets:

Noiseinduced bifurcation/tipping. In Eq. (3), when f(x)=αx, α describes the linear stability of the system. If α is negative, the system is stable as perturbations from the noise term η(t) are damped. If α is positive, the system is unstable as perturbations from the noise term η(t) are not damped. Thus, if a noise perturbation causes α to change sign (i.e., a bifurcation), the system undergoes a transition to a different state. Such stability properties have been previously explored in the context of ice sheet dynamics, where loss of ice sheet stability through marine ice sheet instability or other bifurcations may be caused by variability in climate forcing (Mulder et al., 2018; Christian et al., 2022; Sergienko and Haseloff, 2023).

Multiplicative noise. In Eq. (3), when g(x) is any function that is not even about the fixed point x_{*} ($\frac{\partial f}{\partial x}{}_{x={x}_{*}}=\mathrm{0}$), i.e., $g({x}_{*}\mathit{\eta})\ne g({x}_{*}+\mathit{\eta})$, this describes any system where the amplitude of noise perturbations depends on the system state, causing the entire noise term g(x)η(t) to have a nonzero mean. Physically, such multiplicative noise arises in systems where there is noise in a term that depends on system state. This has previously been explored in the context of simple glacier models (Mantelli et al., 2016; Robel et al., 2018; Mikkelsen et al., 2018).

Nonlinear or asymmetric noise. If β≠1 (excluding the trivial case where β=0) or if the underlying noise process has a nonzero mean (i.e., the distribution of noise is intrinsically asymmetric), then the noise term will cause drift in the mean system state. Because most canonical stochastic models assume that the noise term is linear and sampled from a Gaussian distribution, this potential cause of noiseinduced drift has received considerably less attention in the literature (although it is discussed in detail by Horsthemke and Lefever, 1984). Glacier ice is a viscous nonNewtonian fluid, meaning that glacier flow speed exhibits a strong nonlinear sensitivity to many different types of forcing (Glen, 1955; Millstein et al., 2022). Robel et al. (2018) previously considered this source of noiseinduced drift in the context of ice shelf buttressing, but many other processes related to ice flow may exhibit similar nonlinear noiseinduced drift.
To understand the role of these different potential causes of noiseinduced drift in ice sheet dynamics, we consider several highly idealized stochastic ensembles. In each simulation, we use StISSM to simulate ice velocity and thickness evolution of a single marineterminating glacier in a rectangular channel of uniform width, without floating ice. Model configuration choices such as the stress balance approximation and the basal sliding law are identical to the Greenland ensemble described in the previous section but with a spatially uniform basal friction coefficient (C^{2}). In all configurations, an initial deterministic steady state is obtained by holding all forcing variables constant and running the simulation until the total ice mass of the glacier changes by less than 0.05 % in 200 years. In each idealized stochastic ensemble, calving rate is drawn from a Gaussian distribution (i.e., white noise) with a mean equivalent to the initial deterministic calving rate and standard deviation equal to onethird of the mean. We perform ensemble simulations of 30 members each, running for 2000 years.
3.1 Noiseinduced bifurcation/tipping
Figure 3 shows the results of three idealized stochastic ensembles, all of which have the same background prograde slope of 0.004 in bed topography. In the first stochastic ensemble (Fig. 3a–b), the bed topography includes a single sinusoidal bump 100 m in height in bed topography at the initial terminus position. Once stochastic calving begins, 95 % of the ensemble members start retreating past the bump within the first 140 years of the simulation. The second ensemble (Fig. 3c–d) is identical to the first, except without a bump in bed topography, and the steadystate calving rate used in the spinup is adjusted to maintain a similar terminus position. Though the initial glacier state is not identical due to the difference in bed topography, it is sufficiently similar for us not to attribute the subsequent behavior to a different glacier state. Instead of retreating, all ensemble members advance in response to stochastic calving. The different response to stochastic forcing between these two ensembles indicates that the ensemblemean retreat in the first ensemble is caused by the presence of the bump in bed topography, which adds a wellunderstood bifurcation to the system dynamics related to a positive feedback in ice flow with bed depth. This provides a simple example of mechanism no. 1 identified above, i.e., noiseinduced bifurcation/tipping.
When a noiseinduced bifurcation drives drift in the mean state, the rate of drift depends on the amplitude of stochastic variability up to the amplitude of variability necessary to drive all ensemble members across the bifurcation with high probability. In Fig. 3a–b, this “saturation” of the drift rate is occurring as all ensemble members eventually cross the bifurcation. Further increasing the amplitude of the variability will not be able to drive more ensemble members through the bifurcation, though they might reach it faster near the beginning of the simulation, causing faster initial drift of the mean state. As a point of comparison, in the full Greenland Ice Sheet ensemble discussed in Sect. 2, the magnitude of ensemblemean retreat at Sermeq Kujalleq shows a clear dependence between low (Fig. 2b) and high (Fig. 2a) amplitude of calving variability. In this case, further increasing the amplitude of variability may cause some ensemble members to retreat past the second bed peak, thus increasing the extent of noiseinduced drift.
3.2 Multiplicative noise
Noiseinduced tipping is clearly not the only mechanism causing the drift seen in the more realistic simulations discussed in the prior section since drift still occurs even in the absence of a bifurcation in system dynamics. Multiplicative noise (mechanism no. 2) may explain this drift in the second stochastic ensemble as variability at the calving front perturbs the nearterminus thickness, causing variations in effective pressure and ultimately velocity through the Budd sliding law (Eqs. 1–2). This particular sliding law includes a linear dependence of basal friction on effective pressure and therefore ice thickness although there are other nonlinearities elsewhere, which may play a role in generating drift. Since the variable that is being perturbed stochastically is linearly related to ice flow and the nonlinearities arise elsewhere in the ice sheet dynamical equations, this is considered to be multiplicative noise similar to g(x) being multiplied by η(t) in Eq. (3). To investigate this possibility, we consider a stochastic ensemble (Fig. 3e–f) in which the effective pressure dependence is removed from Eq. (1), effectively introducing a sliding law linear in sliding velocity only. In this case, drift still occurs, indicating that multiplicative noise through evolving effective pressure is unlikely to be the only mechanism causing the drift. Though ice sheet dynamics involve the complex interplay of many factors, the lack of other obvious multiplicative feedbacks likely to cause a significant asymmetry in the variability in terminus thickness or velocity strongly indicates the drift seen in these two ensembles is mainly caused by a different mechanism, i.e., nonlinear noise (mechanism no. 3 above).
3.3 Nonlinear noise
Though there are many sources of nonlinearity in ice sheet dynamics, the fact that only stochasticity in calving causes drift in the Greenland ensemble of the previous section indicates that it is some nonlinear process specific to the glacier terminus which leads to noiseinduced drift in the absence of a bifurcation. Here, we give mathematical explanations for the drift in response to stochastic variability in the terminus position, which applies to tidewater glaciers and glaciers with floating ice shelves.
In all stochastic simulations considered in this study, the mean of the rate of calving at the terminus (u_{c}) does not change, and so any changes in the timeaveraged terminus position must be the result of changes in mean ice flow velocity towards the terminus (u_{f}). For a tidewater glacier, like that simulated in Fig. 3, u_{f} is determined by the momentum balance at the terminus:
where h is the terminus thickness, b is the water depth, ρ_{i} is the ice density, $\mathit{\lambda}=\frac{{\mathit{\rho}}_{\mathrm{w}}}{{\mathit{\rho}}_{\mathrm{i}}}$ is the ratio of water to ice density, g is the gravitational acceleration, A is the depthintegrated Glen's flow law rate factor, and n is the Glen's flow law exponent. Perturbations to the mean terminus position may cause perturbations to the glacier thickness and bed depth at the terminus which can be included through a Reynolds decomposition ($h=\langle h\rangle +{h}^{\prime}$ and $b=\langle b\rangle +{b}^{\prime}$), where all variables enclosed by 〈〉 are timeaveraged and perturbed variables are denoted by ^{′}. All perturbed variables are drawn from a Gaussian distribution with a zero mean. Including these decomposed expressions into the above momentum balance and simplifying yield an expression for the strain rate at the terminus in terms of perturbations.
The quadratic term in this expression is expanded, and we separate terms with only the mean state in their numerator from those including perturbations in their numerator.
We perform a Taylor expansion on the resulting expression in terms of the exponent n, keeping in mind that terms involving perturbations will generally be smaller than terms involving only the mean state. Thus, terms depending on higher powers of h^{′} and b^{′} can be neglected, and we only keep the first two terms of the expansion (i.e., linearize).
We rearrange this expression to emphasize the relative influences of the mean state and perturbations.
To understand the effect of perturbations on the glacier mean state, we take a time average of this expression, which eliminates terms that are linear in a perturbation variable because they have a mean of zero.
Note that in the above step, terms which include perturbations as a sum in the denominator are linearized through a Taylor series expansion before the average is taken, leaving only the terms involving the mean state, 〈h〉. If the perturbation terms are drawn from a Gaussian distribution with variance σ^{2}, then terms involving the square of the perturbation are drawn from a gamma distribution, $\mathrm{\Gamma}\left(\frac{\mathrm{1}}{\mathrm{2}},\mathrm{2}{\mathit{\sigma}}^{\mathrm{2}}\right)$, which has a nonzero mean equal to σ^{2}. Thus, the rate of drift depends on how large $n\mathit{\lambda}{\mathit{\sigma}}_{\mathrm{b}}^{\mathrm{2}}$ is relative to 〈h〉^{2}−λ〈b〉^{2}, where ${\mathit{\sigma}}_{\mathrm{b}}^{\mathrm{2}}$ is the variance of the perturbations in bed depth due to perturbations in the ice front position. The sign of this leadingorder term causing the drift is negative, causing a decrease in the nearterminus strain rate and a net positive mass balance near the terminus, driving advance. While we might expect that σ_{b}≪〈b〉, if the bed topography (b_{x}) is steep, then σ_{b}=b_{x}σ_{L} (where σ_{L} is the standard deviation of variability in terminus position) could be a nonnegligible fraction of 〈b〉, causing appreciable drift. Also, if the terminus is at or near flotation, then $\langle h{\rangle}^{\mathrm{2}}\mathit{\lambda}\langle b{\rangle}^{\mathrm{2}}\approx {\mathit{\lambda}}^{\mathrm{2}}\langle b{\rangle}^{\mathrm{2}}\mathit{\lambda}\langle b{\rangle}^{\mathrm{2}}\approx \mathrm{0.1}\langle b{\rangle}^{\mathrm{2}}$ and the denominator of the above expression would be sufficiently small to admit nonnegligible drift. The simulations in Fig. 3c–f do exhibit such thickening and advance of the initially grounded terminus. Given that both steep bed topography and nearflotation termini are common in Greenland, we may expect this effect to be common, though we do not simulate any cases of ensemblemean glacier advance in the more realistic Greenland ensemble (Fig. 3c).
For a glacier with a floating ice shelf, the calving front is not grounded and so the momentum balance does not depend on the bed depth, making the above analysis not applicable. We rather consider the effect of buttressing from the floating ice shelf on the velocity of ice through the grounding line. Haseloff and Sergienko (2018) perform an asymptotic analysis to derive an approximation for the ice flow velocity (u_{g}) through a strongly buttressed grounding line:
where Λ is a parameter governing lateral shear stress within the ice and L_{s} is the ice shelf length. This expression assumes that ice loss occurs entirely through ablation at the calving front and lateral shear stresses increase linearly across the ice shelf. We consider stochastic calving at the calving front of the floating ice, causing Gaussian, zeromean perturbations to the ice shelf length (${L}_{\mathrm{s}}=\langle {L}_{\mathrm{s}}\rangle +{L}_{\mathrm{s}}^{\prime}$). We insert this Reynolds decomposition in the above expression for grounding line flux and take the Taylor expansion of the resulting expression to get
From this expression, we neglect higherorder terms and rearrange to resemble the original flux expression:
Taking the time average, the term which is linear in ${L}_{\mathrm{s}}^{\prime}$ vanishes, leaving
The ${{L}_{\mathrm{s}}^{\prime}}^{\mathrm{2}}$ term is drawn from the $\mathrm{\Gamma}\left(\frac{\mathrm{1}}{\mathrm{2}},\mathrm{2}{\mathit{\sigma}}_{{L}_{\mathrm{s}}}^{\mathrm{2}}\right)$ distribution, which has a nonzero mean equal to ${\mathit{\sigma}}_{{L}_{\mathrm{s}}}^{\mathrm{2}}$. Thus, when $\frac{n(n+\mathrm{1}){\mathit{\sigma}}_{{L}_{\mathrm{s}}}^{\mathrm{2}}}{\mathrm{2}}$ is nonnegligible compared to 〈L_{s}〉^{2}, the timeaveraged ice flow velocity through the grounding line is increased by stochastic calving, which causes the grounding line to retreat. We can note here that different assumptions can be made about the form of lateral shear stress variation across the floating ice shelf or the dominant source of mass loss and, in general, that the rate of ice flow through the grounding line will be nonlinear in terms of the ice shelf length (Haseloff and Sergienko, 2018), causing the sort of nonlinear noiseinduced drift discussed here.
3.4 Attributing causes of drift
Returning to the more realistic Greenland Ice Sheet ensemble simulations (Fig. 2a–b, d–e), we conclude that in most glaciers for which strong noiseinduced drift is simulated, there are easily identifiable bed topography features indicating that noiseinduced bifurcations are the most common cause of noiseinduced drift (as previously argued in Christian et al., 2022). Conversely, there are no tidewater glaciers in this realistic ensemble exhibiting ensembleaverage terminus advance due to nonlinearity in hydrostatic stress terms discussed in the previous section. This is likely because glaciers tend to stabilize at peaks in bed topography (Robel et al., 2022), making it more likely that the sudden onset of stochastic calving would lead to a retreat from noiseinduced bifurcation rather than sustained advance due to the nonlinear noise mechanism. In contrast, during the earliest stage of Petermann Glacier's retreat (Fig. 2e), the bed is entirely prograde and yet ensemblemean retreat still occurs. At the time of our study, Petermann Glacier is one of only two glaciers in Greenland with a buttressing ice shelf remaining. Thus, the mechanism of drift due to nonlinearities in buttressing, discussed in the previous section, is likely responsible for the early stages of strong retreat of the Petermann grounding line before reaching a bed peak after which a noiseinduced bifurcation over a bed peak likely also plays an important role in the simulated ensemblemean retreat.
We also briefly note that Lauritzen et al. (2023) find that variability in surface temperature can cause noiseinduced drift through a positive degreeday (PDD) model for SMB, though they do not speculate on the cause of this drift (or refer to it as such). It is likely that the strong nonlinearities in their PDD model are the cause of the noiseinduced drift they find in their results, as their simulations do not appear to include bifurcations in SMB or sources of multiplicative noise. Regardless of the precise mechanism of noiseinduced drift in different model configurations, our simulations show that there are a range of different mechanisms intrinsic to ice sheet dynamics that cause noiseinduced drift to be an expected and essential aspect of ice sheet evolution and therefore of realistic model simulations. We purposely adopt a conservative approach to applying stochastic forcing directly to terms in the mass conservation equations of the ice sheet model, but we expect that stochastic variability in climatic and glaciological processes drives noiseinduced drift through many different mechanisms in more realistic ice sheet simulations.
The Greenland ensemble simulations in this study exhibit a tendency for noiseinduced retreat and ice loss. Thus, the spinup of an ice sheet model without variability in forcing is likely to lead to a modeled ice sheet that is biased compared to observations of real ice sheets, which are naturally subject to variable forcing and resulting noiseinduced drift. Such a mismatch is typically reduced by tuning or optimizing model parameter values, including those related to ice sliding, viscosity, calving, and ocean melt, through inversion. However, calibrating a parameter to minimize model–observation mismatch arising due to processes not represented in the model may introduce compensating errors in the modeled state. Ice sheet models that tune one parameter to reduce biases in other parameters have been shown to have substantially biased sensitivity to changes in forcing (Berends et al., 2023).
Many contemporary projections of future ice sheet evolution omit variability in forcing for transient projections due to challenges related to modeling ocean circulation near ice sheets or the lack of output from climate models far into the future (Golledge et al., 2015; DeConto et al., 2021). Such an omission may lead the modeled ice sheet sensitivity to future changes to be biased, as noiseinduced retreat is an important and realistic component of the forced response. As discussed in the prior section and prior studies (Christian et al., 2022), in the absence of variability, many glaciers may not cross important thresholds to rapid retreat, and thus their projected response to climate forcing would be considerably less than is likely in reality. Additionally, potential future changes in the amplitude of variability (e.g., Bintanja et al., 2020) could increase the likelihood of crossing noiseinduced bifurcations and amplify the impacts of statedependent and nonlinear noise. Such effects cannot be captured if variability in forcing is omitted entirely.
Other contemporary projections of future ice sheet evolution (e.g., many of the models participating in the recent ISMIP6 intercomparisons; Goelzer et al., 2020; Seroussi et al., 2020) start from a calibrated initial state and then simulate the freerunning ice sheet state in response to forcing including variability. In such a simulation design, the sudden onset of variability could introduce a transient noiseinduced drift. If the drift causes ice loss, as in the ensembles described in Sect. 2, this would cause the projected ice sheet sensitivity to forcing to be too high. Other recent modeling studies use a calibrated initial state but then recalibrate the ice sheet sensitivity to a changing mean climate with historical observations of ice sheet change (e.g., Nias et al., 2019; DeConto et al., 2021). In such a case, the calibrated sensitivity to changes in the mean climate would be too low due to the spurious influence of noiseinduced drift following the sudden onset of variability in the model. Similarly, the practice of removing control simulations, with forcing held constant to diagnose ice sheet sensitivity to forcing (Seroussi et al., 2020; Goelzer et al., 2020), may introduce bias due to the lack of noiseinduced drift in control simulations.
Noiseinduced drift in ice sheets should not only be thought of as a source of bias in models. Real ice sheets are subject to stochastic variability in many processes, thus meaning that their state, whether steady or not, includes the effect of noiseinduced drift. The potential ice sheet model biases identified here all result from an incomplete representation of these real sources of variability within climate or glaciological processes. To eliminate or lessen these biases in ice sheet models, we recommend two possible solutions for initializing ice sheets model simulations: (1) initializing directly from the observed ice sheet state without relaxation, even when the ice sheet is out of balance or (2) including internal variability in the forcing of ice sheet models during spinup. The first proposed solution recognizes that the observed state of ice sheets in the real world, subject to variability, should implicitly include the tendency resulting from noiseinduced drift. Ice sheet modelers may prefer using such a solution as it requires less computational resources; however, data assimilation methods for accurately reproducing observed nonsteady ice sheet states are still a nascent area of development (Goldberg and Heimbach, 2013; Choi et al., 2023). The second suggested solution is likely to be necessary if an initial steady state for a simulation is desired and observations of ice sheet state and tendency are not available, as in most simulations starting prior to the satellite era. Improving both glaciological process models (e.g., hydrology and calving) and the efficiency of coupling to climate models should also yield improvements in the complete and accurate representation of variability. Finally, stochastic ice sheet modeling (e.g., StISSM; Verjans et al., 2022) provides a parallel approach to accurately including variability within ice sheet models in a computationally efficient manner.
StISSM is an opensource largescale stochastic ice sheet model that is currently included as part of the public release of ISSM. The public SVN repository for the ISSM code can be found at https://issm.ess.uci.edu/svn/issm/issm/trunk (Larour et al., 2024) and downloaded using the username “anon” and password “anon”. The documentation of the code version used here is available at https://issm.jpl.nasa.gov/documentation/ (last access: 30 October 2023). Scripts for the Greenland spinup simulation follow Verjans et al. (2022) and are available at https://doi.org/10.5281/zenodo.7347470 (Verjans, 2022).
AAR conceived the study, developed the theory, and led the writing of the paper. VV performed the model simulations, postprocessed the model output, and contributed to writing. AAA contributed to the mathematical analysis and writing.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We thank John Erich Christian and Judith Berner for fruitful discussions during the conception of this study. We thank Johanna Beckmann for pointing out an error in the original preprint.
This research has been supported by the HeisingSimons Foundation (grant no. 20201965).
This paper was edited by Alexander Robinson and reviewed by two anonymous referees.
Berends, C. J., van de Wal, R. S. W., van den Akker, T., and Lipscomb, W. H.: Compensating errors in inversions for subglacial bed roughness: same steady state, different dynamic response, The Cryosphere, 17, 1585–1600, https://doi.org/10.5194/tc1715852023, 2023. a
Bintanja, R., van der Wiel, K., Van der Linden, E., Reusen, J., Bogerd, L., Krikken, F., and Selten, F.: Strong future increases in Arctic precipitation variability linked to poleward moisture transport, Science Advances, 6, eaax6869, https://doi.org/10.1126/sciadv.aax6869, 2020. a
Budd, W., Keage, P., and Blundy, N.: Empirical studies of ice sliding, J. Glaciol., 23, 157–170, 1979. a
Choi, Y., Seroussi, H., Morlighem, M., Schlegel, N.J., and Gardner, A.: Impact of timedependent data assimilation on ice flow model initialization and projections: a case study of Kjer Glacier, Greenland, The Cryosphere, 17, 5499–5517, https://doi.org/10.5194/tc1754992023, 2023. a
Christian, J. E., Robel, A. A., and Catania, G.: A probabilistic framework for quantifying the role of anthropogenic climate change in marineterminating glacier retreats, The Cryosphere, 16, 2725–2743, https://doi.org/10.5194/tc1627252022, 2022. a, b, c
DeConto, R. M., Pollard, D., Alley, R. B., Velicogna, I., Gasson, E., Gomez, N., Sadai, S., Condron, A., Gilford, D. M., Ashe, E. L., Kopp, R. E., Li, D., and Dutton, A.: The Paris Climate Agreement and future sealevel rise from Antarctica, Nature, 593, 83–89, 2021. a, b, c
Ettema, J., van den Broeke, M. R., van Meijgaard, E., van de Berg, W. J., Bamber, J. L., Box, J. E., and Bales, R. C.: Higher surface mass balance of the Greenland ice sheet revealed by highresolution climate modeling, Geophys. Res. Lett., 36, L12501, https://doi.org/10.1029/2009GL038110, 2009. a
Felikson, D., Nowicki, S., Nias, I., Morlighem, M., and Seroussi, H.: Seasonal Tidewater Glacier Terminus Oscillations Bias MultiDecadal Projections of Ice Mass Change, J. Geophys. Res.Earth, 127, e2021JF006249, https://doi.org/10.1029/2021JF006249, 2022. a
Fettweis, X., Hofer, S., KrebsKanzow, U., Amory, C., Aoki, T., Berends, C. J., Born, A., Box, J. E., Delhasse, A., Fujita, K., Gierz, P., Goelzer, H., Hanna, E., Hashimoto, A., Huybrechts, P., Kapsch, M.L., King, M. D., Kittel, C., Lang, C., Langen, P. L., Lenaerts, J. T. M., Liston, G. E., Lohmann, G., Mernild, S. H., Mikolajewicz, U., Modali, K., Mottram, R. H., Niwano, M., Noël, B., Ryan, J. C., Smith, A., Streffing, J., Tedesco, M., van de Berg, W. J., van den Broeke, M., van de Wal, R. S. W., van Kampenhout, L., Wilton, D., Wouters, B., Ziemen, F., and Zolles, T.: GrSMBMIP: intercomparison of the modelled 1980–2012 surface mass balance over the Greenland Ice Sheet, The Cryosphere, 14, 3935–3958, https://doi.org/10.5194/tc1439352020, 2020. a
Glen, J. W.: The creep of polycrystalline ice, P. R. Soc. AMath. Phy., 228, 519–538, 1955. a
Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., AbeOuchi, A., Shepherd, A., Simon, E., Agosta, C., Alexander, P., Aschwanden, A., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Golledge, N. R., Greve, R., Humbert, A., Huybrechts, P., Le clec'h, S., Lee, V., Leguy, G., Little, C., Lowry, D. P., Morlighem, M., Nias, I., Quiquet, A., Rückamp, M., Schlegel, N.J., Slater, D. A., Smith, R. S., Straneo, F., Tarasov, L., van de Wal, R., and van den Broeke, M.: The future sealevel contribution of the Greenland ice sheet: a multimodel ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096, https://doi.org/10.5194/tc1430712020, 2020. a, b, c
Goldberg, D. N. and Heimbach, P.: Parameter and state estimation with a timedependent adjoint marine ice sheet model, The Cryosphere, 7, 1659–1678, https://doi.org/10.5194/tc716592013, 2013. a
Golledge, N. R., Kowalewski, D. E., Naish, T. R., Levy, R. H., Fogwill, C. J., and Gasson, E. G.: The multimillennial Antarctic commitment to future sealevel rise, Nature, 526, 421–425, 2015. a, b
Hanna, E., Huybrechts, P., Cappelen, J., Steffen, K., Bales, R. C., Burgess, E., McConnell, J. R., Peder Steffensen, J., Van den Broeke, M., Wake, L., Bigg, G., Griffiths, M., and Savas, D.: Greenland Ice Sheet surface mass balance 1870 to 2010 based on Twentieth Century Reanalysis, and links with global climate forcing, J. Geophys. Res.Atmos., 116, D24121, https://doi.org/10.1029/2011JD016387, 2011. a
Haseloff, M. and Sergienko, O. V.: The effect of buttressing on grounding line dynamics, J. Glaciol., 64, 417–431, 2018. a, b
Hindmarsh, R. C. and Le Meur, E.: Dynamical processes involved in the retreat of marine ice sheets, J. Glaciol., 47, 271–282, 2001. a
Horsthemke, W. and Lefever, R.: Noiseinduced transitions in physics, chemistry, and biology, Springer, ISBN 9780387113593, 1984. a, b, c
Kloeden, P. E. and Platen, E.: Numerical Solutions of Stochastic Differential Equations, ISBN 9783642081071, 1995. a
Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res.Earth, 117, F01022, https://doi.org/10.1029/2011JF002140, 2012. a
Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Ice sheet and Sea Level System Model, ISSM [data set], https://issm.ess.uci.edu/svn/issm/issm/trunk, last access: 24 May 2024.
Lauritzen, M., Aðalgeirsdóttir, Guðfinna, R. N., Grinsted, A., Noel, B., and Hvidberg, C. S.: The influence of interannual temperature variability on the Greenland Ice Sheet volume, Ann. Glaciol., 1–8, https://doi.org/10.1017/aog.2023.53, 2023. a, b, c, d
MacAyeal, D. R.: Largescale ice flow over a viscous basal sediment: Theory and application to ice stream B, Antarctica, J. Geophys. Res.Sol. Ea., 94, 4071–4087, 1989. a
Mantelli, E., Bertagni, M. B., and Ridolfi, L.: Stochastic ice stream dynamics, P. Natl. Acad. Sci. USA, 113, E4594–E4600, 2016. a
Mikkelsen, T. B., Grinsted, A., and Ditlevsen, P.: Influence of temperature fluctuations on equilibrium ice sheet volume, The Cryosphere, 12, 39–47, https://doi.org/10.5194/tc12392018, 2018. a, b, c, d
Millstein, J. D., Minchew, B. M., and Pegler, S. S.: Ice viscosity is more sensitive to stress than commonly assumed, Commun. Earth Environ., 3, 57, https://doi.org/10.1038/s4324702200385x, 2022. a
Mulder, T., Baars, S., Wubs, F., and Dijkstra, H.: Stochastic marine ice sheet variability, J. Fluid Mech., 843, 748–777, 2018. a
Nguyen, A. T., Kwok, R., and Menemenlis, D.: Source and pathway of the Western Arctic upper halocline in a dataconstrained coupled ocean and sea ice model, J. Phys. Oceanogr., 42, 802–823, 2012. a
Nias, I. J., Cornford, S. L., Edwards, T. L., Gourmelen, N., and Payne, A. J.: Assessing uncertainty in the dynamical ice response to ocean warming in the Amundsen Sea Embayment, West Antarctica, Geophys. Res. Lett., 46, 11253–11260, 2019. a
Nye, J. F.: The response of glaciers and icesheets to seasonal and climatic changes, P. R. Soc. Lond. A Mat., 256, 559–584, 1960. a
Oerlemans, J. and Van Der Veen, C. J.: Ice sheets and climate, Vol. 21, Springer, https://doi.org/10.1007/9789400963252, 1984. a
Otosaka, I. N., Shepherd, A., Ivins, E. R., Schlegel, N.J., Amory, C., van den Broeke, M. R., Horwath, M., Joughin, I., King, M. D., Krinner, G., Nowicki, S., Payne, A. J., Rignot, E., Scambos, T., Simon, K. M., Smith, B. E., Sørensen, L. S., Velicogna, I., Whitehouse, P. L., A, G., Agosta, C., Ahlstrøm, A. P., Blazquez, A., Colgan, W., Engdahl, M. E., Fettweis, X., Forsberg, R., Gallée, H., Gardner, A., Gilbert, L., Gourmelen, N., Groh, A., Gunter, B. C., Harig, C., Helm, V., Khan, S. A., Kittel, C., Konrad, H., Langen, P. L., Lecavalier, B. S., Liang, C.C., Loomis, B. D., McMillan, M., Melini, D., Mernild, S. H., Mottram, R., Mouginot, J., Nilsson, J., Noël, B., Pattle, M. E., Peltier, W. R., Pie, N., Roca, M., Sasgen, I., Save, H. V., Seo, K.W., Scheuchl, B., Schrama, E. J. O., Schröder, L., Simonsen, S. B., Slater, T., Spada, G., Sutterley, T. C., Vishwakarma, B. D., van Wessem, J. M., Wiese, D., van der Wal, W., and Wouters, B.: Mass balance of the Greenland and Antarctic ice sheets from 1992 to 2020, Earth Syst. Sci. Data, 15, 1597–1616, https://doi.org/10.5194/essd1515972023, 2023. a
Penland, C.: Noise out of chaos and why it won't go away, B. Am. Meteorol. Soc., 84, 921–925, 2003. 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, 123, 2205–2227, 2018. a, b, c, d
Robel, A. A., Pegler, S. S., Catania, G., Felikson, D., and Simkins, L. M.: Ambiguous stability of glaciers at bed peaks, J. Glaciol., 68, 1177–1184, 2022. a
Roe, G. H. and O’Neal, M. A.: The response of glaciers to intrinsic climate variability: observations and models of lateHolocene variations in the Pacific Northwest, J. Glaciol., 55, 839–854, https://doi.org/10.3189/002214309790152438, 2009. a
Sergienko, O. and Haseloff, M.: “Stable” and “unstable” are not useful descriptions of marine ice sheets in the Earth's climate system, J. Glaciol., 69, 1483–1499, https://doi.org/10.1017/jog.2023.40, 2023. a
Seroussi, H., Nowicki, S., Payne, A. J., Goelzer, H., Lipscomb, W. H., AbeOuchi, A., Agosta, C., Albrecht, T., AsayDavis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., GaltonFenzi, B. K., Gladstone, R., Golledge, N. R., Gregory, J. M., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huybrechts, P., Jourdain, N. C., Kleiner, T., Larour, E., Leguy, G. R., Lowry, D. P., Little, C. M., Morlighem, M., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Reese, R., Schlegel, N.J., Shepherd, A., Simon, E., Smith, R. S., Straneo, F., Sun, S., Trusel, L. D., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., Zhao, C., Zhang, T., and Zwinger, T.: ISMIP6 Antarctica: a multimodel ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070, https://doi.org/10.5194/tc1430332020, 2020. a, b
Tsai, V. C., Stewart, A. L., and Thompson, A. F.: Marine icesheet profiles and stability under Coulomb basal conditions, J. Glaciol., 61, 205–215, 2015. a
Ultee, L., Robel, A. A., and Castruccio, S.: A stochastic parameterization of ice sheet surface mass balance for the Stochastic IceSheet and SeaLevel System Model (StISSM v1.0), Geosci. Model Dev., 17, 1041–1057, https://doi.org/10.5194/gmd1710412024, 2024. a, b
Verjans, V.: Data for “The Stochastic IceSheet and SeaLevel System Model v1.0 (StISSM v1.0)” by Verjans et al., Zenodo [data set], https://doi.org/10.5281/zenodo.7347470, 2022.
Verjans, V., Robel, A. A., Seroussi, H., Ultee, L., and Thompson, A. F.: The Stochastic IceSheet and SeaLevel System Model v1.0 (StISSM v1.0), Geosci. Model Dev., 15, 8269–8293, https://doi.org/10.5194/gmd1582692022, 2022. a, b, c, d, e
Verjans, V., Robel, A., Thompson, A. F., and Seroussi, H.: Bias correction and statistical modeling of variable oceanic forcing of Greenland outlet glaciers, J. Adv. Model. Earth Sy., 15, e2023MS003610, https://doi.org/10.1029/2023MS003610, 2023. a, b
Wood, R. R., Lehner, F., Pendergrass, A. G., and Schlunegger, S.: Changes in precipitation variability across time scales in multiple global climate model large ensembles, Environ. Res. Lett., 16, 084022, https://doi.org/10.1088/17489326/ac10dd, 2021. a
Zwally, H. J., Giovinetto, M., Beckley, M., and Saba, J.: Antarctic and Greenland Drainage Systems, https://earth.gsfc.nasa.gov/cryo/data/polaraltimetry/antarcticandgreenlanddrainagesystems (last access: 24 May 2024), 2012. a, b