Articles | Volume 15, issue 22
Geosci. Model Dev., 15, 8269–8293, 2022
Geosci. Model Dev., 15, 8269–8293, 2022
Model description paper
18 Nov 2022
Model description paper | 18 Nov 2022

The Stochastic Ice-Sheet and Sea-Level System Model v1.0 (StISSM v1.0)

The Stochastic Ice-Sheet and Sea-Level System Model v1.0 (StISSM v1.0)
Vincent Verjans1, Alexander A. Robel1, Helene Seroussi2, Lizz Ultee3, and Andrew F. Thompson4 Vincent Verjans et al.
  • 1School of Earth and Atmospheric Sciences, Georgia Institute of Technology, Atlanta, GA, USA
  • 2Thayer School of Engineering, Dartmouth College, Hanover, NH, USA
  • 3Department of Earth and Climate Sciences, Middlebury College, Middlebury, VT, USA
  • 4Environmental Science and Engineering, California Institute of Technology, Pasadena, CA 91125, USA

Correspondence: Vincent Verjans (


We introduce the first version of the Stochastic Ice-sheet and Sea-level System Model (StISSM v1.0), which adds stochastic parameterizations within a state-of-the-art large-scale ice sheet model. In StISSM v1.0, stochastic parameterizations target climatic fields with internal variability, as well as glaciological processes exhibiting variability that cannot be resolved at the spatiotemporal resolution of ice sheet models: calving and subglacial hydrology. Because both climate and unresolved glaciological processes include internal variability, stochastic parameterizations allow StISSM v1.0 to account for the impacts of their high-frequency variability on ice dynamics and on the long-term evolution of modeled glaciers and ice sheets. StISSM v1.0 additionally includes statistical models to represent surface mass balance and oceanic forcing as autoregressive processes. Such models, once appropriately calibrated, allow users to sample irreducible uncertainty in climate prediction without the need for computationally expensive ensembles from climate models. When combined together, these novel features of StISSM v1.0 enable quantification of irreducible uncertainty in ice sheet model simulations and of ice sheet sensitivity to noisy forcings. We detail the implementation strategy of StISSM v1.0, evaluate its capabilities in idealized model experiments, demonstrate its applicability at the scale of a Greenland ice sheet simulation, and highlight priorities for future developments. Results from our test experiments demonstrate the complexity of ice sheet response to variability, such as asymmetric and/or non-zero mean responses to symmetric, zero-mean imposed variability. They also show differing levels of projection uncertainty for stochastic variability in different processes. These features are in line with results from stochastic experiments in climate and ocean models, as well as with the theoretical expected behavior of noise-forced non-linear systems.

1 Introduction

Process-based numerical ice sheet models (ISMs) are the principal tool for projections of future mass balance of the Greenland and Antarctic ice sheets and their future contribution to sea-level rise. They simulate gravity-driven ice flow, given some climatic forcing and boundary conditions at the surface, basal, and lateral boundaries. In the past decade, a number of physically based ISMs aimed at simulating large-scale ice sheet dynamics and projecting future sea-level rise have been developed and/or substantially improved (e.g., Larour et al.2012; Cornford et al.2013; Gagliardini et al.2013; Pattyn2017; Hoffman et al.2018; Lipscomb et al.2019; Berends et al.2021). Recently, important advances have been made in including key physical processes (Bondzio et al.2016), data assimilation (Goldberg and Sergienko2011), adaptive grid refinement (dos Santos et al.2019), and coupling with external climatic models (Gladstone et al.2021), among other developments. In this context, there has been a growing interest in performing model intercomparison experiments, in which different models simulate ice sheet evolution over a given set of climatic projections. Such studies capture a range of possible future behaviors of the Greenland and Antarctic ice sheets, depending on the level of global warming, and quantify the uncertainty associated with discrepancies between ISMs. These recent efforts have culminated in the Ice Sheet Model Intercomparison Project for the Coupled Model Intercomparison Project 6 (ISMIP6) (Nowicki et al.2020; Goelzer et al.2020b; Seroussi et al.2020). ISMIP6 provides the basis for estimates of future sea-level rise from the ice sheets in the IPCC Assessment Report 6 (AR6; Masson-Delmotte et al.2021).

Different ISM responses to climatic forcing stem from a variety of ISM differences related to selection of physical processes, approximation of ice flow equations, model resolution, initial conditions and geometry, numerical methods, and parameterization of poorly constrained processes (e.g., Goelzer et al.2018; Seroussi et al.2019). The latter source is arguably the most investigated source of uncertainty in ISMs because it pertains to influential parameters of which direct observations are often non-existent at the ice sheet scale (Pollard and DeConto2012). Well-studied examples include the basal sliding coefficient and the ice viscosity factor; these parametric uncertainties in ice sliding and viscosity are the focus of numerous inverse technique applications (e.g., Morlighem et al.2010; Petra et al.2012; Pollard and DeConto2012; Perego et al.2014), as well as uncertainty quantification studies (e.g., Schlegel et al.2018; Aschwanden al.2019; Bulthuis et al.2019).

In addition to discrepancies between ISMs, another large uncertainty in ISM projections is attributable to future atmospheric and oceanic forcings (Nowicki and Seroussi2018; Pattyn et al.2018). Realizations of possible forcings are provided by global climate models (GCMs), which can be further downscaled with regional climate models (RCMs). However, running GCMs or RCMs requires substantial computational resources, prior to running the ISM. This limits the number of available climatic forcing scenarios that resolve the features required for ISM forcing in model ensemble studies. For example, for Antarctic simulations under a high-emission scenario, ISMIP6 used six climatic forcing scenarios, each generated from a different GCM. However, small round-off level differences in GCM initial conditions can, by themselves, yield a large spread in projected climatic fields (Kay et al.2015; Maher et al.2019). This is often referred to as internal climate variability, aleatoric climatic uncertainty, or irreducible climatic uncertainty. As a consequence, for any choice, configuration, and parameterization of GCMs and ISMs, internal climate variability is an additional and unavoidable source of uncertainty. Internal climate variability is not only a model feature, but it has also been extensively validated through observational climatic records (Mitchell1976; Chylek et al.2012; McKinnon et al.2017). It has been demonstrated that the response of ice sheet dynamics to climatic forcing is complex and non-linear (Huybrechts et al.2011; Goelzer et al.2013; Fyke et al.2018; Robel et al.2019), highlighting the importance of understanding and quantifying ice sheet sensitivity to internal climate variability (Christian et al.2022).

Furthermore, some glaciological processes exhibit variability on small spatiotemporal scales and thus cannot be resolved in current ISMs. Examples of such processes include iceberg calving, ice fracturing, and hydrology (Bassis2011; Hewitt2013; Albrecht and Levermann2014; Kingslake2015). Parameterizing such small-scale processes as deterministic forcings on the ice sheet dynamic evolution involves simplifying assumptions, which often neglect internal variability in these processes. A more accurate representation of these processes requires the inclusion of variability, rather than a constant forcing or parameter value. Stochastic parameterization has been successfully used in ocean models and climate models for the last several decades to address similar limitations, in subgrid-scale mixing processes for example (Hasselmann1976; Farrell and Ioannou1995; Porta Mana and Zanna2014; Berner et al.2017).

Studies with simple idealized models have demonstrated that mountain glaciers and marine-terminating glaciers are sensitive to internal variability within the glacier and the climate systems (Oerlemans2000; Roe and Baker2014; Robel et al.2014, 2018; Christian et al.2020, 2022). The magnitude of variability not only affects the magnitude of fluctuations in glacier length and volume, but also glacier mean state, known as noise-induced drift (Hindmarsh and Le Meur2001; Mikkelsen et al.2018; Robel et al.2018). Glaciers can also exhibit long-term fluctuations even when forced exclusively by inter-annual climatic fluctuations, due to their long memory of past forcing (Roe and O'Neal2009; Robel et al.2018). Mantelli et al. (2016) showed that natural climate variability can push ice streams away from a state of stable behavior as well as lead to stable behaviors that do not exist in the absence of variability. Such theoretical studies have prompted interest in ISM simulations accounting for climate variability. Consequently, other studies have specifically demonstrated asymmetric and non-linear responses of West Antarctic glaciers to ocean forcing variability (Snow et al.2017; Hoffman et al.2019; Robel et al.2019). Their results show that the glacier response to internal variability is complex and modulated by ice dynamics, geometric glacier configuration, and statistics of the variability imposed.

Tsai et al. (2017, 2020) used large-scale ISM simulations to evaluate the impact of internal climate variability on Greenland- and Antarctic-wide projections. Their approach used direct forcings from 40 to 50 members from coarse GCM ensemble runs and showed that the spread in ice sheet response due to internal climatic variability amounts to a significant fraction of the mean response. Nonetheless, this approach faces some limitations. Firstly, GCMs produce coarse-resolution outputs (∼100 km grid scale, monthly steps) and cannot resolve processes with strong spatial gradients. Secondly, their outputs do not have a one-to-one correspondence to inputs required for ISMs and thus require assumptions for such conversions. These limitations can be addressed by forcing a RCM at its boundaries with the climatic fields from a GCM in order to dynamically interpolate processes at a higher resolution, a method called dynamical downscaling. However, dynamically downscaling each GCM ensemble member is computationally impractical.

In this study, we describe the first large-scale stochastic ISM, the Stochastic Ice-sheet and Sea-level System Model (StISSM v1.0). StISSM v1.0 adds stochastic capabilities to the Ice-sheet and Sea-level System Model (ISSM; Larour et al.2012). The stochastic parameterizations target two different sources of variability: first, the variability internal to the ocean–atmosphere system, but external to ice sheets, and second, inherent variability in some ice sheet processes, such as calving and subglacial hydrology. We refer to this collection of processes as “forcings with internal variability”. The goal is to provide a straightforward and efficient tool to the ice sheet modeling community for assessing irreducible uncertainty in ice sheet projections. StISSM v1.0 has several advantages over more ad hoc approaches to forcing models with variable fields: it gives the user the choice of which processes exhibit stochasticity, it does not require direct forcing from ensemble runs of GCMs, it allows for different spatiotemporal autocorrelations and correlations between processes, and it exploits parallelization to efficiently facilitate large-ensemble simulations. Statistics that determine the magnitude and the spatiotemporal dimensions of variability should be provided by users and thus constrained from theory, observations, or other model simulations (Bassis2011; Chylek et al.2012; Castruccio et al.2014; Christensen2020; Hu and Castruccio2021). The latter option shows that StISSM v1.0 does not completely eliminate the need to run large climate model ensembles to constrain climatic variability. However, statistical models can be calibrated to a small number of climate model runs and implemented in StISSM v1.0 in order to reproduce the characteristics of climate model output (Castruccio and Stein2013; Castruccio et al.2014; Bao et al.2016). We emphasize that the purpose of StISSM v1.0 is not to fully explore parametric uncertainty in ice sheet modeling, which is a separate task and the subject of much ongoing research (e.g., Schlegel et al.2018; Aschwanden al.2019; Bulthuis et al.2022). Instead, it is the first integrated computational tool that focuses on the impacts of internal variability on large-scale ISM simulations, with an emphasis on usability.

2 Methods

The new stochastic capabilities are implemented within the core of the source code of ISSM. We refer readers to Larour et al. (2012) for a general description of ISSM. Usage of stochasticity is optional: if turned off, ISSM simulations are fully deterministic.

2.1 Stochastic fields

Stochastic variability can be applied to a number of variables in ISSM, independently or with intervariable correlation. The stochastic variables implemented for v1.0 of StISSM encompass both climatic forcings and unresolved glaciological processes: surface mass balance (SMB), ocean forcing, calving, and subglacial water pressure. These variables were prioritized for the implementation of stochasticity because they are known to be subject to internal climate variability (SMB and ocean forcing), and/or they reflect the impact of unresolved small-scale processes (calving and subglacial water pressure) (Ribergaard et al.2008; Bassis2011; Hewitt2013; Fettweis et al.2020). Our modeling framework ensures that, in the future, stochastic variability can easily be implemented for other variables and parameters of ISSM. In StISSM v1.0, the model schemes that parameterize these variables selected for stochasticity are prescribed SMB, autoregressive SMB, prescribed calving rate, prescribed floating ice melt rate, depth-dependent parameterization of floating ice melt, autoregressive depth-dependent parameterization of floating ice melt, autoregressive thermal forcing for terminus melting, and subglacial water pressure. We give here a brief description of the different stochastic variables implemented in StISSM v1.0 and their parameterization.

The notion of “prescribed” means that the values of a variable are explicitly provided by the user, as opposed to calculated within ISSM. They can be prescribed as either single values or varying in space and/or time. Turning on stochasticity for such variables in StISSM v1.0 implies that Gaussian white noise is added to the prescribed values at a user-defined temporal frequency in the simulation. By definition, Gaussian white noise has zero mean and is uncorrelated in time. In StISSM v1.0, a generic variable y can be modeled with additive Gaussian white noise. If y has a prescribed temporal mean value y, its value at any model time step tk is represented as

(1) y t k = y + ϵ y t k ϵ y N 0 , σ y 2 ,

where ϵy(tk) is the noise term added to y at time step tk, and σy is the standard deviation of ϵy, which must be provided by the user. As explained in Sect. 2.2, σy is fixed in time but can be variable in space, hence allowing ϵy to vary across the model domain.

A common depth-dependent parameterization of floating ice melt rate uses a piecewise linear function in depth (e.g., Favier et al.2014; Seroussi et al.2014):

(2) m fl ( z ) = m fl , up ,  if  z z up m fl , up + z - z up z dp - z up m fl , dp - m fl , up ,  if  z up > z > z dp m fl , dp ,  if  z z dp ,

where mfl(z) is the melt rate (meters ice equivalent per year, m ice eq. yr−1) at a given vertical level z (m), constrained by mfl,up at zup and mfl,dp at zdp. In this parameterization, noise is added to the mean value of the floating ice melt rate at depth, mfl,dp. As such, mfl,dp follows Eq. (1), where mfl,dp is substituted for y. The rationale for stochastic variability in mfl,dp is, for example, to represent variability in deep-water temperature below Antarctic ice shelves, which is known to have a strong impact on ice dynamics (Jenkins et al.2018).

The deterministic subglacial water pressure, pw, can be calculated in different ways in ISSM. First, it can be forced to 0 over the entire domain (pw = 0), assuming the absence of subglacial water. Second, it can be computed as

(3) p w = ρ w g ( z sl - z b ) ,

where ρw is the density of water (kg m−3), g is gravity (m s−2), zsl is the sea-level elevation (m), and zb is the elevation of the base of the ice column (m). Third, it can be computed as in Eq. (3) but set to 0 if negative. In each of these three cases, applying stochasticity adds Gaussian noise to the mean value pw, which corresponds to using Eq. (1) with pw substituted for y. The water pressure is subsequently used in the model to calculate the effective pressure Neff:

(4) N eff = ρ i g H - p w ,

where ρi is the density of solid ice, and H is the thickness of the ice column. Neff is used in various parameterizations of the basal sliding speed available in ISSM (Brondex et al.2019), allowing ϵpwtk to influence ice dynamics.

Thermal forcing, TF (K), quantifies the excess ocean temperature with respect to the freezing point of water at the interface between the front of grounded outlet glaciers and the ocean. TF thus enters in the parameterization of melt rates at the terminus of outlet glaciers, mtrm. In ISSM, mtrm follows the formulation of Rignot et al. (2016):

(5) m trm = ( A h w q sg α + B ) TF β ,

where qsg is the subglacial water flux (m d−1); h is water depth (m); and A, B, α, and β are calibration parameters. Equation (5) applies to the front of grounded outlet glaciers rather than under ice shelves and is thus more representative of Greenland than Antarctic conditions. Variability in ocean temperatures around Greenland has a large impact on outlet glacier dynamics (e.g., Straneo and Heimbach2013; Wood et al.2021), motivating our choice to prioritize stochastic variability in TF. Thus, in StISSM v1.0, TF can be modeled as an autoregressive process, as detailed in Sect. 2.3.

2.2 Numerics and spatiotemporality of stochasticity

Allowing all the ϵy(tk) terms introduced in Sect. 2.1 to vary in time and space implies that the spatial dimensions of stochasticity and the stochastic time step should be specified; these concepts are explained in this section.

The stochastic time step corresponds to the temporal frequency at which new noise terms for the stochastic variables are computed. The only restriction on the choice of the stochastic time step is that it cannot be smaller than the main time step of the numerical model simulation, i.e., the time step used by ISSM. It is important to specify the stochastic time step separately from the simulation time step because the variability in a time series depends not only on the amplitude of the noise imposed but also on the temporal frequency at which the noise is imposed. As such, if the stochastic time step was simply set equal to the simulation time step, changing the simulation time step would modify the variability imposed by the forcings to the ice sheet. When the stochastic time step is larger than the numerical model time step, then the noise term is not changed on every model time step. Integrating the noise in this fashion, and requiring the provided noise parameters to be self-consistent with the stochastic time step, means that the ice sheet responds to a forcing with characteristics of variability unaltered by other numerical considerations. Thus, the noisy forcing frequency and amplitude remain independent of the numerical model time stepping scheme, and the latter does not influence the sensitivity of the ice sheet to stochastic variability. At this stage, StISSM v1.0 uses an identical stochastic time step for all variables modeled with additive Gaussian white noise (Eq. 1), but implementing different stochastic time steps is a possible avenue for future development.

The spatial dimensions of stochasticity account for the number of sub-domains of the computational domain that share the same noise terms. The stochastic fluctuations are uniformly applied in each separate sub-domain. For example, a domain could be separated into individual glacier catchments. The number of sub-domains is prescribed during the parameterization of the model and can be as large as the number of mesh elements in the domain.

StISSM v1.0 computes all noise terms according to a Gaussian distribution. The Gaussian noise can have different correlation features in space, in time, and between variables. While future work will focus on better constraining statistical distributions of variability in the processes of interest, many geophysical processes fluctuate with a Gaussian distribution when integrated over time (e.g., Hasselmann1976). The distribution of the noise is multivariate if stochasticity is applied to several variables or if a variable has a spatial dimensionality greater than one:

(6) ϵ t k N ( 0 , Σ ) ,

where ϵt is a vector of size of the total stochastic dimensionality, dtotst, and Σ is the global covariance matrix of size dtotst×dtotst; dtotst is the sum of the number of stochastic variables times their respective spatial dimension. In this way, covariance entries can be prescribed between all the sub-domains and between any variable, thus including the different climatic forcing variables and the glaciological processes. The covariance matrix prescribed must be a valid covariance matrix, i.e., positive semi-definite with the variances σ2 along the diagonal. It is the responsibility of the user to prescribe covariance entries that suit their interest, which can be challenging, depending on the specifics of the problem investigated (e.g., Hu and Castruccio2021). Examples of relevant covariances include between SMB and subglacial water pressure, between floating ice melt rates and calving, between TF in different sub-domains, between SMB and TF, etc. The simplest choice is to set all non-diagonal entries to 0, which corresponds to statistical independence between all variables and all sub-domains.

From the stochastic time step, each simulation time step is determined as being a stochastic model step or not. At a stochastic model step, new ϵ(tk) terms are generated from Eq. (6). Otherwise, ϵ(tk) values from the previous step are re-used. At a stochastic model step, all ϵ(tk) terms are computed simultaneously and before the solvers of the governing equations of ISSM. Because ϵ(tk) terms are specific to sub-domains, an ϵ(tk) term computed for a given sub-domain is subsequently assigned to the nodal points of all the elements belonging to the sub-domain.

Stochasticity has been implemented in the most general way possible, such that developing stochasticity for a new variable would only require reproducing the code from another variable, with minimal adaptations needed for variable names and for potential specificities of the new variable. Moreover, the stochastic noise generation is mostly implemented in a separate module, thus causing minimal interference to developments of any other aspect of ISSM. All the stochastic schemes are implemented in the C++ source code of ISSM and are integral parts of the core of the model, but the schemes are not called if stochasticity is not required by the user. The random number generator implemented in ISSM is the commonly used linear congruential generator, which is a recursive algorithm with the advantages of being fast and easy to implement (Knuth1998). For the sake of reproducibility of results or troubleshooting, a randomness flag can optionally be set to false during the configuration of a simulation of StISSM v1.0. In this case, each random number generation uses a seed that is a deterministic function of the time step.

2.3 Autoregressive schemes

In an autoregressive (AR) process, the variable of interest evolves in time and depends linearly on its own previous values. AR models are a powerful tool in climatic time series analysis because they are discretized versions of differential equations. They capture characteristic timescales of geophysical processes, and they have been shown to characterize many complex climatic variables (Hasselmann1988; von Storch and Zwiers1999; Mudelsee2014). Long-term records of accumulation and atmospheric and oceanic temperatures in Greenland and Antarctica are commonly represented by AR processes (e.g., Roe and Steig2004; Thomas et al.2009; Mikkelsen et al.2018; Rosier et al.2021). For all these reasons, we have implemented AR capabilities in StISSM v1.0, but more complex time series models can be implemented in future developments. We use the notation η to represent a variable governed by a stochastic AR process, in contrast to y, which is governed by an additive Gaussian white noise process (see Eq. 1). A general AR model of order p for an autocorrelated variable η at a time step tk is then given as

(7) η t k = μ t k + i = 1 p φ i η t k - i + ϵ η t k η t k = η t k - μ t k μ t k = β 0 + β 1 t ,

where ϵη is a Gaussian noise term uncorrelated in time, and φi represents lag coefficients. The deterministic term μ(tk) allows a non-zero intercept term, β0, to be prescribed to η, as well as a linear trend, β1. The deterministic μ is removed from η to yield the terms ηtk-i. The order p of an AR process allows multiple degrees of freedom influencing η to be captured, and that may act on different timescales (von Storch and Zwiers1999). Using a higher-order AR model thus allows more complicated temporal variability to be captured but implies the risk of overfitting if the calibration time series are too short.

We have implemented AR options in StISSM v1.0 for the three climate-related variables mentioned in Sect. 2.1: (i) SMB; (ii) deep-water melt rate, mfl,dp, when using the depth-dependent parameterization Eq. (2); and (iii) thermal forcing for terminus melting, TF (Eq. 5). Thus, any of these variables can be computed with an AR model, following Eq. (7). AR capabilities can be turned on for a single variable or for multiple variables. The order p and all the coefficients β0, β1, and φi are prescribed by the user and are specific to the variable chosen. The latter parameters are all fixed in time but can vary in space as detailed below.

All the variables computed via an AR model have their specific spatial dimensions and temporal setting. The spatial dimensions work in the same way as for the generic stochasticity (Sect. 2.2). The coefficients β0, β1, and φi are thus specific to each sub-domain, as are the ϵη terms computed in Eq. (7). The sub-domains of the AR variables need to be specified and do not need to be identical to the sub-domains of the other stochastic variables. Noise terms for Eq. (7) are generated for each sub-domain of the AR variables. The global covariance matrix must include covariance terms between AR sub-domains and with the sub-domains of other stochastic variables. The AR time setting corresponds to t in Eq. (7). A variable η(tk) is recomputed following Eq. (7) at the frequency determined by tk and remains constant in between if the ISSM simulation time step is shorter than successive tk steps. The AR frequency cannot be higher than the stochastic frequency such that a new ϵη term is available every time η(tk) is recomputed via Eq. (7). Throughout a simulation, StISSM v1.0 retains in memory only those terms needed to compute AR variables at the present simulation time step, the number of which depends on the order p in Eq. (7). It should be noted that η (Eq. 7) and y (Eq. 1) are both the realization of a random process used as an ISSM variable, with the former being an autoregressive stochastic process and the latter an additive Gaussian white noise process.

The SMB AR scheme optionally allows for dynamic SMB–elevation feedback through prescribed altitudinal gradients of SMB. Such gradients, called SMB lapse rates, relate SMB values at individual mesh elements to varying elevation in space and/or time and are regularly applied for Greenland ice sheet simulations (Edwards et al.2014; Goelzer et al.2020a). Each sub-domain of the SMB AR scheme has its specific set of lapse rates and their corresponding elevation ranges in which they apply. There is no limit on the number of different lapse rates per sub-domain, as long as each lapse rate has a corresponding elevation range. The lapse rates serve to adjust SMB at individual mesh elements as a function of their elevation, with respect to the sub-domain reference SMB value calculated via Eq. (7). Furthermore, they allow for dynamic SMB feedbacks in time as ice thickens or thins throughout a simulation.

2.4 Ensemble simulations

In order to sample the component of irreducible uncertainty due to internal and climate variability, StISSM v1.0 allows for ensemble runs of the ice sheet model with stochastic parameterizations activated. Each ensemble is characterized by a selection of stochastic variables and a given configuration of the stochasticity (Sect. 2.1 and 2.2). The number of simulations, referred to as ensemble members, is chosen by the user; each member is then characterized by a unique stochastic realization. All the simulations for the different members can be run in parallel, allowing for efficient simulations. The runs of the different members are executed on different nodes, and each separate member run can further be parallelized on different processors using the usual ISSM parallelization capabilities (Larour et al.2012). StISSM v1.0 allows any of the possible output variables of ISSM to be saved and at a user-specified frequency in order to manage output size. Output variables can be scalar (e.g., total ice mass) or multi-dimensional (e.g., ice thickness at each mesh element) fields. From the outputs, ensemble statistics can be computed (example codes are provided; see the “Code and data availability” section). Log files are automatically generated for debugging purposes, as is usual for ISSM runs. The implementation of ensemble simulations in StISSM v1.0 is straightforward, as illustrated in Algorithm 1. The procedure for launching simulations does depend on the system and scheduler being used. StISSM v1.0 includes examples of ensemble launchers for common task schedulers in computer clusters.

Algorithm 1Parallelization procedure for ensemble runs.

choose stochastic variables, output variables, output frequency
configure stochasticity
n number of ensemble members
for ii=0 to n do
nameii name of simulation ii
launch simulation of nameii
end for
(n simulations run in parallel)
for ii=0 to n do
retrieve results of nameii
end for
3 Model experiments

We perform three sets of experiments to test and demonstrate the new capabilities of StISSM v1.0. The first set simulates a marine-terminating glacier with geometry taken from the benchmark configuration of the Marine Ice Sheet Model Intercomparison Project (MISMIP+, as described in Asay-Davis et al.2016). The second set simulates a quarter of an idealized circular ice sheet (IQIS) with a fast-flowing ice stream. The third simulates the Greenland ice sheet (GrIS). We detail the configuration of the three demonstration experiments in this section.

The MISMIP+ configuration is a well-known and thoroughly tested configuration. While this is the first study applying stochasticity to MISMIP+, our results can be compared to prior studies with a large range of different ISMs (Asay-Davis et al.2016; Seroussi and Morlighem2018; Cornford et al.2020). The IQIS design resembles that of a real ice sheet, but with an idealized setup. This configuration allows us to investigate the role of different types of stochastic forcings, without complications from more realistic setups. Finally, the GrIS simulations demonstrate the applicability of StISSM v1.0 to a realistic model configuration on an entire ice sheet, which is a necessary consideration for our ice sheet modeling objectives.


Our MISMIP+ experiment follows the description and parameterizations of the MISMIP+ Ice1r design in Asay-Davis et al. (2016) and Experiment 2 of Seroussi and Morlighem (2018). The domain spans 640 km in the x direction and 80 km in the y direction, with a fixed ice front at x=640 km (Fig. 1). Along the centerline y=40 km, the bed topography is 150 m below sea level at x=0 km and progressively decreases in the x direction, but with a bed bump between x=390 and 506 km (Fig. 1). The spin-up run starts from a thin (10 to 100 m thickness) glacier over the entire domain, is fully deterministic, has a horizontal resolution of 1 km, and uses the shallow-shelf approximation (MacAyeal1989). The glacier progressively builds up through a constant positive SMB set to 0.3 m ice eq. yr−1 and is uniform over the domain. For the melting rates applied to floating ice, we use the depth-dependent melting parameterization Eq. (2). We fix the values mfl,up=0 m ice eq. yr−1, zup=-50 m, mfl,dp=1 m ice eq. yr−1, and zdp=-500 m. Similarly to Seroussi and Morlighem (2018), we use a Weertman-type sliding law (Weertman1957):

(8) τ b = - C W | | u b | | 1 m - 1 u b ,

where τb is the basal stress (Pa), ub is the basal velocity (m yr−1), m=3, and CW is a basal drag coefficient taken equal to 1.0×104 Pa m1/3 yr1/3. A steady state is reached after 19 000 years, at which point SMB is balanced by the melt of floating ice and ice flow out of the model domain. The time step used during the spin-up is 1/2 year over the first 15 000 years and refined to 1/4 year over the last 4000 years. We use this steady state as an initial state for transient experiments. Under such relatively low-melt conditions, the initial state is characterized by a thick ice-shelf, buttressing and stabilizing the grounded part of the glacier, and the grounding line along the center flowline is located within a bed trough at x=397.0 km (Fig. 1). The relative changes in ice mass and in grounded ice area over the last 1000 years of the spin-up are both <0.15 %.

Figure 1Initial steady-state configuration for the MISMIP+ experiments. (a) Ice thickness in (x,y) plane, (b) vertical profile along y=40 km.


From this steady state, we perform five ensembles of transient simulations. The transient experiments are performed over a period of 500 years, with a time step of 1/4 years and with the same spatial resolution of 1 km as in the spin-up run. Each ensemble consists of 500 individual simulations, referred to as members. Each ensemble member has the mean melt rate signal of mfl,dp=1 m ice eq. yr−1 but applies stochastic perturbations ϵmfl,dp (Eqs. 1, 2, 7). While mfl,dp is equal to the forcing melt rate in the spin-up, the stochastic perturbations will cause deviations of the glacier from its initial deterministic steady state. Ensembles differ in the statistics of the stochastic perturbations, while members of the same ensemble only differ by their unique realization of random noise. Through our five ensembles, we evaluate the sensitivity of the model configuration to persistence of ocean melt variability in time. All the ensembles represent mfl,dp as an order 1 (i.e., p=1 in Eq. 7) AR process, commonly written AR(1). We apply no trend in mfl,dp (i.e., β1=0 in Eq. 7), and we set the constant term of the AR(1) model to the mean melt rate (i.e., β0=mfl,dp=1 m ice eq. yr−1 in Eq. 7) and use a stochastic time step of 1 year, i.e., annual variability in the stochastic perturbations. In an AR(1) model, the autocorrelation coefficient ρ1 is equivalent to φ1 in Eq. (7) and can be converted to a more intuitive characteristic timescale, τ, via (e.g., Burke and Roe 2014)

(9) τ = Δ t 1 - ρ 1 ,

where Δt is the time discretization of the stochastic time series, here set to 1 year. We test a range of five values of τ: 1, 2.5, 10, 20, and 50 years. Note that τ=1 year corresponds to annual uncorrelated white noise (equivalent to Eq. 1). We denote these five ensembles as MISMIPτ=1, MISMIPτ=2.5, MISMIPτ=10, MISMIPτ=20, and MISMIPτ=50, respectively. For all ensembles, the standard deviation of the Gaussian noise term is set to 1/3 of the deterministic mean mfl,dp, σmfl,dp = 1/3 m ice eq. yr−1:

(10) ϵ m fl , dp N 0 , σ m fl , dp 2 for  MISMIP τ = 1 , 2.5 , 10 , 20 , 50 .

In all the MISMIP+ transient experiments, the ϵmfl,dp perturbations are uniformly applied over the entire model domain.

3.2 Idealized Quarter Ice Sheet

The Idealized Quarter Ice Sheet (IQIS) experiment is performed on a square domain of 750 km×750 km and meant to represent a quarter of a circular marine ice sheet in an idealized configuration (Fig. 2). We impose Dirichlet conditions with values of 0 m yr−1 for the x-direction velocities at y=0 m and y-direction velocities at x=0 m to represent ice divides. By symmetry, the four quarters of the fully circular ice sheet would each be identical; thus we only simulate one quarter to save computational expense. The bed topography decreases linearly along the radius, from 10 m above sea level at the origin (lower-left corner) to 20 m below sea level at x=750 km, y=750 km (upper-right corner), corresponding to a very gradual bed slope in the radial direction of 4×10-5. The friction law used for parameterizing basal sliding is the Budd friction law (Budd et al.1984):

(11) τ b = - C B u b N eff ,

where CB is a basal drag coefficient set to 1600 m−1 yr. We lower CB to 400 m−1 yr on a 57 km wide strip starting 200 km away from the origin, simulating the presence of a fast-flowing ice stream (Fig. 2).

Figure 2Initial steady-state configuration for the IQIS experiments. (a) Ice thickness, (b) ice velocities.


The deterministic steady state of the ice sheet is constructed via a balance between constant positive SMB over the domain (0.4 m ice eq. yr−1) and frontal ablation at the grounded ice–ocean front. We split the frontal ablation between a calving rate flux, crate, and a frontal melt rate, mtrm, following Eq. (5) with qsg=0 m d−1 and TF=4 K. Note that crate and mtrm have the same effect: ice ablation at the terminus. The purpose of using these two terms is to allow for separate stochastic perturbations in calving and TF during the transient experiments. We set crate quadratically increasing along the radius, and this radial dependence is enforced by a multiplicative factor, cratefac, constant across the domain:

(12) c rate = c rate fac rad 750 × 10 3 2 ,

where “rad” denotes the radius (m), and cratefac is set to 690 m ice eq. yr−1. This radial dependence allows the front to reach a stable position, where frontal ablation balances total SMB. The spin-up simulation is run for 79 000 years. In the final stage of the spin-up to a steady state (last 10 000 years), we use a horizontal resolution varying between 15 km in the interior and 800 m at the front, a time step of 1/4 year, and the mono-layer higher-order stress approximation. This stress approximation effectively captures the same longitudinal extension and vertically shearing flow as higher-order flow approximations, but without needing vertical model layers (dos Santos et al.2022). The relative changes in total mass and area of the ice sheet over the last 1000 spin-up years are both <4×10-3 %. The final steady-state ice thicknesses and velocities of IQIS are shown in Fig. 2. It has an ice-covered area of ∼31 000 km2 and an ice mass of 350 715 Gt, equivalent to 1/5 of the area and 1/8 of the mass of the Greenland ice sheet.

We use this steady state as an initial state for transient ensemble experiments with stochasticity. We perform four sets of ensemble experiments (Table 1). Each of the four sets is characterized by activating stochastic fluctuations for specific forcings, while the mean forcings remain the same as in the spin-up. Thus, SMB = 0.4 m ice eq. yr−1, cratefac=690 m ice eq. yr−1 (Eq. 12), qsg=0 m d−1, and TF=4 K (Eq. 5). As such, our experiments quantify how fluctuations in different processes around a given mean forcing can drive the IQIS away from its steady state.

Table 1Configuration of the IQIS transient ensemble experiments. τ=1 year corresponds to annual noise uncorrelated in time. σcrate is equal to 1/3 of the mean of Eq. (12) evaluated at the front location. Note that all σ values are taken as 1/3 of the corresponding mean forcing. The noise correlation between variables (ϵSMB and ϵcrate in the exterior sub-domain in IQIS3) is represented by rvariables, and rdom is the correlation in ϵSMB between the interior (rad≤300 km) and exterior (rad>300 km) parts of the domain.

Download Print Version | Download XLSX

The configuration of the IQIS transient experiments is detailed in Table 1. The experiments are designed to yield meaningful comparisons between ensembles with variability in different variables. IQIS1 and IQIS2 show the relative strength of SMB versus calving fluctuations on our experimental design. IQIS3 shows the impact of combining the SMB and calving stochastic perturbations. Finally, IQIS4 shows how decadal persistence of oceanic forcing at the ice front can impact ice sheet dynamics.

The ensemble experiments consist of 500 members and a simulation period of 500 years, and the spatial and temporal resolutions are kept identical to the final spin-up configuration. The stochastic time step is set to 1 year, such that the fluctuations imposed have an annual sampling frequency. The first set, IQIS1, applies annual white noise in SMB (see Eq. 1). The standard deviation of the noise amplitude, σSMB, is taken as 1/3 of the mean (Table 1), comparable to the SMB relative inter-annual variability in Greenland (Fettweis et al.2020). Here, we separate the domain into interior (rad≤300 km) and exterior (rad>300 km) sub-domains. The correlation in ϵSMB between both parts is set to 0.6. The second set, IQIS2, applies annual white noise in calving rates, ϵcrate (see Eq. 1). We choose the same 1/3 ratio of noise-to-signal as for SMB in IQIS1 for better comparability. Because calving rates vary in space (Eq. 12), we take σcrate as 1/3 of the mean calving rates at the ice front at the end of the spin-up phase (Table 1). The third set, IQIS3, combines the stochastic forcings of IQIS1 and IQIS2 in a single ensemble experiment (Table 1). It uses a negative correlation between calving and SMB in the exterior sub-domain to represent, for example, the impact of increased surface melt on hydrofracturing (Benn et al.2007). The last set of IQIS experiments, IQIS4, models TF (Eq. 5) as an AR(1) process (Eq. 7) without a trend and with the same mean as in the spin-up (4 K) (i.e., β0=4 and β1=0 in Eq. 7). We set σTF to 1/3 of the mean signal for ease of comparisons with the other experiments, but here we apply a decadal noise persistence (Table 1). In Eq. (5), we use the parametrization of Rignot et al. (2016) and set A=3×10-4, B=0.15, α=0.39, and β=1.18.

3.3 Greenland ice sheet

To demonstrate that StISSM v1.0 is readily applicable at the scale of ice sheet simulations, we simulate the evolution of the GrIS with stochastic SMB and ocean forcings. The configuration uses an initial state matched to observations but is spun up to reach a deterministic steady state before launching the transient experiments with stochasticity applied. The initial state uses the bed topography, the ice thickness and the ice mask from BedMachine v4 (Morlighem et al.2017), the ice velocity field from Joughin et al. (2017), the geothermal heat flux from Shapiro and Ritzwoller (2004), and surface temperatures from Ettema et al. (2009). We solve for a thermal steady state in three dimensions, with 10 vertical layers, to compute vertical temperature profiles and then calculate the ice rheology field from the depth-averaged temperatures following Cuffey and Paterson (2010). We invert for basal friction coefficients CB in the Budd sliding law (Eq. 11). We use a linear regression of CB on bed elevation to extrapolate the field in areas covered by less than 500 m of ice. This avoids spurious patterns in CB in marginal areas where observed velocity gradients are large and also allows extrapolation of CB in ice-free areas where the ice sheet can extend during model simulations.

After this initialization, we perform a deterministic spin-up in order to reach a GrIS configuration in a steady state. We emphasize that the purpose of our simulations is not to predict future ice mass balance of the GrIS, and we do not argue that the real GrIS is in a steady state; our goal with this spin-up is to have a steady baseline against which to compare a transient ensemble. We separate the GrIS in 19 different basins following an existing delineation (Zwally et al.2012). The computation of SMB is basin-specific and calibrated to the mean 1961–1990 modeled SMB field of RACMO2 (van Angelen et al.2014). We fit piecewise linear functions of elevation with two breakpoints to the SMB data in order to derive three SMB lapse rates per basin, although only two lapse rates (i.e., a single breakpoint) are used if the fitting yields an unrealistic lapse rate over a narrow elevation range (Table 2). The reference SMB in each basin corresponds to the basin-specific piecewise linear function evaluated at the mean basin elevation (Table 2). Stochasticity is turned off during the spin-up, but SMB values of any mesh element do change in time because they are adjusted as thickening and thinning patterns affect the elevation of the elements.

Table 2Climatic forcing in each basin for the GrIS simulations. Basin numbers correspond to the delineation in Zwally et al. (2012), shown in Fig. (3). Piecewise linear functions of SMB to elevation are fitted to the RACMO2 SMB in each basin. The fit uses two breakpoints to derive three lapse rates. If the fit results in one of the lapse rate values being unrealistic and applying over a narrow elevation range, we use a fit with a single breakpoint and two lapse rates. The reference SMB corresponds to the piecewise linear function evaluated at the mean elevation. Ocean TF is the TF value used in Eq. (5) to model frontal melt during the second phase of the spin-up at outlet glaciers where we parameterize ocean melt.

Download Print Version | Download XLSX

The spin-up itself is separated in two different phases, both of them using a weekly time step and the two-dimensional shallow-shelf approximation. In the first phase, we fix the ice sheet margin positions and implement free-flux boundary conditions at the ice margins, meaning that boundary ice fluxes adjust to incoming fluxes to keep margins fixed in space. We use a spatial resolution ranging from 25 km in the slowest-flowing areas to 2 km in the fastest-flowing areas. During this first spin-up phase, the modeled ice sheet adjusts to the SMB field until it reaches a steady state. The steady state of this first phase requires a dynamic equilibrium between ice flow and the SMB field and thus takes 30 000 years.

In the second phase of the spin-up, we allow for moving margins at 11 of the major outlet glaciers of the GrIS, where we parameterize ocean melt (Fig. 3). Resolving these 11 outlet glaciers and migration of their termini requires finer meshing close to their termini. Frontal ablation at the ice–ocean boundaries is applied through a combination of melting at the termini (Eq. 5) and calving. Values of TF are specific to each basin and taken as the approximate 1990–2018 average values reported by Wood et al. (2021) (Table 2). Calving values are prescribed as constant for each of the 11 glaciers and taken to minimize the departure from the steady state reached in the first phase (Table 3). We refine the mesh resolution at the 11 outlet glaciers down to 800 m while keeping the fixed front and free-flux boundary conditions for the other outlet glaciers. We select only 11 moving glaciers to avoid the computational load associated with refined horizontal resolution at all outlet glaciers. Our simulation set-up represents the impact of ocean forcing on the GrIS reasonably well since ice discharge is dominated by a small (<20) number of glaciers (Enderlin et al.2014).

Table 3Calving rates applied at the 11 outlet glaciers where terminus migration is simulated and ocean melt parameterized. Each Roman numeral is associated with a corresponding glacier, shown in Fig. 3. Total frontal ablation results from both calving and melt at the terminus, which is governed by basin-specific TF values given in Table 2.

Download Print Version | Download XLSX

Comparing the simulated GrIS state at the end of the spin-up to observations, the total ice mass and ice-covered area are ∼11 800 Gt (0.4 %) and ∼71 000 km2 (3.5 %) lower, respectively. The main differences are an interior thickening in the south and an inland retreat in the north (see Appendix A). Velocity magnitudes and patterns remain consistent with observed velocities (Fig. 3). The five simulated outlet glaciers in the north are retreated compared to observations, except for glacier II (Petermann). Over the last 200 years of spin-up, the changes in total ice mass and ice-covered area are +24.9 Gt (+0.001 %) and −316 km2 (−0.016 %). More details concerning the GrIS steady state are given in Appendix A.

Figure 3Steady-state ice velocities at the end of the second phase of the GrIS spin-up. Arabic numerals (black) show the individual basins. Roman numerals (cyan) show the outlet glaciers where ice front movement is simulated and ocean melt parameterized.

We use the GrIS final steady state as an initial state for our transient experiments with stochasticity turned on, which will cause deviations from the steady state. We perform a single transient ensemble of 200 members over 500 years, with a stochastic time step set to 1 year, thus representing annual fluctuations. While this initial state is very close to a steady state, we still perform a deterministic control run of 500 years to quantify the amount of deterministic model drift, which is minimal. In the stochastic transient ensemble, we apply stochastic fluctuations in the climate forcing fields SMB and TF. We represent both of these forcings as AR(1) processes (Eq. 7). The means of the forcings (i.e., β0 in Eq. 7) are kept equal to the constant values applied during the spin-up (Table 2), and we do not impose any trend (i.e., β1=0 in Eq. 7). We set τ=3 years for SMB (following Mikkelsen et al.2018) and τ=10 years for TF. Each basin represents an individual spatial dimension for the covariance matrix (Fig. 3). Since our simulation set-up uses 19 basins and 2 variables with stochastic variability (SMB and TF), the covariance matrix is of dimensions 38×38. The standard deviations of the noise terms are all set to 1/3 of the mean forcing:

(13) σ SMB , b = 1 3 SMB ref , b for  b = 1 , , 19 σ TF , b = 1 3 TF b for  b = 1 , , 19 ,

where the basin-specific reference SMB, SMBref,b, and TFb forcings are given in Table 2. To prescribe all the correlations involving ϵSMB,b and ϵTF,b, and between all basins, we use simple and reasonably realistic assumptions: (i) a variable is more strongly correlated with values in the neighboring basins than in the non-neighboring basins; (ii) ϵSMB and ϵTF are negatively correlated; and (iii) between different basins, ϵSMB or ϵTF variables are more strongly correlated with themselves than with the other variable. Under the constraint of using a positive semi-definite covariance matrix, all correlation absolute values range between 0.4 and 0.6 (see Appendix B for details). In the transient simulations, SMB lapse rates (Table 2), calving rates (Table 3), the weekly model time step, and the spatial resolution are kept identical to the second phase of the spin-up.

4 Results

In this section, we analyze the results of the transient MISMIP+, IQIS, and GrIS experiments in terms of total ice mass (Gt) evolution. While our analyses focus on a variable summed over the entire domain, we note that localized changes in ice thickness and/or ice extent occur and may be larger than the global patterns.


At the start of the MISMIP+ transient experiments, the total mass is 39 097 Gt, and the grounding line position is at x=397 km, stable and within a bed trough (Fig. 1). Quantitative results are presented in Table 4 and Fig. 4, while the evolution of individual ensemble members and of the probability density functions (PDFs) are displayed in Fig. 5. We base our analysis on the leading-order statistical moments of the final ice mass PDF in each ensemble, focusing on the mean, standard deviation (σE), and skewness of the ensemble (Table 4, Fig. 4). We use the Shapiro–Wilk test to evaluate if the final ice mass PDFs are consistent with a normal distribution (Shapiro and Wilk1965). This test measures the fit between standard normal quantiles and the ordered and standardized ice mass values of the ensemble. It has also been shown to be more powerful than many commonly used normality tests (Razali and Wah2011). According to the Shapiro–Wilk test at the 5 % significance level, the final PDF of each ensemble is consistent with a normal distribution (Table 4). In all the ensembles, the mean final ice mass is larger than the initial mass, but by less than 1σE. The PDFs of the final glacier state have not yet converged after 500 years, as both the mean and σE of all ensembles still show statistically significant positive trends over the last 50 years of simulation (p values <0.001).

Table 4Statistics of the final ice mass distributions for the five MISMIP+ transient ensembles. The relative change is calculated with respect to the initial ice mass (39 097 Gt).

Download Print Version | Download XLSX

Figure 4Evolution throughout the transient experiments of (a) the standard deviation and (b) the skewness in total ice mass for the five MISMIP+ transient ensembles. (c) Boxplots of the final ice mass distributions. Red dots indicate ice masses beyond 1.5 times the interquartile range from the quartile.


Figure 5Left: change in ice mass throughout the transient experiments for each ensemble member in the five MISMIP+ transient ensembles. Thick lines show the ensemble means. The right y axis shows change relative to the initial ice mass. Right: corresponding PDFs of the ice mass change after 25, 50, 100, and 500 years of simulations.


The cause of mean gains in ice mass for the MISMIP+ transient ensembles relates to the initial grounding line position in a bed trough (Fig. 1) and to the form of the melt forcing imposed. Due to the depth-dependent melt parameterization (Eq. 2), the initial geometric configuration causes an asymmetry in melt rate response: both advancing and retreating glaciers move their grounding line depth to higher elevations, hence decreasing their respective mean melt rate at the grounding line. This effect limits further retreat of retreating glaciers, while it favors further advance of advancing glaciers. Moreover, the bed slope itself is asymmetric, being steeper upstream than downstream of the trough (Fig. 1). As a consequence, the mean final glaciers are more advanced and with a larger ice mass than in the initial state. These results demonstrate that noisy forcing with zero mean can drive non-zero mean response of the glacier, which is known as noise-induced drift (e.g., Hindmarsh and Le Meur2001; Mikkelsen et al.2018; Robel et al.2018). Here, the noise-induced drift can be associated with asymmetries in geometrical configurations and forcings and with non-linearities in ice dynamics. In addition to the factors mentioned above, there are other non-linearities at play that could contribute to the noise-induced drift: the floating area affected by melt varies as the grounding line migrates; lateral stresses exerted from the grounded parts on the domain sides depend on the geometry; and ice shelf buttressing depends non-linearly on ice shelf length, ice shelf thickness, and ice thickness at the grounding line (Haseloff and Sergienko2018). The individual contributions of all these different factors on the noise-induced drift are difficult to disentangle.

Furthermore, the mean gain in ice mass is highest for the ensembles with higher τ values: MISMIPτ=20 and MISMIPτ=50. The spread in final ice mass, quantified by σE, increases sublinearly with τ. For example, compared to MISMIPτ=1, MISMIPτ=10 yields a σE 8.9 times larger, while MISMIPτ=50 yields a σE 36.9 times larger. The skewness of lower characteristic timescales (τ≤10 years) shows more variability over time (Fig. 4); short time persistence of the noise in melt rates frequently causes temporary glacier excursions that deviate strongly from the ensemble mean. In contrast, skewness of the higher characteristic timescale ensembles, MISMIPτ=20 and MISMIPτ=50, is smoother. With long time persistence of the noise applied, unusual retreats or advances of the MISMIP+ glacier are slower to develop but also to recover back towards the ensemble mean. While MISMIPτ=20 and MISMIPτ=50 show the largest spread in ice mass, their skewness is close to 0 throughout the 500 years, indicating that large and small glacier outliers are approximately equally likely and of similar magnitude in terms of mass difference with respect to the ensemble mean.

Roe and Baker (2016) derived an expression relating variability in glacier length, σL(τ), to τ:

(14) σ L ( τ ) = σ L ( τ = Δ t ) 2 τ Δ t - 1 1 / 2 ,

where we use Δt=1 year, and σL(τt) corresponds to glacier length variability under white noise forcing. This expression was derived from an analytical three-stage mountain glacier model and found to match results of a numerical flowline model (Roe and Baker2014, 2016). Equation (14) predicts that the variability amplitude in glacier length increases with the square root of the characteristic timescale of the noise. In our experiments, the relationship between τ and the standard deviation of the grounding line position is sublinear but does not follow the square root dependence predicted by Eq. (14). The absence of match with the theory predicted by Roe and Baker (2016) possibly illustrates that lateral shearing cannot be neglected in the MISMIP+ configuration, that Eq. (14) does not hold on irregular bed slopes, and/or that Eq. (14) is only suited for glaciers where mass loss is controlled by SMB rather than by melt beneath buttressing floating ice. Ultimately, even idealized marine-terminating glaciers, such as the MISMIP+ configuration, include a much wider range of processes than the simple mountain glacier considered by Roe and Baker (2016).

4.2 Idealized Quarter Ice Sheet

The initial state of the IQIS is in equilibrium; thus any deviation from the initial state is attributable to the stochastic fluctuations imposed in our IQIS1, IQIS2, IQIS3, and IQIS4 transient experiments (Table 1). The initial ice mass is 350 715 Gt. Table 5 and Fig. 6 show the quantitative results of each ensemble experiment, and Fig. 7 displays the evolution of each ensemble member and of the PDFs over time. According to the Shapiro–Wilk normality test (Shapiro and Wilk1965), all the final ice mass distributions are consistent with normality at the 5 % significance level (Table 5). As in the MISMIP+ experiments, the ice mass distributions have not reached a statistical steady state after 500 years of simulations, as both the mean and the σE values of all ensembles still show statistically significant trends in the last 50 years (p values <0.001) (see also Fig. 6). For IQIS1 (i.e., stochastic SMB) and IQIS4 (i.e., stochastic TF), the mean final ice mass is −7.6 % and +5.4 % of their respective σE away from the initial equilibrium ice mass. These small deviations of the mean state show that the stochastic fluctuations applied in these two ensembles do not cause significant noise-induced drift after 500 years.

Table 5Statistics of the final ice mass distributions for the four IQIS ensembles. The relative change is calculated with respect to the initial ice mass (350 715 Gt).

Download Print Version | Download XLSX

Figure 6Evolution throughout the transient experiments of (a) the standard deviation and (b) the skewness in total ice mass for the four IQIS ensembles. (c) Boxplots of the final ice mass distributions. The dashed line shows the initial ice mass, from the deterministic steady state. Red dots indicate ice masses beyond 1.5 times the interquartile range from the quartile. In labels, variables between brackets denote variables with stochastic variability imposed (Table 1).


Figure 7Left: change in ice mass throughout the transient experiments for each ensemble member in the four IQIS ensembles. In labels, variables between brackets denote variables with stochastic variability imposed (Table 1). Thick lines show the ensemble means. The right y axis shows change relative to the initial ice mass. Right: corresponding PDFs of the ice mass change after 50, 150, 300, and 500 years of simulations.


In contrast, the results from IQIS2 (i.e., stochastic crate) demonstrate strong evidence of noise-induced drift caused by stochastic white noise fluctuations in crate. The mean final ice mass is 2035 Gt larger than the initial mass, which corresponds to >3.5σE and to >0.5 % of the total ice mass. Since crate cannot be negative at any given time step, crate is set to 0 when ϵcrate pushes it below 0. This causes a slight asymmetry in the distribution of the ϵcrate applied. However, this effect is negligible as only 33 of the 500 ensemble members have crate=0 at any time step, and all ensemble members have a final ice mass larger than the initial mass. The latter is true even for ensemble members for which the stochastic realization of their ϵcrate time series causes a total calving flux larger than the total calving flux of the deterministic scenario without stochastic perturbations. As such, it is the stochastic nature of our calving forcing and the ice sheet response to the stochastic forcing that cause the emergence of noise-induced drift, and not the total cumulative calving flux.

The perturbations ϵcrate and ϵTF, through their influence on mtrm (Eq. 5), both impose stochastic noise on the frontal ablation. One could therefore expect similar noise-induced drift to appear in IQIS4, but there is no statistically significant evidence for this in our results after 500 years of transient simulation. Separate tests (not shown) demonstrate (i) that noise-induced drift appears in IQIS4 if τ is reduced to 1 year (i.e., annual white noise) and (ii) that noise-induced drift in IQIS4 becomes statistically significant if more time is allowed for convergence of the PDF while keeping τ=10 years. Thus, the longer persistence time of the noise (τ=10 years) requires longer timescales for the PDF to statistically converge and thus for the noise-induced drift to be significant with respect to the ensemble spread.

The results of IQIS3 (i.e., stochastic SMB and crate) show a strong noise-induced drift of approximately the same magnitude as IQIS2 (Table 5). This is consistent with the same stochastic variability applied on crate in both ensembles (Table 1). However, the additional stochasticity on SMB increases the spread, and the initial steady-state ice mass is within the final interquartile range, but still >1σE lower than the final mean. The σE of IQIS3 is very close to the sum of those of IQIS1 and IQIS2 and is 98 Gt larger than if the latter were added in quadrature with the additional covariance factor accounting for the correlation between ϵSMB and ϵcrate. As such, variability in the forcings does not have a one-to-one correspondence with variability in the final ice sheet mass.

Applying decadal variability in ocean thermal forcing (IQIS4) leads to the highest spread in the final ice mass distribution. Furthermore, IQIS4 also shows the strongest skewness throughout the 500-year period, and it is consistently negative (Fig. 6). The negative skew is driven by a larger number of outliers in the lower tail compared to the higher tail (Fig. 6). As such, time persistence in the stochastic ocean forcing causes more scenarios of extreme mass loss compared to the mean response, and the evolution of the skew contrasts with that of the MISMIP+ experiments. This confirms that response to stochastic forcing is not only asymmetric but also that the asymmetry depends on the ice sheet state.

4.3 Greenland ice sheet

Results of the GrIS ensemble with correlated stochastic variability in SMB and TF forcings are shown in Fig. 8 and Table 6. We compare the variability in the ensemble to the drift in the deterministic control run. The initial ice mass is 2.743×106 Gt, and the deterministic drift over the 500 years is small but non-zero, up to +678 Gt (+0.02 %) at the end of the 500-year run (Table 6). In contrast, the standard deviation in the final ice mass of the stochastic ensemble is 8907 Gt, thus >13 times larger than the deterministic drift. The mean final ice mass is 3175 Gt (−0.1 %) lower than the initial mass, which is thus 13σE away from the final mean, and remains within the final interquartile range (Fig. 8 and Table 6).

Table 6Statistics of the ice mass distributions of the GrIS ensemble after 125, 250, 375, and 500 years of simulation. The relative change and deterministic drift are calculated with respect to the initial ice mass (2 743 269 Gt).

Download Print Version | Download XLSX

The skew in ice mass varies between positive and negative phases, while the ensemble mean is strongly decreasing and the ensemble spread strongly increasing after 500 years (Fig. 9 and Table 6). This shows that the distribution is still far from a statistical steady state, as both the ensemble mean and σE still show significant trends over the last 50 years of simulation (p values <0.001). We note that throughout the 500 years, the ensemble ice mass distribution remains consistent with normality (Table 6). Among the 200 ensemble members, the highest and lowest final ice masses are +2.45 and −3.23σE away from the mean, respectively, hence contributing to the final negative skew. This also means that the maximum range of final ice mass in our ensemble amounts to 50 700 Gt, i.e, 1.8 % of the modeled GrIS. In the same way, after 500 years, the ±2σ range of our ensemble is 1.3 % of the GrIS mass (Table 6). As illustrated by the distribution of the final ice mass, the GrIS exhibits a non-zero mean response to noisy forcings (Fig. 8), despite the imposed climatic fluctuations being symmetric around 0 and its initial state being close to equilibrium. While the initial mass is within the ±1σ range of the mean final mass, only 37 % of the ensemble members show a mass increase (Fig. 8). The persistence of an approximately linear decrease in the mean over the 500-year period suggests that the mean of the converged state would be even lower over the thousands of years likely necessary to reach a statistical steady state.

Figure 8(a) Change in ice mass throughout the transient experiment for each ensemble member in the GrIS ensemble. The thick line shows the ensemble mean. The right y axis shows change relative to the initial ice mass. (b) Boxplot of the final ice mass distribution. Red dots indicate ice masses beyond 1.5 times the interquartile range from the quartile. (c) PDF of the ice mass change after 50, 150, 300, and 500 years of simulations.


Figure 9Evolution throughout the transient experiment of (red) the standard deviation and (green) the skewness in total ice mass for the GrIS ensemble.


5 Discussion

Stochastic modeling is well established in climate modeling (e.g., Porta Mana and Zanna2014; Berner et al.2017). This modeling approach is based on a rigorous mathematical framework, originating from stochastic differential equations and statistical physics (Majda et al.2001; Franzke et al.2015). In climate models, stochastic parameterization of internal variability and unresolved processes has been shown to improve the skill of probabilistic forecasts, reduce systematic model errors, capture regime transitions, and modify the modeled response to external forcing (Berner et al.2017; Palmer2019). StISSM v1.0 represents the first attempt to include stochastic parameterizations in large-scale ISMs. Our results for simplified and idealized model experiments show that features similar to those observed in stochastic climate modeling occur in large-scale ISMs, such as noise-induced drift and a modified mean response to external forcing. These features are in agreement with previous simple model experiments (e.g., Mikkelsen et al.2018; Robel et al.2018). Our results also reveal that simple noise terms propagate to ice sheet evolution uncertainty in a complex way because of the high degree of non-linearity in ice sheet dynamics and may be more nuanced than theory (Roe and Baker2016; Robel et al.2019). As such, model simulations are required to quantify the response of any particular glacier or ice sheet configuration to climate and internal variability, and this response cannot be trivially estimated a priori.

Irreducible uncertainty is not quantified in current ice sheet model intercomparison projects (Goelzer et al.2020b; Seroussi et al.2020). Stochastic parameterizations facilitate quantification of the irreducible uncertainty component in ice sheet projections, which is an integral part of any model prediction. This component of uncertainty is expected to be larger in systems with potential dynamical instabilities, in reality more pronounced than shown in our idealized experiments (Robel et al.2019).

StISSM v1.0 allows for stochasticity in variables which exhibit internal variability. The features of spatiotemporal correlation can be prescribed, as well as intervariable correlations. Our model experiments show that the stochastic parameterizations implemented are functional and can be used at ice sheet scale. We have aimed at making StISSM v1.0 as user-friendly as possible, in such a way that any user familiar with ISSM should find the use of StISSM v1.0 straightforward. Ensemble runs and parallelization allow for adequate sampling of irreducible uncertainty in model simulations. In general, a StISSM v1.0 simulation run with stochastic parameterizations uses additional computational resources that are negligible compared to a corresponding deterministic simulation. As in methods for estimating the role of parameter uncertainty in ice sheet evolution (e.g., Schlegel et al.2018; Aschwanden al.2019; Bulthuis et al.2019), the main computational expense comes in running many simulations. The autoregressive modeling capabilities implemented offer a computationally fast way to generate climatic forcings with prescribed statistical properties, thus sampling natural climate variability without the need to run a costly climate model for each ensemble member. In ISMs, variability in some glaciological processes such as calving and supra-, intra-, and subglacial water movement are particularly difficult to simulate; accurately resolving these processes remains elusive. Stochastic parameterizations of such unresolved processes provide a way forward to better account for their impacts on large-scale and long-term ice dynamics. In this first version of a stochastic ISM, we have implemented simple forms of stochastic processes and statistical generators of climate forcing: additive Gaussian white noise and autoregressive time series models, respectively. This lays the groundwork for future, more sophisticated schemes specifically calibrated to represent the details of variability in glaciological and climatic processes. In particular, priorities are to implement seasonality in the statistical models, incorporate more complete time series models (e.g., autoregressive moving average, ARMA), and to allow for other forms of noise forcing in order to represent non-Gaussianity in components of the climate and ice sheet systems (e.g., Perron and Sura2013).

A practical question that arises concerns the number of members needed per ensemble. Here, we have used 500 members for the MISMIP+ and IQIS ensembles and 200 for the GrIS ensemble to limit computational expense. As the number of members increases, the statistics of the ensemble progressively converge to the statistics of the true underlying distribution. In other words, results from ensembles with increasingly more members converge to the results of an ensemble with infinitely many members. Convergence plots of the statistics of interest, such as the final mean, standard deviation, and skew in our case, show their progressive convergence and are a useful tool in evaluating the number of members needed. We show such an analysis of our results in Appendix C, demonstrating adequate convergence of the ensemble statistics with 100 to 150 members in our experiments. It must be kept in mind that the number of members needed to represent the true distribution depends not only on non-linearities in the system that cause larger variability between members, but also on the timescale of stochastic variability imposed and on the statistics of interest. For example, correctly estimating the 99th percentile of the distribution requires a larger ensemble than for estimating the mean, as the effective number of members influencing this statistic is smaller.

With this first version of a large-scale stochastic ISM, several future research priorities can be identified. First, the statistics of variability in the unresolved processes need to be evaluated in order for the stochastic parameterizations to capture their impacts on ice dynamics accurately. This could be achieved via theoretical, observational, and/or high-fidelity modeling studies (e.g., Bassis2011; Hewitt2013; Emetc et al.2018; Christensen2020; Åström et al.2021; Hu and Castruccio2021). We anticipate providing workflows for estimating these parameters in future studies. Such workflows will provide calibrated covariance matrices and autoregressive parameters ready to use for StISSM v1.0 simulations. Better representing the statistics of stochastic forcings may necessitate future improvements in StISSM to allow for non-Gaussian variability or other statistical time series models, as mentioned above. Second, there is a computational limitation: sampling irreducible uncertainty requires ensemble runs with a statistically representative number of ensemble members. As ice-sheet-wide ISM simulations remain computationally expensive, ensemble runs thus require a compromise between number of members and model resolution or number of simulated processes. Third, while we have implemented stochastic capabilities for subglacial water pressure and verified their functionality in test experiments, this study does not discuss simulations with stochasticity imposed on this variable. Realistic benchmark model experiments first require work on choosing an appropriate sliding law and constraining water pressure variability in time. Similarly, while the autoregressive models for climate forcing are flexible approximations of climate variability, they need to be calibrated against climate records or model outputs (Chylek et al.2012; Castruccio et al.2014). Finally, a major challenge is the validation of a stochastic parameterization. Current ISMs calibrate model parameters to match observational datasets that are much shorter than the response timescale of ice sheets. As such they may compensate for not representing the effects of variability in some forcings by using biased parameter values. For example, the noise-induced drift effect in our IQIS results using variability in calving rates could be compensated for by a biased estimate of the commonly used ice tensile strength parameter, which is generally tuned to values much lower than minimal values from field- and laboratory-derived measurements (Petrovic 2003; Amaral et al.2020; Ultee et al.2020; Choi et al.2021; Hillebrand et al.2022).

Our results show that ice sheets need a long period of time to converge to a statistical steady state in the presence of noisy forcing. This raises the question of how zero-mean variability during spin-up could affect the ensuing modeled initial state of an ice sheet. We have used realistic inter-annual and decadal variability in SMB and ocean forcing for our GrIS experiment, reasonably representative of the noisy climate under which the Greenland ice sheet evolves. After 500 years of our synthetic experiment, this resulted in a ±2σ spread equivalent to 1.3 % of the initial mass and a decrease in the mean mass of −0.1 %. Spin-up procedures assuming that climate forcing is constant in time, or without high-frequency noise, lead to subsequent changes in transient experiments that are partially caused by the ice sheet adapting to the abrupt presence of higher-frequency variability.

6 Conclusions

This study has described the development, implementation, and testing of StISSM v1.0, the first stochastic large-scale ice sheet model. Variables with implemented stochastic parameterizations in this first version encompass climate forcing and glaciological processes that are unresolved at the spatiotemporal resolution of ice sheet models: SMB, ocean forcing, calving, and subglacial water pressure. Stochastic climate forcing captures the irreducible uncertainty in climate predictions and how it translates into projected ice sheet mass balance uncertainty. Using stochastic parameterizations for unresolved glaciological processes facilitates the quantification of the impacts of internal variability in such processes on ice dynamics. StISSM v1.0 also includes built-in statistical models for generation of stochastic variability in SMB and oceanic forcing, represented as autoregressive processes. The statistics of the stochastic variability and of the autoregressive climate models are prescribed by users and can thus be adjusted to particular user needs.

We have tested the stochastic capabilities in idealized, synthetic model experiments. These tests have demonstrated that the stochastic parameterizations are functional, and can be upscaled to realistic ice sheet configuration. Our results show that stochastic forcings cause responses of the ice sheet system in line with those observed in stochastic climate and ocean model experiments. For example, stochastic forcing causes not only variability in the final state, but also non-zero tendencies in the response, noise-induced drift, and long timescales needed for ice sheet state convergence. Even in the simple experiments proposed here, the features of the response are complex and cannot be quantified without running ensemble simulations. The response of a particular system is sensitive to the type of forcing, to the geometric configuration, and to the intrinsic non-linearity of ice dynamics. Our results thus raise important questions about representing fluctuating processes with constant deterministic parameterizations, about neglecting high-frequency climatic noise, and about ice sheet model initialization performed without imposing variability.

Our strategy for the development of StISSM v1.0 allows for potential future extensions of stochastic capabilities to other variables in a straightforward manner. In the future, calibration work will be needed to constrain the statistical models for climate forcing, as well as the variability in unresolved glaciological processes such as calving and hydrology. Such an effort will require combining observations, theory, and results from high-fidelity model experiments to understand the internal spatiotemporal variability in processes of interest. Our implementation allows for any spatial, temporal, and intervariable correlation features. StISSM v1.0 thus provides a robust modeling framework to quantify the impacts of forcings with internal variability on ice sheet mass balance.

Appendix A: GrIS spin-up state

In this section, we briefly provide additional details concerning the GrIS configuration at the end of the spin-up (Sect. 3.3). The fields of ice thickness and ice velocities are displayed in Fig. A1c and d. These can be compared to the initial GrIS configuration matched to observational datasets (Sect. 3.3) in Fig. A1a and b. The differences in ice thickness and velocities between both configurations are shown in Fig. A1e and f. We note three main differences: (i) a major thinning at the margins of the Northeast Greenland Ice Stream (NEGIS), (ii) thickening in the southwest, and (iii) velocity changes at the simulated outlet glaciers typically ranging between −200 and +200 m yr−1. The thickness differences across the GrIS are mostly caused by the geometrical adjustment of the ice sheet away from its observed state, which in turn drives SMB feedback processes via the lapse rates applied (Table 2). In contrast, geometrical changes due to frontal ablation and migration allowed in the second phase of the spin-up (see Sect. 3.3) are small because the calving rates have been calibrated for this purpose (Table 3). At the margins of the NEGIS, thinning is caused by an imbalance between ice flow and initial thickness, which in turn causes the SMB to decrease. There is a noticeable increase in ice velocities just upstream of the margin of NEGIS (Fig. A1f), but with a minor impact on ice flow due to the strong thinning in that area (Fig. A1e).

Figure A1Initial configuration from observational datasets: (a) ice thickness and (b) ice velocities. Configuration of the steady state at the end of the spin-up: (c) ice thickness and (d) ice velocities. Total differences between the steady-state and the initial configurations: (e) ice thickness (c–a) and (f) ice velocities (d–b).

Figure A2Differences in (a) ice thickness and (b) ice velocities over the last 200 years of the spin-up to a GrIS steady state. Arabic numerals (black) show the individual basins. Roman numerals (cyan) show the outlet glaciers where ice front movement is simulated and ocean melt parameterized.

In Fig. A2, we show the changes in ice thickness and velocities over the last 200 years of spin-up. Note that the second phase of the spin-up, in which we resolve the 11 outlet glaciers, has a total duration of 1000 years; Fig. A2 displays the changes between year 800 and year 1000 of that phase. The changes are relatively small, with maximum magnitudes in thickness change and velocity change of ∼40 m (rate of 0.2 m yr−1) and ∼35 m yr−1 (rate of 1.8 m yr−2), respectively. Only glacier XI (Heilprin) still shows small patterns of frontal migration and thinning after 800 years of simulation with frontal migration activated (Fig. A2). The drift due to imperfect equilibrium of the initial state is quantified with a control run (Sect. 3.3). We consider this drift sufficiently small (Sect. 4.3), and thus with negligible impact on the transient results driven by stochasticity. A larger drift would obscure the interpretation of results because it would require the assumption that changes due to drift and stochastic variability in forcings can be separated linearly.

Appendix B: Correlation matrix for the GrIS transient ensemble

The correlation matrix relating noise terms for both SMB and TF in all basins of the GrIS stochastic transient runs is specified by Eq. (B1). In Eq. (B1), subscripts i and j denote basins that constitute individual spatial dimensions in the covariance matrix. From the correlation matrix given by Eq. (B1), the covariance matrix is subsequently obtained by left- and right-multiplying it with the diagonal matrix of the standard deviations (see Eq. 13).

(B1) r ϵ SMB , i , ϵ SMB , j = 0.6  if  i  is neighbor of  j r ϵ TF , i , ϵ TF , j = 0.6  if  i  is neighbor of  j r ϵ SMB , i , ϵ TF , j = - 0.6  if  i = j r ϵ SMB , i , ϵ SMB , j = 0.5  if  i  is not a neighbor of  j r ϵ TF , i , ϵ TF , j = 0.5  if  i  is not a neighbor of  j r ϵ SMB , i , ϵ TF , j = - 0.4  if  i j .
Appendix C: Convergence of the statistics

Figure C1Convergence of the (a) mean, (b) standard deviation, and (c) skew in final ice mass for the MISMIP+ transient ensembles as a function of number of ensemble members. For (b), we show the standard deviation relative to the standard deviation of the full ensemble with 500 members.


Figure C2Convergence of the (a) mean, (b) standard deviation, and (c) skew in final ice mass for the IQIS transient ensembles as a function of number of ensemble members. For (b), we show the standard deviation relative to the standard deviation of the full ensemble with 500 members.


Figure C3Convergence of the (a) mean, (b) standard deviation, and (c) skew in final ice mass for the GrIS transient ensemble as a function of number of ensemble members. For (b), we show the standard deviation relative to the standard deviation of the full ensemble with 200 members.


In this section, we analyze how the statistics of interest converge as the number of members per ensemble increases. The statistics of interest are the mean, the standard deviation, and the skew in final ice mass. The analyses of convergence for the MISMIP+, IQIS, and GrIS ensembles are shown in Figs. C1, C2, and C3, respectively. To generate these figures, an initial random sample of five ensemble members is selected, and the statistics for this sample are computed. Iteratively, we add a random sample of five members to this initial sample and compute the statistics at each iteration. The process is performed until all the ensemble members are included in the analysis. A qualitative evaluation of the figures shows that statistics exhibit variability generally until 100 to 150 members are included. Beyond 150 members, the statistics show adequate convergence. We show σE relative to σE of the full corresponding ensemble because variability in absolute σE values depends on the full ensemble σE. As mentioned in Sect. 5, we emphasize that the number of members needed for convergence will vary depending on the ice sheet system investigated, on the timescale of stochastic variability, and on the statistics of interest.

Code and data availability

The stochastic schemes evaluated here are currently implemented in the public release of ISSM. The code can be downloaded, compiled, and executed following the instructions available on the ISSM website: (last access: 21 May 2022). The public SVN repository for the ISSM code can also be found directly at (Larour et al.2020) and downloaded using username “anon” and password “anon”. The version of the code for this study, corresponding to ISSM release 4.19, is SVN version tag number 27017. The documentation of the code version used here is available at (last access: 21 May 2022). The ISSM code, including all the StISSM v1.0 capabilities, is also available as a Zenodo dataset: (ISSM Team2022). The simulation results, the scripts to reproduce all the figures, and the scripts to perform statistical analyses as well as all the input files and preprocessing, run control, and postprocessing scripts to reproduce the simulations are available as a Zenodo dataset: (Verjans2022).

Author contributions

VV led the model development, performed the model experiments, and led the writing of the manuscript. AAR supervised the work. HS contributed to the model development. VV, AAR, HS, and AFT conceived the study. LU and VV have worked on the formulation and calibration of the statistical models. All authors provided comments and suggested edits to the manuscript. All authors, and this study, are part of the Stochastic Ice Sheet Project, aimed at understanding ice sheet sensitivity to variability.

Competing interests

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 in published maps and institutional affiliations.


Helene Seroussi was also funded by the NSF Navigating the New Arctic program. Computing resources were provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology, Atlanta. We acknowledge HPC assistance from Fang (Cherry) Liu. We thank Mathieu Morlighem and Justin Quinn for providing helpful advice about ISSM. We thank Kevin Bulthuis for the implementation of the random number generator. We thank all developers of ISSM for their continuing work on model development. We acknowledge the two anonymous reviewers for their constructive comments that helped improve the quality of the manuscript. Vincent Verjans thanks John Christian for his interest in the study and for insightful discussions about ice sheet sensitivity to variability.

Financial support

This research has been supported by the Heising-Simons Foundation (grant no. 2020-1965).

Review statement

This paper was edited by Christopher Horvat and reviewed by two anonymous referees.


Albrecht, T. and Levermann, A.: Fracture-induced softening for large-scale ice dynamics, The Cryosphere, 8, 587–605,, 2014. a

Amaral, T., Bartholomaus, T. C., and Enderlin, E. M.: Evaluation of Iceberg Calving Models Against Observations From Greenland Outlet Glaciers, J. Geophys. Res.-Earth, 125, 1–29,, 2020. a

Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone, R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland, P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.: Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1), Geosci. Model Dev., 9, 2471–2497,, 2016. a, b, c

Aschwanden, A., Fahnestock, M. A., Truffer, M., Brinkerhoff, D. J., Hock, R., Khroulev, C., Mottram, R., and Khan, S. A.: Contribution of the Greenland Ice Sheet to sea level over the next millennium, Sci. Adv., 5, 6,, 2019. a, b, c

Åström, J., Cook, S., Enderlin, E., Sutherland, D., Mazur, A., and Glasser, N.: Fragmentation theory reveals processes controlling iceberg size distributions, J. Glaciol., 67, 603–612,, 2021. a

Bao, J., McInerney, D. J., and Stein, M. L.: A spatial-dependent model for climate emulation, Environmetrics, 27, 396–408,, 2016. a

Bassis, J. N.: The statistical physics of iceberg calving and the emergence of universal calving laws, J. Glaciol., 57, 3–16, 2011. a, b, c, d

Benn, D. I., Warren, C. R., and Mottram, R. H.: Calving processes and the dynamics of calving glaciers, Earth-Sci. Rev., 82, 143–179,, 2007. a

Berends, C. J., Goelzer, H., and van de Wal, R. S. W.: The Utrecht Finite Volume Ice-Sheet Model: UFEMISM (version 1.0), Geosci. Model Dev., 14, 2443–2470,, 2021. a

Berner, J., Achatz, U., Batté, L., Bengtsson, L., Cámara, A. D. L., Christensen, H. M., Colangeli, M., Coleman, D. R., Crommelin, D., Dolaptchiev, S. I., and Franzke, C. L.: Stochastic parameterization: Toward a new view of weather and climate models, B. Am. Meteorol. Soc., 98, 565–588, 2017. a, b, c

Bondzio, J. H., Seroussi, H., Morlighem, M., Kleiner, T., Rückamp, M., Humbert, A., and Larour, E. Y.: Modelling calving front dynamics using a level-set method: application to Jakobshavn Isbræ, West Greenland, The Cryosphere, 10, 497–510,, 2016. a

Brondex, J., Gillet-Chaulet, F., and Gagliardini, O.: Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law, The Cryosphere, 13, 177–195,, 2019. a

Budd, W. F., Jenssen, D., and Smith, I. N.: A three-dimensional time-dependent model of three West Antarctic ice-sheet, Ann. Glaciol., 5, 29–36, 1984. a

Bulthuis, K., Arnst, M., Sun, S., and Pattyn, F.: Uncertainty quantification of the multi-centennial response of the Antarctic ice sheet to climate change, The Cryosphere, 13, 1349–1380,, 2019. a, b

Bulthuis, K. and Larour, E.: Implementation of a Gaussian Markov random field sampler for forward uncertainty quantification in the Ice-sheet and Sea-level System Model v4.19, Geosci. Model Dev., 15, 1195–1217,, 2022. a

Burke, E. E. and Roe, G. H.: The absence of memory in the climatic forcing of glaciers, Clim. Dynam., 42, 1335–1346,, 2014. a

Castruccio, S. and Stein, M. L.: Global space-time models for climate ensembles, Ann. Appl. Stat., 7, 1593–1611,, 2013. a

Castruccio, S., McInerney, D. J., Stein, M. L., Crouch, F. L., Jacob, R. L., and Moyer, E. J.: Statistical emulation of climate model projections based on precomputed GCM runs, J. Climate, 27, 1829–1844,, 2014. a, b, c

Christian, J. E., Robel, A. A., Proistosescu, C., Roe, G., Koutnik, M., and Christianson, K.: The contrasting response of outlet glaciers to interior and ocean forcing, The Cryosphere, 14, 2515–2535,, 2020. a

Christian, J. E., Robel, A. A., and Catania, G.: A probabilistic framework for quantifying the role of anthropogenic climate change in marine-terminating glacier retreats, The Cryosphere, 16, 2725–2743,, 2022. a, b

Choi, Y., Morlighem, M., Rignot, E., and Wood, M.: Ice dynamics will remain a primary driver of Greenland ice sheet mass loss over the next century, Commun. Earth. Environ., 2, 26,, 2021. a

Christensen, H. M.: Constraining stochastic parametrisation schemes using high-resolution simulations, Q. J. Roy. Meteor. Soc., 146, 938–962,, 2020. a, b

Chylek, P., Folland, C., Frankcombe, L., Dijkstra, H., Lesins, G., and Dubey, M.: Greenland ice core evidence for spatial and temporal variability of the Atlantic Multidecadal Oscillation, Geophys. Res. Lett., 39, L09705,, 2012. a, b, c

Cornford, S. L., Martin, D. F., Graves, D. T., Ranken, D. F., Le Brocq, A. M., Gladstone, R. M., Payne, A. J., Ng, E. G., and Lipscomb, W. H.: Adaptive mesh, finite volume modeling of marine ice sheets, J. Comput. Phys., 232, 529–549,, 2013. a

Cornford, S. L., Seroussi, H., Asay-Davis, X. S., Gudmundsson, G. H., Arthern, R., Borstad, C., Christmann, J., Dias dos Santos, T., Feldmann, J., Goldberg, D., Hoffman, M. J., Humbert, A., Kleiner, T., Leguy, G., Lipscomb, W. H., Merino, N., Durand, G., Morlighem, M., Pollard, D., Rückamp, M., Williams, C. R., and Yu, H.: Results of the third Marine Ice Sheet Model Intercomparison Project (MISMIP+), The Cryosphere, 14, 2283–2301,, 2020. a

Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, Butterworth-Heinemann, Oxford,, 2010. a

dos Santos, T. D., Morlighem, M., Seroussi, H., Devloo, P. R. B., and Simões, J. C.: Implementation and performance of adaptive mesh refinement in the Ice Sheet System Model (ISSM v4.14), Geosci. Model Dev., 12, 215–232,, 2019. a

Dias dos Santos, T., Morlighem, M., and Brinkerhoff, D.: A new vertically integrated MOno-Layer Higher-Order (MOLHO) ice flow model, The Cryosphere, 16, 179–195,, 2022. a

Edwards, T. L., Fettweis, X., Gagliardini, O., Gillet-Chaulet, F., Goelzer, H., Gregory, J. M., Hoffman, M., Huybrechts, P., Payne, A. J., Perego, M., Price, S., Quiquet, A., and Ritz, C.: Probabilistic parameterisation of the surface mass balance–elevation feedback in regional climate model simulations of the Greenland ice sheet, The Cryosphere, 8, 181–194,, 2014. a

Emetc, V., Tregoning, P., Morlighem, M., Borstad, C., and Sambridge, M.: A statistical fracture model for Antarctic ice shelves and glaciers, The Cryosphere, 12, 3187–3213,, 2018. a

Enderlin, E. M., Howat, I. M., Jeong, S., Noh, M.-J., van Angelen, J. H., and van den Broeke, M. R.: An improved mass budget for the Greenland ice sheet, Geophys. Res. Lett., 41, 866–872,, 2014. a

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 high-resolution climate modeling, Geophys. Res. Lett., 36, L12501,, 2009. a

Farrell, B. F. and Ioannou, P. J.: Stochastic dynamics of the midlatitude atmospheric jet, J. Atmos. Sci., 52, 1642–1656,<1642:SDOTMA>2.0.CO;2, 1995. a

Favier, L., Durand, G., Cornford, S. L., Gudmundsson, G. H., Gagliardini, O., Gillet-Chaulet, F., Zwinger, T., Payne, A. J., and Le Brocq, A. M.: Retreat of Pine Island Glacier controlled by marine ice-sheet instability, Nat. Clim. Change, 4, 117–121,, 2014. a

Fettweis, X., Hofer, S., Krebs-Kanzow, 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,, 2020. a, b

Franzke, C., Berner J., Lucarini V., OKane, T. J., and Williams P. D.: Stochastic climate theory and modeling, WIRES Clim. Change, 6, 6378,, 2015. a

Fyke, J., Sergienko, O., Löfverström, M., Price, S., and Lenaerts, J. T.: An overview of interactions and feedbacks between ice sheets and the Earth system, Rev. Geophys., 56, 361–408, 2018. a

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318,, 2013. a

Gladstone, R., Galton-Fenzi, B., Gwyther, D., Zhou, Q., Hattermann, T., Zhao, C., Jong, L., Xia, Y., Guo, X., Petrakopoulos, K., Zwinger, T., Shapero, D., and Moore, J.: The Framework For Ice Sheet–Ocean Coupling (FISOC) V1.1, Geosci. Model Dev., 14, 889–905,, 2021. a

Goelzer, H., Huybrechts, P., Fürst, J. J., Andersen, M. L., Edwards, T. L., Fettweis, X., Nick, F. M., Payne, A. J., and Shannon, S. R.: Sensitivity of Greenland ice sheet projections to model formulations, J. Glaciol., 59, 733–749,, 2013. a

Goelzer, H., Nowicki, S., Edwards, T., Beckley, M., Abe-Ouchi, A., Aschwanden, A., Calov, R., Gagliardini, O., Gillet-Chaulet, F., Golledge, N. R., Gregory, J., Greve, R., Humbert, A., Huybrechts, P., Kennedy, J. H., Larour, E., Lipscomb, W. H., Le clec'h, S., Lee, V., Morlighem, M., Pattyn, F., Payne, A. J., Rodehacke, C., Rückamp, M., Saito, F., Schlegel, N., Seroussi, H., Shepherd, A., Sun, S., van de Wal, R., and Ziemen, F. A.: Design and results of the ice sheet model initialisation experiments initMIP-Greenland: an ISMIP6 intercomparison, The Cryosphere, 12, 1433–1460,, 2018. a

Goelzer, H., Noël, B. P. Y., Edwards, T. L., Fettweis, X., Gregory, J. M., Lipscomb, W. H., van de Wal, R. S. W., and van den Broeke, M. R.: Remapping of Greenland ice sheet surface mass balance anomalies for large ensemble sea-level change projections, The Cryosphere, 14, 1747–1762,, 2020a. a

Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., Abe-Ouchi, 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 sea-level contribution of the Greenland ice sheet: a multi-model ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096,, 2020b. a, b

Goldberg, D. N. and Sergienko, O. V.: Data assimilation using a hybrid ice flow model, The Cryosphere, 5, 315–327,, 2011. a

Haseloff, M. and Sergienko, O.: The effect of buttressing on grounding line dynamics, J. Glaciol., 64, 417–431,, 2018. a

Hasselmann, K.: Stochastic Climate Models. 1. Theory, Tellus, 28, 473–485,, 1976. a, b

Hasselmann, K.: PIPs and POPs: The reduction of complex dynamical systems using principal interaction and oscillation patterns, J. Geophys. Res., 93, 11015–11021, 1988. a

Hewitt, I.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sc. Lett., 371, 16–25, 2013. a, b, c

Hillebrand, T. R., Hoffman, M. J., Perego, M., Price, S. F., and Howat, I. M.: The contribution of Humboldt Glacier, North Greenland, to sea-level rise through 2100 constrained by recent observations of speedup and retreat, The Cryosphere Discuss. [preprint],, in review, 2022. a

Hindmarsh, R. C. A. and Le Meur, E.: Dynamical processes involved in the retreat of marine ice sheets, J. Glaciol., 47, 271–282, 2001. a, b

Hoffman, M. J., Perego, M., Price, S. F., Lipscomb, W. H., Zhang, T., Jacobsen, D., Tezaur, I., Salinger, A. G., Tuminaro, R., and Bertagna, L.: MPAS-Albany Land Ice (MALI): a variable-resolution ice sheet model for Earth system modeling using Voronoi grids, Geosci. Model Dev., 11, 3747–3780,, 2018. a

Hoffman, M. J., Asay-Davis, X., Price, S. F., Fyke, J., and Perego, M.: Effect of Subshelf Melt Variability on Sea Level Rise Contribution From Thwaites Glacier, Antarctica, J. Geophys. Res.-Earth, 124, 2798–2822,, 2019. a

Hu, W. and Castruccio, S.: Approximating the Internal Variability of Bias-Corrected Global Temperature Projections with Spatial Stochastic Generators, J. Climate, 34, 8409–8418,, 2021. a, b, c

Huybrechts, P., Goelzer, H., Janssens, I., Driesschaert, E., Fichefet, T., Goosse, H., and Loutre, M. F.: Response of the Greenland and Antarctic Ice Sheets to Multi-Millennial Greenhouse Warming in the Earth System Model of Intermediate Complexity LOVECLIM, Surv. Geophys., 32, 397–416,, 2011. a

ISSM Team: Ice-sheet and Sea-level System Model source code, v4.21 r27238 (4.21), Zenodo [code],, 2022. a

Jenkins, A., Shoosmith, D., Dutrieux, P., Jacobs, S., Kim, T. W., Lee, S. H., Ha, H. K., and Stammerjohn, S.: West Antarctic Ice Sheet retreat in the Amundsen Sea driven by decadal oceanic variability, Nat. Geosci., 11, 733–738,, 2018. a

Joughin, I. A. N., Smith, B. E., and Howat, I. M.: A complete map of Greenland ice velocity derived from satellite data collected over 20 years, J. Glaciol., 64, 1–11,, 2017. a

Kay, J. E., Deser, C., Phillips, A., Mai, A., Hannay, C., Strand, G., Arblaster, J. M., Bates, S., Danabasoglu, G., Edwards, J., Holland, M., Kushner, P., Lamarque, J.-F., Lawrence, D., Lindsay, K., Middleton, A., Munoz, E., Neale, R., Oleson, K., Polvani, L., and Vertenstein, M.: The Community Earth System Model (CESM) large ensemble project: A community resource for studying climate change in the presence of internal climate variability, B. Am. Meteorol. Soc., 96, 1333–1349,, 2015. a

Kingslake, J.: Chaotic dynamics of a glaciohydraulic model, J. Glaciol., 61, 493–502,, 2015. a

Knuth, D.: The art of computer programming, vol. 1, Fundamental Algorithms, 3rd edn., Addison-Wesley, Reading MA, ISBN 0-201-89683-4, 1998. 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., 117, F01022,, 2012. a, b, c, d

Larour, E., Morlighem, M., and Seroussi, H.: Ice-Sheet and Sea-Level System Model, svn repository [data set], (last access: 21 May 2022), 2020. a

Lipscomb, W. H., Price, S. F., Hoffman, M. J., Leguy, G. R., Bennett, A. R., Bradley, S. L., Evans, K. J., Fyke, J. G., Kennedy, J. H., Perego, M., Ranken, D. M., Sacks, W. J., Salinger, A. G., Vargo, L. J., and Worley, P. H.: Description and evaluation of the Community Ice Sheet Model (CISM) v2.1, Geosci. Model Dev., 12, 387–424,, 2019. a

MacAyeal, D. R.: Large-scale ice flow over a viscous basal sediment: Theory and application to Ice Stream B, Antarctica, J. Geophys. Res., 94, 4071–4087, 1989. a

Maher, N., Milinski, S., Suarez-Gutierrez, L., Botzet, M., Kornblueh, L., Takano, Y., Kröger, J., Ghosh, R., Hedemann, C., Li, C., Li, H., Manzini, E., Notz, D., Putrasahan, D., Boysen, L., Claussen, M., Ilyina, T., Olonscheck, D., Raddatz, T., Stevens, B., and Marotzke, J.: The Max Planck Institute Grand Ensemble – enabling the exploration of climate system variability, J. Adv. Model. Earth Syst., 11, 2050–2069,, 2019. a

Majda, A. J., Timofeyev, I., and Eijnden, E. V.: A mathematical framework for stochastic climate models, Commun. Pur. Appl. Math., 54, 891–974,, 2001. a

Mantelli, E., Bertagni, M. B., and Ridolfi, L.: Stochastic ice stream dynamics, P. Natl. Acad. Sci. USA, 113, E4594–E4600,, 2016. a

Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B. (Eds.): Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge Univerity Press, Cambridge,, 2021. a

McKinnon, K. A., Poppick, A., Dunn-Sigouin, E., and Deser, C.: An `Observational Large Ensemble' to compare observed and modeled temperature trend uncertainty due to internal variability, J. Climate, 30, 7585–7598,, 2017. a

Mikkelsen, T. B., Grinsted, A., and Ditlevsen, P.: Influence of temperature fluctuations on equilibrium ice sheet volume, The Cryosphere, 12, 39–47,, 2018. a, b, c, d, e

Mitchell, J. M.: An overview of climate variability and its causal mechanisms, Quaternary Res., 6, 481–493, 1976. a

Morlighem, M., Rignot, E., Seroussi, H., Larour, E., Dhia, H. B., and Aubry, D.: Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica, Geophys. Res. Lett., 37, 6,, 2010. a

Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., Noël, B. P. Y., Cofaigh, C. Ó., Palmer, S., Rysgaard, S., Seroussi, H., Siegert, M. J., Slabon, P., Straneo, F., van den Broeke, M. R., Weinrebe, W., Wood, M., and Zinglersen, K. B.: BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multi-beam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, 11051–11061,, 2017. a

Mudelsee, M.: Climate Time Series Analysis: Classical Statistical and Bootstrap Methods, Springer, New York, USA,, 2014. a

Nowicki, S. and Seroussi, H.: Projections of future sea level contributions from the Greenland and Antarctic Ice Sheets: Challenges beyond dynamical ice sheet modeling, Oceanography, 31, 109–117,, 2018. a

Nowicki, S., Goelzer, H., Seroussi, H., Payne, A. J., Lipscomb, W. H., Abe-Ouchi, A., Agosta, C., Alexander, P., Asay-Davis, X. S., Barthel, A., Bracegirdle, T. J., Cullather, R., Felikson, D., Fettweis, X., Gregory, J. M., Hattermann, T., Jourdain, N. C., Kuipers Munneke, P., Larour, E., Little, C. M., Morlighem, M., Nias, I., Shepherd, A., Simon, E., Slater, D., Smith, R. S., Straneo, F., Trusel, L. D., van den Broeke, M. R., and van de Wal, R.: Experimental protocol for sea level projections from ISMIP6 stand-alone ice sheet models, The Cryosphere, 14, 2331–2368,, 2020. a

Oerlemans, J.: Holocene glacier fluctuations: is the current rate of retreat exceptional?, Ann. Glaciol., 31, 39–44,, 2000. a

Palmer, T. N.: Stochastic weather and climate models, Nat. Rev. Phys., 1, 463–471,, 2019. a

Pattyn, F.: Sea-level response to melting of Antarctic ice shelves on multi-centennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0), The Cryosphere, 11, 1851–1878,, 2017. a

Pattyn, F., Ritz, C., Hanna, E., Asay-Davis, X., DeConto, R., Durand, G., Favier, L., Fettweis, X., Goelzer, H., Golledge, N. R., Munneke, P. K., Lenaerts, J. T. M., Nowicki, S., Payne, A. J., Robinson, A., Seroussi, H., Trusel, L. D., and van den Broeke, M.: The Greenland and Antarctic ice sheets under 1.5 C global warming, Nat. Clim. Change, 8, 1053–1061,, 2018. a

Perego, M., Price, S., and Stadler, G.: Optimal initial conditions for coupling ice sheet models to Earth system models, J. Geophys. Res., 119, 1894–1917,, 2014. a

Perron, M. and Sura, P.: Climatology of non-Gaussian atmospheric statistics, J. Climate, 26, 1063–1083,, 2013. a

Petra, N., Zhu, H., Stadler, G. Hughes, T. J. R., and Ghattas, O.: An inexact Gauss–Newton method for inversion of basal sliding and rheology parameters in a nonlinear Stokes ice sheet model, J. Glaciol., 58, 889–903, 2012. a

Petrovic, J. J.: Review Mechanical properties of ice and snow, J. Mater. Sci., 38, 1–6, 2003. a

Pollard, D. and DeConto, R. M.: A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica, The Cryosphere, 6, 953–971,, 2012. a, b

Porta Mana, P. and Zanna, L.: Toward a stochastic parameterization of ocean mesoscale eddies, Ocean Model., 79, 1–20, 2014. a, b

Razali, N. M., and Wah, Y. B.: Power comparisons of Shapiro-Wilk, Kolmogorov-Smirnov, Lilliefors and Anderson-Darling tests. J. Statist. Model. Anal., 2, 21–33, 2011. a

Ribergaard, M. H., Olsen, S. M., and Mortensen, J.: Oceanographic Investigations off West Greenland 2007, Danish Metrological Institute Centre for Ocean and Ice (DMI), Copenhagen, 48 pp., (last access: 13 May 2022), 2008. a

Rignot, E., Xu, Y., Menemenlis, D., Mouginot, J., Scheuchl, B., Li, X., Morlighem, M., Seroussi, H., den Broeke, M. V., Fenty, I., Cai, C., An, L., and Fleurian, B. D.: Modeling of ocean-induced ice melt rates of five west Greenland glaciers over the past two decades, Geophys. Res. Lett., 43, 6374–6382,, 2016. a, b

Robel, A. A., Schoof, C., and Tziperman, E.: Rapid grounding line migration induced by internal ice stream variability, J. Geophys. Res.-Earth Surf., 119, 2430–2447,, 2014. a

Robel, A. A., Roe, G. H., and Haseloff, M.: Response of Marine-Terminating Glaciers to Forcing: Time Scales, Sensitivities, Instabilities, and Stochastic Dynamics, J. Geophys. Res.-Earth Surf., 123, 2205–2227,, 2018. a, b, c, d, e

Robel, A. A., Seroussi, H., and Roe, G. H.: Marine ice sheet instability amplifies and skews uncertainty in projections of future sea-level rise, P. Natl. Acad. Sci. USA, 116, 14887–14892,, 2019. a, b, c, d

Roe, G. H. and Baker, M. B.: Glacier response to climate perturbations: An accurate linear geometric model, J. Glaciol., 60, 670–684,, 2014. a, b

Roe, G. H. and Baker, M. B.: The response of glaciers to climatic persistence, J. Glaciol., 62, 440–450, 2016. a, b, c, d, e

Roe, G. H. and O'Neal, M. A.: The response of glaciers to intrinsic climate variability: Observations and models of late-holocene variations in the Pacific northwest, J. Glaciol., 55, 839–854,, 2009. a

Roe, G. H. and Steig, E. J.: Characterization of Millennial-scale Climate Variability, J. Climate, 17, 1929–1944,<1929:COMCV>2.0.CO;2, 2004. a

Rosier, S. H. R., Reese, R., Donges, J. F., De Rydt, J., Gudmundsson, G. H., and Winkelmann, R.: The tipping points and early warning indicators for Pine Island Glacier, West Antarctica, The Cryosphere, 15, 1501–1516,, 2021. a

Schlegel, N.-J., Seroussi, H., Schodlok, M. P., Larour, E. Y., Boening, C., Limonadi, D., Watkins, M. M., Morlighem, M., and van den Broeke, M. R.: Exploration of Antarctic Ice Sheet 100-year contribution to sea level rise and associated model uncertainties using the ISSM framework, The Cryosphere, 12, 3511–3534,, 2018. a, b, c

Seroussi, H. and Morlighem, M.: Representation of basal melting at the grounding line in ice flow models, The Cryosphere, 12, 3085–3096,, 2018. a, b, c

Seroussi, H., Morlighem, M., Rignot, E., Mouginot, J., Larour, E., Schodlok, M., and Khazendar, A.: Sensitivity of the dynamics of Pine Island Glacier, West Antarctica, to climate forcing for the next 50 years, The Cryosphere, 8, 1699–1710,, 2014. a

Seroussi, H., Nowicki, S., Simon, E., Abe-Ouchi, A., Albrecht, T., Brondex, J., Cornford, S., Dumas, C., Gillet-Chaulet, F., Goelzer, H., Golledge, N. R., Gregory, J. M., Greve, R., Hoffman, M. J., Humbert, A., Huybrechts, P., Kleiner, T., Larour, E., Leguy, G., Lipscomb, W. H., Lowry, D., Mengel, M., Morlighem, M., Pattyn, F., Payne, A. J., Pollard, D., Price, S. F., Quiquet, A., Reerink, T. J., Reese, R., Rodehacke, C. B., Schlegel, N.-J., Shepherd, A., Sun, S., Sutter, J., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., and Zhang, T.: initMIP-Antarctica: an ice sheet model initialization experiment of ISMIP6, The Cryosphere, 13, 1441–1471,, 2019. a

Seroussi, H., Nowicki, S., Payne, A. J., Goelzer, H., Lipscomb, W. H., Abe-Ouchi, A., Agosta, C., Albrecht, T., Asay-Davis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., Galton-Fenzi, 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 multi-model ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070,, 2020. a, b

Shapiro, N. M. and Ritzwoller, M. H.: Inferring surface heat flux distributions guided by a global seismic model: particular application to Antarctica, Earth Planet. Sc. Lett., 223, 213–224,, 2004. a

Shapiro, S. S. and Wilk, M. B.: An analysis of variance test for normality (complete samples), Biometrika, 52, 591–611, 1965. a, b

Snow, K., Goldberg, D. N., Holland, P. R., Jordan, J. R., Arthern, R. J., and Jenkins, A.: The response of ice sheets to climate variability, Geophys. Res. Lett., 44, 11878–11885, 2017.  a

Straneo, F. and Heimbach, P.: North Atlantic warming and the retreat of Greenland's outlet glaciers, Nature, 504, 36–43, 2013. a

Thomas, E. R., Dennis, P. F., Bracegirdle, T. J., and Franzke, C.: Ice core evidence for significant 100-year regional warming on the Antarctic Peninsula, Geophys. Res. Lett., 36, L20704,, 2009. a

Tsai, C. Y., Forest, C. E., and Pollard, D.: Assessing the contribution of internal climate variability to anthropogenic changes in ice sheet volume, Geophys. Res. Lett., 44, 6261–6268,, 2017. a

Tsai, C. Y., Forest, C. E., and Pollard, D.: The role of internal climate variability in projecting Antarctica's contribution to future sea-level rise, Clim. Dynam., 55, 1875–1892,, 2020. a

Ultee, L., Meyer, C., and Minchew, B.: Tensile strength of glacial ice deduced from observations of the 2015 eastern Skaftá cauldron collapse, Vatnajökull ice cap, Iceland, J. Glaciol., 66, 1024–1033,, 2020. a

van Angelen, J. H., van den Broeke, M. R., Wouters, B., and Lenaerts, J. T. M.: Contemporary (1960–2012) Evolution of the Climate and Surface Mass Balance of the Greenland Ice Sheet, Surv. Geophys., 35, 1155–1174,, 2014. a

Verjans, V.: Data for “The Stochastic Ice-Sheet and Sea-Level System Model v1.0 (StISSM v1.0)” by Verjans et al., Zenodo [data set],, 2022. a

von Storch, H. and Zwiers, F.: Statistical Analysis in Climate Research, Cambridge University Press, Cambridge, ISBN 0521450713, 1999. a, b

Weertman, J.: On the sliding of glaciers, J. Glaciol., 3, 33–38, 1957. a

Wood, M., Rignot, E., Fenty, I., An, L., Bjørk, A., van den Broeke, M., Cai, C., Kane, E., Menemenlis, D., Millan, R., Morlighem, M., Mouginot, J., Noël, B., Scheuchl, B., Velicogna, I., Willis, J. K., and Zhang, H.: Ocean forcing drives glacier retreat in Greenland, Sci. Adv., 7, eaba7282,, 2021. a, b

Zwally, H. J., Giovinetto, M. B., Beckley, M. A., and Saba, J. L.: Antarctic and Greenland drainage systems, GSFC cryospheric sciences laboratory, (last access: 22 March 2022), 2012. a, b

Short summary
We describe the development of the first large-scale ice sheet model that accounts for stochasticity in a range of processes. Stochasticity allows the impacts of inherently uncertain processes on ice sheets to be represented. This includes climatic uncertainty, as the climate is inherently chaotic. Furthermore, stochastic capabilities also encompass poorly constrained glaciological processes that display strong variability at fine spatiotemporal scales. We present the model and test experiments.