the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Implementation of a Gaussian Markov random field sampler for forward uncertainty quantification in the Ice-sheet and Sea-level System Model v4.19

### Kevin Bulthuis

### Eric Larour

Assessing the impact of uncertainties in ice-sheet models is a major and challenging issue that needs to be faced by the ice-sheet community to provide more robust and reliable model-based projections of ice-sheet mass balance. In recent years, uncertainty quantification (UQ) has been increasingly used to characterize and explore uncertainty in ice-sheet models and improve the robustness of their projections. A typical UQ analysis first involves the (probabilistic) characterization of the sources of uncertainty, followed by the propagation and sensitivity analysis of these sources of uncertainty. Previous studies concerned with UQ in ice-sheet models have generally focused on the last two steps but have paid relatively little attention to the preliminary and critical step of the characterization of uncertainty. Sources of uncertainty in ice-sheet models, like uncertainties in ice-sheet geometry or surface mass balance, typically vary in space and potentially in time. For that reason, they are more adequately described as spatio-(temporal) random fields, which account naturally for spatial (and temporal) correlation. As a means of improving the characterization of the sources of uncertainties for forward UQ analysis within the Ice-sheet and Sea-level System Model (ISSM), we present in this paper a stochastic sampler for Gaussian random fields with Matérn covariance function. The class of Matérn covariance functions provides a flexible model able to capture statistical dependence between locations with different degrees of spatial correlation or smoothness properties. The implementation of this stochastic sampler is based on a notable explicit link between Gaussian random fields with Matérn covariance function and a certain stochastic partial differential equation. Discretization of this stochastic partial differential equation by the finite-element method results in a sparse, scalable and computationally efficient representation known as a Gaussian Markov random field. In addition, spatio-temporal samples can be generated by combining an autoregressive temporal model and the Matérn field. The implementation is tested on a set of synthetic experiments to verify that it captures the desired spatial and temporal correlations well. Finally, we illustrate the interest of this stochastic sampler for forward UQ analysis in an application concerned with assessing the impact of various sources of uncertainties on the Pine Island Glacier, West Antarctica. We find that larger spatial and temporal correlations lengths will both likely result in increased uncertainty in the projections.

©2021. California Institute of Technology. Government sponsorship acknowledged.

Despite large improvements in ice-sheet modelling in recent years, model-based estimates of ice-sheet mass balance remain characterized by large uncertainty. The main sources of uncertainty are associated with limitations related to poorly modelled physical processes, the model resolution, poorly constrained initial conditions, uncertainties in external climate forcing (e.g. surface mass balance and ocean-induced melting), or uncertain input data such as the ice sheet geometry (e.g. bedrock topography and surface elevation) or boundary conditions (e.g. basal friction and geothermal heat flux). In order to provide more robust and reliable model-based estimates of ice-sheet mass balance, we therefore need to understand how model outputs are affected by or sensitive to input parameters.

To this aim, uncertainty quantification (UQ) methods have become a powerful and popular tool to deduce the impact of sources of uncertainty on ice-sheet projections (propagation of uncertainty) or to ascertain and rank the impact of each source of uncertainty on the projection uncertainty (sensitivity analysis) (Bulthuis et al., 2019, 2020; Edwards et al., 2019, 2021; Larour et al., 2012b; Ritz et al., 2015; Schlegel et al., 2018). While the aforementioned studies have mainly focused on the propagation and sensitivity analysis of sources of uncertainty in ice-sheet models, relatively little attention has been paid to the preliminary and critical characterization and description of these sources of uncertainty. Uncertainties in input parameters are generally lumped into uncertain (lumped) parameters that represent one or several sources of parametric and/or unmodelled dynamics uncertainty (Bulthuis et al., 2019; Edwards et al., 2019, 2021; Ritz et al., 2015). These uncertain lumped parameters then appear in parameterizations and reduced-order models of complex physical processes such as basal friction or ocean-induced melting. However, such a characterization does not reflect the fact that most sources of uncertainty in numerical ice-sheet models generally vary in space and potentially in time. With the aim of sampling spatially distributed input variables, Larour et al. (2012b) implemented an approach based on a partitioning of the computational domain, with later applications in Larour et al. (2012a), Schlegel et al. (2013, 2015, 2018), and Schlegel and Larour (2019). In this approach the computational domain is divided into a number of partitions fixed by the user or randomly determined based on the CHACO mesh partitioner (Hendrickson and Leland, 1995). For each partition, a statistical distribution is specified for the field being sampled; thus, sampling is not carried out over the entire field but over the multiple partitions. The number of partitions controls the spatial scale of the random perturbation: the larger the number of the partitions, the more local the perturbation.

In this paper, we propose characterizing uncertain spatially varying input parameters as spatial random fields, that is, an infinite set of random variables indexed by the spatial coordinate. More specifically, we focus on the class of Gaussian random fields, which provides a popular statistical model to represent stochastic phenomena in engineering, spatial analysis, and geostatistics (e.g. kriging interpolation) (Cressie, 1993; Schabenberger and Gotway, 2005). They have also been employed in glaciology in a number of studies, including Isaac et al. (2015), Babaniyi et al. (2021) and Brinkerhoff (2021). Gaussian random fields are entirely defined by their mean (or trend) function and their covariance function. The mean function is typically used to capture the dependance of a stochastic phenomenon at large scales, while the covariance function specifies the spatial dependency between pairs of locations. Among the families of covariance functions, we focus on the Matérn family, which provides a popular and practical choice to represent spatial correlations in geostatistics due to its high flexibility (Stein, 1999). In particular, the Matérn family involves a set of parameters that control the correlation length and smoothness of the realizations of the random field.

The goal of this paper is to introduce the implementation of a stochastic sampler for Gaussian Matérn random fields for forward uncertainty quantification within the Ice-sheet and Sea-level System Model (ISSM) (Larour et al., 2012c) and to illustrate the value of this sampler for forward UQ analysis in ice-sheet models. Typical methods for sampling from Gaussian random fields, for instance methods based on a factorization of the covariance matrix, can be computationally intensive for large-scale problems or not appropriate for unstructured meshes, as can be encountered in applications in glaciology. Here, we implement a Gaussian Matérn random field sampler based on an explicit link between Gaussian random fields with Matérn covariance function and a stochastic partial differential equation (SPDE) (Lindgren et al., 2011). Discretization of this SPDE with the finite-element method (FEM) results in a sparse, scalable and computationally efficient representation of a Gaussian random field known as a Gaussian Markov random field (GRMF). The SPDE approach for Gaussian random fields has already been considered in glaciology, though that was in the context of defining a prior distribution for Bayesian inverse problem in glaciology (Isaac et al., 2015; Petra et al., 2014). In addition, spatio-temporal samples can be generated by combining a temporal autoregressive model and the Matérn random field (Cameletti et al., 2012).

This paper is organized as follows. In Sect. 2, we review the link between Gaussian random fields with Matérn covariance function and a stochastic partial differential equation (SPDE) and describe its discretization with the FEM. In Sect. 3, we test the implementation of the SPDE on a set of synthetic experiments to verify that it captures the desired spatial and temporal correlations well. In Sect. 4, we provide an illustration concerned with assessing the impact of various sources of uncertainties on the Pine Island Glacier (PIG), West Antarctica. Finally, we provide an overall discussion of the results and the approach in Sect. 5.

In this work, we model a spatially varying uncertain input parameter as a random field $\left\{x\left(\mathit{s}\right),\mathit{s}\in D\right\}$ indexed over the
computational domain *D*. In the following, we assume this random field to be a Gaussian random field with zero mean.

## 2.1 Matérn covariance function

A Gaussian random field $\left\{x\left(\mathit{s}\right),\mathit{s}\in D\right\}$ with zero mean is totally characterized by its covariance or kernel function
$C(\mathit{s},{\mathit{s}}^{\prime})$, which is defined mathematically by 𝔼[*x*(** s**)

*x*(

*s*^{′})], where 𝔼 denotes the mathematical expectation. The covariance function describes the statistical relationship (or spatial correlation) between two locations,

**and**

*s*

*s*^{′}, in the domain. Roughly speaking, it describes the similarity of the values taken by the random field at different locations in the domain, with locations close to each other more likely to have similar values. In theory, any positive semi-definite function $C(\mathit{s},{\mathit{s}}^{\prime})$ defines a valid covariance function (Rasmussen and Williams, 2006). In practice, the covariance function is often assumed to be a function of the space lag ${\mathit{s}}^{\prime}-\mathit{s}$ (stationary covariance function) or even the spatial distance $\Vert {\mathit{s}}^{\prime}-\mathit{s}\Vert $ (isotropic covariance function).

Among the families of covariance functions, the Matérn family is a popular choice to represent spatial correlation in geostatistics. Applications include spatial modelling of greenhouse gas emissions (Western et al., 2020), temperature or precipitation anomalies (Furrer et al., 2006), soil properties (Minasny and McBratney, 2005), acoustic wave speeds in seismology (Bui-Thanh et al., 2013), and basal friction in glaciology (Isaac et al., 2015; Petra et al., 2014). The Matérn covariance function is defined as

where *K*_{ν} is the modified Bessel function of the second kind (Abramowitz and Stegun, 1970) and order *ν*>0, Γ is the Gamma function,
*σ*^{2} is the variance of the random field, and *κ* is a positive scaling parameter. The scaling parameter *κ* is more naturally interpreted
as a range parameter *ρ*, with the empirically derived relationship $\mathit{\rho}=\sqrt{\mathrm{8}\mathit{\nu}}/\mathit{\kappa}$ corresponding to correlations near 0.1 at the
distance *ρ* (Lindgren et al., 2011). The order *ν* (or smoothness index) controls the smoothness of the realizations; more specifically, the
realizations are ⌈*ν*−1⌉ times differentiable. Therefore, the order parameter *ν* allows for additional flexibility in spatial
modelling, giving the Matérn family considerable practical value (Stein, 1999). As *ν*→∞, the Matérn covariance function
approaches the squared exponential (also called the Gaussian) form, popular in machine learning, whose realizations are infinitely differentiable. For
*ν*=0.5, the Matérn covariance function takes the exponential form, popular in spatial statistics, which produces continuous but
non-differentiable realizations.

### 2.1.1 Stochastic partial differential equation

In this paper, interest lies in generating realizations of Gaussian random fields for uncertainty quantification in ice-sheet models. Gaussian random
fields have already been employed in glaciology in a number of studies including Isaac et al. (2015), Babaniyi et al. (2021) and
Brinkerhoff (2021). Various methods exist to sample Gaussian random fields (see, e.g. Hristopulos, 2020), yet these methods can be
computationally prohibitive for applications in glaciology. For instance, methods based on a factorization of the covariance matrix typically result
in large dense matrices when the number of degrees of freedom is large. Here, we consider a sampling approach that is based on an explicit link
between Gaussian random fields with Matérn covariance function and a stochastic partial differential equation. Indeed, as noted initially
in Whittle (1954) and more recently revisited in Lindgren et al. (2011), realizations of a Gaussian random field *x*(**s**) with Matérn
covariance function can be obtained as the solution to the fractional stochastic partial differential equation

where *κ* and *τ* are positive parameters, $\mathit{\alpha}=\mathit{\nu}+d/\mathrm{2}$ (where *d* is the spatial dimension), and 𝒲(**s**) is a spatial
Gaussian white noise with unit variance. The parameter *κ* controls the range parameter *ρ*, and the parameters *κ* and *τ* control the
variance of the random field. More explicitly, we have the following relationships (see, e.g. Khristenko et al., 2019):

An intuitive interpretation for the rationale behind the SPDE (Eq. 2) can be found in
the context of the filter theory. For *α*=2, the differential operator (*κ*^{2}−Δ) can be interpreted as a low-pass filter or diffusion
operator that smooths out a Gaussian white noise (no spatial correlation) into a Gaussian random field that exhibits spatial correlation. While a
Gaussian white noise is characterized by a flat wavenumber spectrum, a filtered random field exhibits a decaying wavenumber spectrum, with *κ*^{2}
acting as a cutoff parameter.

Equation (2) requires a choice of boundary conditions. Here, we consider either homogeneous Neumann boundary conditions

or Robin boundary conditions

where *β* is the Robin coefficient and ** n** is the external unit vector to the boundary. We mention that the choice of boundary conditions
typically introduces boundary effects that lead to increased or decreased correlation and pointwise variance close to the boundary; see
Daon et al. (2018) and Khristenko et al. (2019), for instance, for a discussion and approaches to reduce the impact of boundary effects. For

*α*=2 and Robin boundary conditions, the choice of $\mathit{\beta}=\left(\mathit{\kappa}\mathit{\tau}\right)/\mathrm{1.42}$ has been found to be a good choice to mitigate boundary effects (Roininen et al., 2014).

The SPDE approach has several advantages compared to other sampling methods. First, Eq. (2) can be solved efficiently by using standard
discretization schemes (e.g. the finite-element method) and sparse matrix linear algebra (see Sect. 2.1.2). Second, the SPDE
approach provides a flexible model to build more general probabilistic models of random fields that do not require us to impose the covariance structure
explicitly; see Simpson et al. (2011) for a discussion. For instance, the parameters *κ* and *τ* in Eq. (2) can vary spatially in
order to represent non-stationary Gaussian random fields with a spatially varying range parameter and marginal variance. Similarly, the Robin
coefficient can vary spatially to mitigate boundary effects (Daon et al., 2018). Other extensions of the SPDE approach include anisotropic Gaussian
models (Lindgren et al., 2011) or even non-Gaussian models (Bolin, 2014). Finally, the SPDE approach can be extended to define Gaussian random fields
with Matérn covariance function on general spatial domains such as the sphere (Lindgren et al., 2011).

The SPDE approach can also be used to define a proper choice of a prior distribution for inverse problems in infinite dimension (Bui-Thanh et al., 2013; Isaac et al., 2015; Petra et al., 2014; Stuart, 2010). However, this approach does not provide any immediate solution to represent the posterior covariance. For general inverse problems, the posterior distribution does not need to be Gaussian even if the prior distribution is a Gaussian random field. In this case, the posterior covariance can be estimated using, for instance, Markov chain Monte Carlo algorithms (Beskos et al., 2017; Petra et al., 2014) or a Laplace approximation of the posterior distribution (Bui-Thanh et al., 2013; Isaac et al., 2015).

### 2.1.2 Numerical implementation

We assume at this stage that *α*=2. The SPDE (Eq. 2) can be solved efficiently by the finite-element
method; see, for instance, Lindgren et al. (2011) and Chu and Guilleminot (2019). The finite-element representation of the solution is written as

where ${\left\{{\mathit{\psi}}_{k}\right\}}_{k=\mathrm{1}}^{N}$ is a finite-element basis (here piecewise linear functions), *N* is the total number of nodes and
${\left\{{x}_{k}\right\}}_{k=\mathrm{1}}^{N}$ is the stochastic nodal values. The discretized weak form of the differential operator in Eq. (2) with
Robin boundary conditions involves the matrix **K** with entries

where the test functions are the same as the basis functions and ∂*D* denotes the boundary of the domain *D*. The finite-element
discretization of the SPDE also involves the mass matrix **M**, whose entries are given by

As shown in Lindgren et al. (2011), the vector ** x** of nodal values of the FEM solution to the SPDE (Eq. 2)
is a Gaussian random vector that satisfies

where the covariance matrix is given by

where the superscript underpins the order 2. Let ${\mathbf{\Sigma}}^{\left(\mathrm{2}\right)}={\mathbf{L}}^{\left(\mathrm{2}\right)}{{\mathbf{L}}^{\left(\mathrm{2}\right)}}^{T}$ be a matrix factorization of
**Σ**^{(2)} (e.g. Cholesky or QR decomposition). Following this, realizations of ** x** are obtained as

where ** ξ** is a random vector whose entries are independent standard Gaussian random variables. Equivalently, one can show that realizations of

**can obtained by solving the following linear system of equations**

*x*More generally, the covariance matrix **Σ**^{(α)} for any arbitrary integer order *α* is given by (Lindgren et al., 2011)

where ${\mathbf{\Sigma}}^{\left(\mathrm{1}\right)}={\mathbf{K}}^{-\mathrm{1}}$ and ${\mathbf{\Sigma}}^{\left(\mathrm{2}\right)}={\mathbf{K}}^{-\mathrm{1}}{\mathbf{MK}}^{-\mathrm{1}}$. The matrix
**L**^{(α)} for any arbitrary integer order *α* is also given by

where ${\mathbf{L}}^{\left(\mathrm{1}\right)}={\mathbf{K}}^{\mathrm{1}/\mathrm{2}}$ and ${\mathbf{L}}^{\left(\mathrm{2}\right)}={\mathbf{K}}^{-\mathrm{1}}{\mathbf{M}}^{\mathrm{1}/\mathrm{2}}$. Therefore, samples for any arbitrary integer
order *α*>2 can be generated in a recursive way. The algorithm for sampling Gaussian random fields with Matérn covariance function based on
the SPDE approach is summarized in Algorithm 1. Although not implemented in our work, we mention the work by Bolin et al. (2018) that
extends the numerical solution of the SPDE (Eq. 2) to non-integer orders.

The bulk of the computational cost is in evaluating the square root of the matrix **M** or **K**. Even if the matrices **M**
and **K** are sparse, their square root does not need to be sparse. In order to speed up the computation and retain sparsity, the mass
matrix **M** (the same for **K** if ${\mathbf{K}}^{\mathrm{1}/\mathrm{2}}$ needs to be evaluated) can be approximated as a diagonal lump mass
matrix $\stackrel{\mathrm{\u0303}}{\mathbf{M}}$ with diagonal entries

Numerically, this approximation leads to a representation of the precision matrix, i.e. the inverse of the covariance matrix, that is sparse (Lindgren et al., 2011); elements in the precision matrix are non-zero only for diagonal elements and neighbouring nodes. Mathematically, such a sparse representation of the precision matrix corresponds to a Gaussian Markov random field (GRMF) (Rue and Held, 2005). This means that the conditional expectation of the random field at a given node given all the other nodes only depends on the neighbouring nodes. The Markov approximation of the Gaussian random field has been shown to have negligible differences with the exact finite-element representation (Bolin and Lindgren, 2009), and its convergence rate is the same as the exact finite-element representation (Lindgren et al., 2011).

### 2.1.3 Extension to spatio-temporal random fields

Temporal correlation between samples (for transient simulation) is represented with a first-order autoregressive model (AR1 process) (Cameletti et al., 2012; Western et al., 2020):

where *t* denotes a discrete time, *ϕ* the temporal correlation (can be positive, negative or null) and *ϵ*_{t}(** s**) an independent
noise term. The first-order autoregressive model specifies that the value of the random field at time

*t*depends linearly on its value at time

*t*−1 and a stochastic term that models unpredictable effects. The parameter

*ϕ*controls the strength (with sign) of the temporal correlation between the values of the random field at times

*t*and

*t*−1. If

*ϕ*=0, there is no dependence between

*x*

_{t}and

*x*

_{t−1}and only the stochastic term contributes to the spatio-temporal random field. By definition, an AR1 model is an example of a Markov process.

At every time step *t*, the noise term *ϵ*_{t}(** s**) is chosen as a Gaussian random field with Matérn covariance function and positive
variance ${\mathit{\sigma}}_{\mathit{\u03f5}}^{\mathrm{2}}$, obtained as the solution of the SPDE (3). If $\left|\mathit{\varphi}\right|<\mathrm{1}$,

*x*

_{t}is a stationary process in time with zero mean and marginal variance

*σ*

^{2}if

*x*

_{0}is a Gaussian random field with zero mean and variance

*σ*

^{2}and the noise variance is chosen as ${\mathit{\sigma}}_{\mathit{\u03f5}}^{\mathrm{2}}={\mathit{\sigma}}^{\mathrm{2}}(\mathrm{1}-{\mathit{\varphi}}^{\mathrm{2}})$. If $\left|\mathit{\varphi}\right|\ge \mathrm{1}$, the process is non-stationary in time, with the case

*ϕ*=1 corresponding to a discrete random walk. The temporal correlation function (or autocorrelation function) is given by

and is therefore stationary in time because it only depends on the time lag *T*. Algorithm (2) provides a pseudo-code to generate spatio-temporal
samples.

The spatio-temporal model combining Eqs. (2) and (16) defines a spatio-temporal Gaussian random field discretized in time. It represents an example of a separable model in which the spatio-temporal covariance function is the product of the spatial and temporal covariance functions. Therefore, the spatial and temporal covariance structures of the model are totally decoupled; the temporal evolution of the input quantity at a given location does not depend on the temporal evolution at other locations. Although it is beyond the scope of this work, non-separable models can also be formulated based on the SPDE approach (Krainski, 2018; Bakka et al., 2020).

As a means of verifying and testing the implementation of our stochastic sampler, we carry out a set of numerical experiments on a synthetic ice
sheet or ice shelf that is 100 km × 100 km in size. The test mesh comprises 8192 elements and 4225 nodes, which is sufficiently large to
avoid any significant discretization error. In these experiments, we consider Matérn random fields with a parameter range of 20 km and a
unit variance. We first verify the spatial covariance structure by solving the SPDE (Eq. 2) for the orders
*α*=2 and 4, then discuss the impact of the boundary conditions on the simulations in a second experiment. We finally perform transient
simulations to verify the temporal covariance structure.

## 3.1 Experiment 1: verification of covariance structure

As a first experiment, we test our implementation by checking that the generated samples are actually drawn from a Gaussian random field with the
desired Matérn covariance function. More specifically, the marginal distributions have to follow Gaussian distributions with zero mean and unit
variance, while the bivariate distributions have to follow bivariate Gaussian distributions with zero mean and covariance given by
Eq. (1). To test for the covariance structure, we generate an ensemble of 10 000 independent and identically distributed
(i.i.d) samples and estimate the marginal and bivariate distributions from these samples. Figures 1
and 2 show the obtained distributions for the locations ${\mathit{s}}_{\mathrm{1}}=(\mathrm{50},\mathrm{50})$ km, ${\mathit{s}}_{\mathrm{2}}=(\mathrm{53.125},\mathrm{50})$ km
and ${\mathit{s}}_{\mathrm{3}}=(\mathrm{59.375},\mathrm{50})$ km for *α*=2 and *α*=4, respectively. In both figures, we observe that the marginal and bivariate
distributions all follow a Gaussian distribution with variances and covariances close to the theoretical expected values. Small discrepancies between
the estimated variances and covariances and their theoretical values can be explained by the numerical approximations (finite-element discretization
and lumped approximations), the impact of boundary conditions (discussed in Experiment 2), and the limited number of samples used for the estimation.

## 3.2 Experiment 2: influence of the boundary conditions

The second experiment investigates the influence of the boundary conditions on the marginal variance of the samples generated following the SPDE
approach. Figure 3 shows the estimated marginal standard deviation for *α*=2 with 10 000 samples considering either
homogeneous Neumann boundary conditions (Fig. 3a) or Robin boundary conditions (Fig. 3b). The value of the Robin coefficient *β* is taken as $\mathit{\kappa}\mathit{\tau}/\mathrm{1.42}$. With homogeneous Neumann boundary conditions, we observe a clear boundary artefact that yields a larger marginal standard deviation close to
the boundary than in the interior. In particular, the standard deviation is 1.4 times larger along the edges of the domain and 2 times larger near the
corners. The boundary artefact becomes almost negligible near a distance of at least *ρ* from the boundary. On the contrary, Robin boundary
conditions with an appropriate choice of the Robin coefficient help mitigate the boundary artefact significantly. Physically, Neumann boundary
conditions can be interpreted as reflecting boundaries, while Robin boundary conditions can be understood as absorbing boundaries, with the Robin
coefficient being an absorption coefficient.

## 3.3 Experiment 3: transient simulation

As a third experiment, we test our transient implementation by checking that the generated samples follow an AR1 process in time. For different values
of the temporal correlation factor *ϕ*, we generate a series of 10 000 spatio-temporal samples and estimate the autocorrelation function at the
centre location $\mathit{s}=(\mathrm{50},\mathrm{50})$ km. Figure 4 represents the estimated autocorrelation functions and provides a
comparison with the theoretical autocorrelation functions given by Eq. (17). For all temporal correlation factors, the estimated
autocorrelation functions are close to their theoretical value. The small discrepancies can also be explained by numerical approximations and the
limited number of samples in the time series.

## 4.1 Experimental setup

We set up a test problem similar to the application in Larour et al. (2012b), in which the authors performed a UQ analysis on the Pine Island Glacier
using a partitioning approach. We refer the reader to this article for additional information about the ice flow model, the datasets, the input
quantities, and the quantities of interest. We consider three spatially varying uncertain input quantities: the ice thickness *H*, the basal drag
coefficient *α* at the ice–bed interface and the ice hardness *B*. The basal drag coefficient is used in a Budd-type friction law, and the ice
hardness is determined from a thermal model (Larour et al., 2012c). As quantities of interest, we consider the mass flux through 13 fluxgates
positioned along the main tributaries of the PIG (see Fig. 5). The mass flux *M*_{i} through gate *i* is given by

where the integral is a line integral along the fluxgate contour 𝒞_{i}, ** v** is the depth-averaged horizontal velocity computed using the
shallow shelf/stream approximation (SSA) (MacAyeal, 1989),

**the downstream unit normal vector to the fluxgate contour and**

*n**ρ*

_{ice}the ice density. All 13 mass fluxes are computed simultaneously for each forward run of ISSM and each sample of the input quantities.

The anisotropic mesh used in this application comprises 2085 elements and 1112 nodes, with higher spatial resolution near the shear margins and the grounding line. The number of nodes is sufficiently small to make both the sampling (tenths of a second per sample) and the solving of the ice-sheet model (a few seconds per run) computationally fast.

## 4.2 Characterization of uncertainty

Following Larour et al. (2012b), we perturb each input quantity with a multiplicative noise proportional to a local (measurement) error margin. More specifically, we represent each source of uncertainty in the following way:

where *q*_{ref} is a reference value for each source of uncertainty, *σ*_{q} is the (spatially distributed) error margin for each source of
uncertainty and *p*_{q}(** s**) is a spatially distributed random perturbation. The reference thickness value is derived from bedrock topography
(Nitsche et al., 2007) and surface elevation (Bamber et al., 2009; Griggs and Bamber, 2009) data. The temperature-dependent reference ice hardness value is derived from
an Arrhenius law (Larour et al., 2012c) and surface temperature data by Comiso (2000). The reference basal drag coefficient is inferred using an
inverse method (Morlighem et al., 2010) with the reference thickness and ice hardness values. For the UQ analysis in Sect. 4.3
and 4.4, we only perturb the ice thickness, and the error margins are based on ground-penetrating radar (GPR) cross-over measurements
from the 2009 Operation IceBridge Campaign and provided by the Center for Remote Sensing of Ice Sheets (CReSIS); see also Larour et al. (2012b) for
further information. For the sensitivity analysis in Sect. 4.5, we specify a spatially uniform 5 % error margin for all three
inputs to compare sensitivities between input variables with same error margins.

Here, we represent the random perturbation *p*_{q} as a Gaussian random field with Matérn covariance function and unit variance (*σ*^{2}=1)
following the SPDE approach. Therefore, each uncertain input quantity is represented as a Gaussian random field with mean *q*_{ref}, Matérn
covariance function and marginal standard deviation proportional to *q*_{ref}(** s**)

*σ*

_{q}(

**). We solve the SPDE (Eq. 2) for**

*s**α*=2 and impose Robin boundary conditions with $\mathit{\beta}=\left(\mathit{\kappa}\mathit{\tau}\right)/\mathrm{1.42}$. We consider different parameter ranges

*ρ*to investigate the impact of the spatial correlation on the results. The parameters

*κ*and

*τ*are determined from Eqs. (3a)–(3c). As an illustration, we represent samples of the random perturbation for

*ρ*= 10 km in Fig. 6,

*ρ*= 30 km in Fig. 7 and

*ρ*= 100 km in Fig. 8. For this application, we have chosen the values of the parameter range arbitrarily without any proper data assimilation or inversion process. The use of an unconstrained parameter range is justified by the fact that we do not seek to provide new probabilistic mass balance estimates for the Pine Island Glacier but seek only to demonstrate the interest of the implemented stochastic sampler for forward UQ analysis.

## 4.3 Propagation of uncertainty

We use ISSM to propagate uncertainty in ice thickness into a probabilistic characterization of the depth-averaged mass flux across fluxgates. We
implement the propagation of uncertainty using Monte Carlo sampling. To this end, we begin by generating an ensemble of *n* i.i.d samples of the
random perturbation by solving the stochastic equation (Eq. 2) for an ensemble of *n* i.i.d samples of the random right-hand side
𝒲(** s**). Following this, we run the ISSM computational model for each sample of the input variable to generate an ensemble of

*n*i.i.d samples of the depth-averaged mass flux across the 13 fluxgates. From these samples, we can estimate the mean and standard deviation of the mass flux across the fluxgates with the statistical mean and standard deviation of the samples. In addition, we can estimate the probability density distribution of the mass flux across the fluxgates by using a normalized histogram of the samples or kernel density estimation (Scott, 2015). In all of our simulations, we found that

*n*=2500 samples were sufficient to ensure reasonably converged statistical estimates (bootstrap error of a few hundredths of a percent for the mean and a few percent for the standard deviation).

Figure 9 shows, for a parameter range of 30 km, the statistical distribution of mass flux through each gate as a normalized histogram of the samples, along with a kernel density estimate of their probability density function. The estimated mean and coefficient of variation (defined as the standard deviation over the mean) are also provided. The statistical distribution of mass flux through each gate can be considered Gaussian (which can be checked formally by carrying out a Kolmogorov–Smirnov test), with the exception of gate 4, for which the distribution is negatively skewed. The estimated mean values are approximately equal to the mass fluxes obtained for the reference ice thickness. The coefficients of variation range from of a few hundredths of percent to a few percent, with the fluxgates close to the ice front exhibiting the largest uncertainty. As already shown in Larour et al. (2012b), the mass flux response to input uncertainties is particularly smooth.

Figure 10 shows the statistical distribution of mass flux through gate 3 for different parameter ranges *ρ*. We also
present results for a white noise perturbation (zero correlation length) and a uniform Gaussian perturbation (infinite correlation length). We
observe that the mean value does not depend on the parameter range while the overall uncertainty (measured as the coefficient of variation) increases
with the parameter range. With the exception of gate 4, similar conclusions can be drawn for all the other gates. The specific behaviour of gate 4 may
result from the non-Gaussianity of the mass flux distribution, and the coefficient of variation is not the most appropriate measure of dispersion for
skewed distributions.

## 4.4 Sensitivity analysis to ice thickness changes

We carry out a global sensitivity analysis to determine where changes in ice thickness most influence uncertainty in mass flux across fluxgates. For a
given fluxgate, we can build a sensitivity map which gives for each location a sensitivity measure (or index) between the random field at this
location and the mass flux across the fluxgate. Here, we consider the coefficient of correlation between the random field and
the mass flux *M*_{i} as a sensitivity measure, i.e.

where Cov denotes the covariance between two random variables and 𝕍 denotes the variance. We mention that other sensitivity measures like the Spearman's or Kendall rank correlations, the distance correlation sensitivity index, or the Hilbert–Schmidt independence criterion sensitivity index (Da Veiga, 2014; De Lozzo and Marrel, 2016) can also be considered as long as these indices can handle dependent variables.

Figure 11 shows a sensitivity map for the mass flux through each gate for a parameter range of 30 km. The coefficients of correlation are estimated from the 2500 samples obtained in Sect. 4.3. A positive (negative) index value indicates a positive (negative) linear correlation between changes in ice thickness and changes in mass flux across a fluxgate, and a zero value indicates no linear correlation. As expected, only a limited region on the glacier has a significant influence on the mass flux across a given fluxgate. This region is usually concentrated close to the fluxgate, though more remote locations can have an influence, as can be seen for instance with fluxgates 6, 7 and 9. These remote locations lie in the main ice stream of the PIG, thus suggesting that uncertainty in ice thickness in the main ice stream can impact ice flow in upstream tributaries. With the exception of gate 4, sensitivity maps show a positive (linear) correlation between changes in ice thickness and changes in mass flux across fluxgates.

Figure 12 shows a sensitivity map of mass flux through gate 3 for different parameter ranges and white noise and uniform perturbations. As expected, increasing the parameter range expands the zone of influence and the magnitude of the correlation. Therefore, the larger the correlation length, the further away from the gate uncertainties in ice thickness can influence mass flux estimates.

## 4.5 Sensitivity analysis to multiple input quantities

We perform another sensitivity analysis by considering multiple sources of uncertainty, namely uncertainty in ice thickness, basal friction and ice
hardness. We assume that these three sources of uncertainty are statistically independent, and we choose a parameter range of 30 km for all of
them. Because the sources of uncertainty are statistically independent, we can use the Sobol or variance-based sensitivity indices
(Sobol', 2001; Saltelli et al., 2008) to characterize the contribution of each source of uncertainty to the uncertainty in the mass flux across each
fluxgate. The first-order Sobol index *S*_{i,j} measuring the contribution of the random field $\left\{{x}_{i}\left(\mathit{s}\right),\mathit{s}\in D\right\}$ to the
uncertainty in the mass flux *M*_{j} is defined as

The Sobol index of a given uncertain input parameter represents the fraction of the variance of the projections explained as stemming from this sole
uncertain parameter. A value of 1 indicates that the entire variance of the projections is explained by this sole uncertain parameter, and a value of 0
indicates that the uncertain parameter has no impact on the projection uncertainty. We estimate these Sobol indices using Monte Carlo sampling
following the estimation method by Janon et al. (2014). We use a total of *ν*=4 × 5000 evaluations of the computational model to estimate the
Sobol indices.

Figure 13 shows the Sobol sensitivity indices for each gate. Though the sum of all three sensitivity indices should be less than 1, the sum of the estimated indices exceeds this value for a few gates (gates 6, 12 and 13). This suggests that the Sobol estimates are not totally converged for some gates. From Fig. 13, we find that the uncertainty in the ice thickness has the largest influence on the uncertainty in the mass flux across fluxgates, accounting for more than 60 % to 80 % of their variance. The second most influential input parameter is the basal friction coefficient (Sobol indices ranging from 5 % to 35 %), except for gates 1 and 13 lying on the ice shelf where there is no basal friction.

## 4.6 Transient analysis

Finally, we present a transient experiment of the Pine Island Glacier over a 10-year period. The Pine Island Glacier is forced by basal melting
underneath the floating portion of the domain and surface mass balance. We impose a basal melting rate of 25 m yr^{−1}, and we consider the
surface mass balance (smb) to be uncertain. We represent the uncertain smb as

where smb_{ref} is a reference value for the smb taken from Vaughan et al. (1999), *σ*_{smb} is a constant 10 % error
margin and *p*_{smb}(** s**,

*t*) a spatio-temporal random perturbation. The spatio-temporal random perturbation is represented as a spatio-temporal Gaussian random field given by Eq. (16). The spatial noise has a parameter range of 30 km. We assume that the smb exhibits a temporal monthly variability, and consequently we discretize the temporal process in time with a 1-month time step. We use the same time step to run the ISSM model forward in time. Figures 14 and 15 show two realizations of the random perturbation for a temporal correlation of 0.5 and 0.95, respectively.

Figure 16 shows the statistical distribution of mass flux through gate 3 as a function of time for a temporal correlation of 0.5. We observe an increase in the mean mass flux through the gate due mainly to the ocean forcing, with the estimated mean values approximately equal to the mass flux in the absence of stochastic smb forcing. Most gates show an increase in their mean mass flux in time, except for gate 4, which shows a decrease in its mean mass flux, and gates 2, 8, 9 and 13, which show a relatively constant mean mass flux in time. However, we observe an increase in the coefficient of variation of the mass flux for all gates, thus suggesting that mass balance estimates become increasingly uncertain in time.

Figure 17 shows the coefficient of variation at the end of the simulation as a function of the
temporal correlation *ϕ* for a few gates. From this figure, it is clear that increased temporal correlation results in increased uncertainty in the mass flux, with
a more significant increase for higher values of the temporal correlation.

Model-based projections of ice-sheet mass balance should all be ideally given in the form of probabilistic projections that provide an assessment of
the impact of uncertainties in boundary conditions, climate forcing, ice sheet geometry or initial conditions. In this context, uncertainty
quantification analysis with ice-sheet models should capture the spatial and temporal variability of these sources of uncertainties. To this end,
random fields are an essential statistical tool for the modelling and analysis of uncertain spatio-temporal processes. Here, we take advantage of the
explicit link between Gaussian random fields with Matérn covariance function and a stochastic partial differential equation to implement a random
field sampler within ISSM. The FEM-oriented implementation of ISSM and its fully object-oriented and highly parallelized architecture allow for a
straightforward and computationally efficient implementation of this SPDE. This stochastic sampler provides an alternative to the pre-existing UQ mesh
partitioning approach implemented within the ISSM-DAKOTA UQ framework. The partitioning approach can be interpreted as the representation of a random
field with perfect correlation for locations in a same partition and zero (Larour et al., 2012b) or non-zero correlation (Larour et al., 2020) for locations
in different partitions. Because the partitions are constant throughout the sampling, realizations strongly depend on the partitioning, which is not
the case with the stochastic sampler we implemented. In the partitioning approach, the number of partitions artificially controls the spatial scale of
the random perturbation. In contrast, spatial correlation is an intrinsic feature of Gaussian Matérn random fields and can be controlled through
the parameter range *ρ*. In addition, random samples generated following the SPDE approach are mesh independent as long as the mesh resolution is
sufficiently fine to capture the desired spatial variability.

The results presented here provide complementary insight into the impact of uncertainties on the PIG. First, our results are consistent with previous results by Larour et al. (2012b) showing that the SSA model is stable to uncertainties in ice thickness. For most gates, the probability distribution of the mass flux can be considered as a Gaussian distribution centred around the reference value. The overall uncertainty is relatively small and does not exceed a few percent of the mean value at most. This suggests that the response of the PIG to changes in ice thickness is relatively smooth even if the ice flow dynamics is non-linear. Therefore, robust projections of mass balance estimates of the PIG are possible even in the presence of uncertainties in ice thickness. Second, our results show that a higher spatial correlation length results in larger uncertainty in mass balance estimates. The correlation length controls the size of the region around each gate in which changes in ice thickness significantly impact mass flux through the gates (see Fig. 12). This conclusion is consistent with Schlegel et al. (2018), who showed that many small partitions would likely underestimate uncertainty in model outputs while the use of a single partition would likely overestimate the uncertainty. Third, our sensitivity analysis shows that ice thickness is the most influential parameter in inducing uncertainty in mass flux estimates. Although this conclusion is consistent with Larour et al. (2012b), both conclusions cannot be directly compared because we carried out a global sensitivity analysis rather than the local sensitivity analysis performed in the latter study. A local sensitivity analysis measures the sensitivity of the mass flux to local changes in input data and assumes that each input variable is independent. This hypothesis is not appropriate for correlated random fields. In contrast, the global sensitivity analysis in Sect. 4.5 allows us to consider the global uncertainty in the random field and only requires the different random fields to be statistically independent. Fourth, our transient simulation shows that uncertainty in mass balance estimates increases in time and for higher temporal correlations.

Finally, we discuss a few limitations of the SPDE approach. First, this approach relies on the assumption that uncertainty in spatially varying input
parameters can be appropriately characterized by a Gaussian random field. Though Gaussian random fields are a practical model for many stochastic
phenomena, this assumption may not be appropriate because input parameters can have an asymmetric distribution or heavy tails or can simply be positive or
bounded. While this may sound restrictive, non-Gaussian random fields can be obtained through non-linear transformations of Gaussian random fields
(Vio et al., 2001; Hristopulos, 2020) or by solving the SPDE with a non-Gaussian white noise (Bolin, 2014). Second, representing spatio-temporal
processes with a separable autoregressive model relies on the assumption that the spatial and temporal variabilities are uncoupled and that the
temporal variability can be simply described with a first-order Markov process. These two assumptions may be too restrictive to describe the complex
spatio-temporal variability of stochastic processes like the surface mass balance (White et al., 2019). Third, the SPDE and autoregressive models
require the specification of model parameters, most notably the spatial and temporal correlation lengths. As shown in Sect. 4,
these parameters have a major impact on the uncertainty in mass balance estimates. Therefore, these parameters should ideally be constrained through
data assimilation methods, modelling efforts, and inverse methods, including Bayesian inference; see, for instance Cameletti et al. (2012) and
Western et al. (2020), for Bayesian spatio-temporal inference with hidden Gaussian random fields. Lastly, we assumed the spatial correlation to be
constant on the PIG. While this assumption might seem reasonable for a single glacier, this may be more questionable for an ice sheet with different
drainage basins. However, this issue may be handled by considering a spatially varying parameter *κ*, as already discussed in
Sect. 2.1.1.

We presented and implemented a Gaussian random field sampler within ISSM for uncertainty quantification analysis of spatially (and temporally) varying uncertain input parameters in ice-sheet models. So far, sampling of spatially varying uncertain input parameters within ISSM has relied on a partitioning of the computational domain and sampling from the partitions. This approach may be limited in representing spatial and temporal correlations. To improve the probabilistic characterization of spatially (and temporally) correlated uncertain input parameters, we proposed to generate realizations of Gaussian random fields with Matérn covariance function. Our implementation relies on an explicit link between Gaussian random fields with Matérn covariance function and a stochastic partial differential equation. This SPDE sampling approach allows for a computationally efficient sampling of random fields using the FEM and provides a flexible mathematical framework to build probabilistic models of spatio-temporal uncertain processes. In addition, model parameters allow us to intrinsically control the spatial and temporal correlations of the realizations and their variance and smoothness.

We applied this stochastic sampler to investigate the impact of spatially (and temporally) varying sources of uncertainties, namely uncertainties in ice thickness, basal friction, ice hardness and surface mass balance, on the Pine Island Glacier. More specifically, we investigated the impact of spatial and temporal correlations on ice-sheet mass balance. We showed that the SSA model is stable to uncertainties. We also showed that larger correlation lengths lead to increased uncertainty in mass balance estimates because the random perturbations have a larger length scale of influence. We found that the most influential source of uncertainties for estimating mass balance is the uncertainty in ice thickness. We also showed in a transient experiment that uncertainty in mass balance estimates increases in time and with higher temporal correlation. Overall, our results demonstrate the need to better constrain the spatial and temporal variability of physical processes impacting ice-sheet dynamics through data assimilation and modelling efforts.

The ISSM code can be downloaded, compiled and executed following the instructions available on the ISSM website https://issm.jpl.nasa.gov/download (last access: 7 September 2021). The public SVN repository for the ISSM code can also be found directly at https://issm.ess.uci.edu/svn/issm/issm/trunk (last access: 17 September 2021) and downloaded using username “anon” and password “anon”. The archived version of the source code used in this paper is made available as part of a Zenodo repository at https://doi.org/10.5281/zenodo.5532775 (Bulthuis and Larour, 2021a). Data to reproduce the results in Sects. 3 and 4 are made available as part of a Zenodo repository at https://doi.org/10.5281/zenodo.5532710 (Bulthuis and Larour, 2021b).

KB and EL discussed the results presented in this paper. KB implemented the stochastic sampler into ISSM and carried out the simulations. KB wrote the bulk of the manuscript with relevant comments from EL.

The contact author has declared that neither they nor their co-author has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We would like to thank the two anonymous referees for their very helpful comments that helped improve the overall quality and readability of the manuscript. Kevin Bulthuis acknowledges support from the NASA Postdoctoral Program (NPP), administered by the Universities Space Research Association (USRA) under contract with the National Aeronautics and Space Administration (NASA). The research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Government sponsorship is acknowledged.

This research has been supported by the Jet Propulsion Laboratory (NASA Postdoctoral Program).

This paper was edited by Alexander Robel and reviewed by two anonymous referees.

Abramowitz, M. and Stegun, I. A.: Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables, 9th edn., Dover Publications, New York, NY, 1970. a

Babaniyi, O., Nicholson, R., Villa, U., and Petra, N.: Inferring the basal sliding coefficient field for the Stokes ice sheet model under rheological uncertainty, The Cryosphere, 15, 1731–1750, https://doi.org/10.5194/tc-15-1731-2021, 2021. a, b

Bakka, H., Krainski, E., Bolin, D., Rue, H., and Lindgren, F.: The diffusion-based extension of the Matérn field to space-time, arXiv [preprint], arXiv:2006.04917, 2020. a

Bamber, J. L., Gomez-Dans, J. L., and Griggs, J. A.: A new 1 km digital elevation model of the Antarctic derived from combined satellite radar and laser data – Part 1: Data and methods, The Cryosphere, 3, 101–111, https://doi.org/10.5194/tc-3-101-2009, 2009. a

Beskos, A., Girolami, M., Lan, S., Farrell, P. E., and Stuart, A. M.: Geometric MCMC for infinite-dimensional inverse problems, J. Comput. Phys., 335, 327–351, https://doi.org/10.1016/j.jcp.2016.12.041, 2017. a

Bolin, D.: Spatial Matérn Fields Driven by Non-Gaussian Noise, Scand. J. Stat., 41, 557–579, https://doi.org/10.1111/sjos.12046, 2014. a, b

Bolin, D. and Lindgren, F.: Wavelet Markov models as efficient alternatives to tapering and convolution fields, Tech. rep., Centre for Mathematical Sciences, Faculty of Engineering, Lund University, 2009. a

Bolin, D., Kirchner, K., and Kovács, M.: Numerical solution of fractional elliptic stochastic PDEs with spatial white noise, IMA J. Numer. Anal., 40, 1051–1073, https://doi.org/10.1093/imanum/dry091, 2018. a

Brinkerhoff, D. J.: Variational Inference at Glacier Scale, arXiv [preprint], arXiv:2108.07263, 2021. a, b

Bui-Thanh, T., Ghattas, O., Martin, J., and Stadler, G.: A Computational Framework for Infinite-Dimensional Bayesian Inverse Problems Part I: The Linearized Case, with Application to Global Seismic Inversion, SIAM J. Sci. Comput., 35, A2494–A2523, https://doi.org/10.1137/12089586x, 2013. a, b, c

Bulthuis, K. and Larour, E.: A new sampling capability for uncertainty quantification in the Ice-sheet and Sea-level System Model v4.19 using Gaussian Markov random fields – Software, Zenodo [code], https://doi.org/10.5281/zenodo.5532775, 2021a. a

Bulthuis, K. and Larour, E.: A new sampling capability for uncertainty quantification in the Ice-sheet and Sea-level System Model v4.19 using Gaussian Markov random fields – Datasets and results, Zenodo [data set], https://doi.org/10.5281/zenodo.5532710, 2021b. 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, https://doi.org/10.5194/tc-13-1349-2019, 2019. a, b

Bulthuis, K., Pattyn, F., and Arnst, M.: A Multifidelity Quantile-Based Approach for Confidence Sets of Random Excursion Sets with Application to Ice-Sheet Dynamics, SIAM/ASA Journal on Uncertainty Quantification, 8, 860–890, https://doi.org/10.1137/19m1280466, 2020. a

Cameletti, M., Lindgren, F., Simpson, D., and Rue, H.: Spatio-temporal modeling of particulate matter concentration through the SPDE approach, ASTA-Adv. Stat. Anal., 97, 109–131, https://doi.org/10.1007/s10182-012-0196-3, 2012. a, b, c

Chu, S. and Guilleminot, J.: Stochastic multiscale modeling with random fields of material properties defined on nonconvex domains, Mech. Res. Commun., 97, 39–45, https://doi.org/10.1016/j.mechrescom.2019.01.008, 2019. a

Comiso, J. C.: Variability and Trends in Antarctic Surface Temperatures from In Situ and Satellite Infrared Measurements, J. Climate, 13, 1674–1696, https://doi.org/10.1175/1520-0442(2000)013<1674:vatias>2.0.co;2, 2000. a

Cressie, N. A. C.: Statistics for Spatial Data, revised edn., John Wiley & Sons, Inc., New York, NY, https://doi.org/10.1002/9781119115151, 1993. a

Da Veiga, S.: Global sensitivity analysis with dependence measures, J. Stat. Comput. Sim., 85, 1283–1305, https://doi.org/10.1080/00949655.2014.945932, 2014. a

Daon, Y. and Stadler, G.: Mitigating the influence of the boundary on PDE-based covariance operators, Inverse Probl. Imag., 12, 1083–1102, https://doi.org/10.3934/ipi.2018045, 2018. a, b

De Lozzo, M. and Marrel, A.: New improvements in the use of dependence measures for sensitivity analysis and screening, J. Stat. Comput. Sim., 86, 3038–3058, https://doi.org/10.1080/00949655.2016.1149854, 2016. a

Edwards, T. L., Brandon, M. A., Durand, G., Edwards, N. R., Golledge, N. R., Holden, P. B., Nias, I. J., Payne, A. J., Ritz, C., and Wernecke, A.: Revisiting Antarctic ice loss due to marine ice-cliff instability, Nature, 566, 58–64, https://doi.org/10.1038/s41586-019-0901-4, 2019. a, b

Edwards, T. L., Nowicki, S., Marzeion, B., Hock, R., Goelzer, H., Seroussi, H., Jourdain, N. C., Slater, D. A., Turner, F. E., Smith, C. J., McKenna, C. M., Simon, E., Abe-Ouchi, A., Gregory, J. M., Larour, E., Lipscomb, W. H., Payne, A. J., Shepherd, A., Agosta, C., Alexander, P., Albrecht, T., Anderson, B., Asay-Davis, X., Aschwanden, A., Barthel, A., Bliss, A., Calov, R., Chambers, C., Champollion, N., Choi, Y., Cullather, R., Cuzzone, J., Dumas, C., Felikson, D., Fettweis, X., Fujita, K., Galton-Fenzi, B. K., Gladstone, R., Golledge, N. R., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huss, M., Huybrechts, P., Immerzeel, W., Kleiner, T., Kraaijenbrink, P., Le clec'h, S., Lee, V., Leguy, G. R., Little, C. M., Lowry, D. P., Malles, J.-H., Martin, D. F., Maussion, F., Morlighem, M., O'Neill, J. F., Nias, I., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Radić, V., Reese, R., Rounce, D. R., Rückamp, M., Sakai, A., Shafer, C., Schlegel, N.-J., Shannon, S., Smith, R. S., Straneo, F., Sun, S., Tarasov, L., Trusel, L. D., Breedam, J. V., van de Wal, R., van den Broeke, M., Winkelmann, R., Zekollari, H., Zhao, C., Zhang, T., and Zwinger, T.: Projected land ice contributions to twenty-first-century sea level rise, Nature, 593, 74–82, https://doi.org/10.1038/s41586-021-03302-y, 2021. a, b

Furrer, R., Genton, M. G., and Nychka, D.: Covariance Tapering for Interpolation of Large Spatial Datasets, J. Comput. Graph. Stat., 15, 502–523, https://doi.org/10.1198/106186006x132178, 2006. a

Griggs, J. A. and Bamber, J. L.: A new 1 km digital elevation model of Antarctica derived from combined radar and laser data – Part 2: Validation and error estimates, The Cryosphere, 3, 113–123, https://doi.org/10.5194/tc-3-113-2009, 2009. a

Hendrickson, B. and Leland, R.: The Chaco user's guide, version 2.0, Technical Report SAND-95-2344, Tech. rep., Sandia National Laboratories, https://cfwebprod.sandia.gov/cfdocs/CompResearch/docs/guide.pdf (last access: 3 September 2021), 1995. a

Hristopulos, D. T.: Random Fields for Spatial Data Modeling: A Primer for Scientists and Engineers, Advances in Geographic Information Science, Springer Netherlands, Dordrecht, the Netherlands, https://doi.org/10.1007/978-94-024-1918-4, 2020. a, b

Isaac, T., Petra, N., Stadler, G., and Ghattas, O.: Scalable and efficient algorithms for the propagation of uncertainty from data through inference to prediction for large-scale problems, with application to flow of the Antarctic ice sheet, J. Comput. Phys., 296, 348–368, https://doi.org/10.1016/j.jcp.2015.04.047, 2015. a, b, c, d, e, f

Janon, A., Klein, T., Lagnoux, A., Nodet, M., and Prieur, C.: Asymptotic normality and efficiency of two Sobol index estimators, ESAIM-Probab. Stat., 18, 342–364, https://doi.org/10.1051/ps/2013040, 2014. a

Khristenko, U., Scarabosio, L., Swierczynski, P., Ullmann, E., and Wohlmuth, B.: Analysis of Boundary Effects on PDE-Based Sampling of Whittle–Matérn Random Fields, SIAM/ASA Journal on Uncertainty Quantification, 7, 948–974, https://doi.org/10.1137/18m1215700, 2019. a, b

Krainski, E. T.: Statistical Analysis of Space-time Date: New Models and Applications, PhD thesis, Norwegian University of Science and Technology, Trondheim, Norway, 2018. a

Larour, E., Morlighem, M., Seroussi, H., Schiermeier, J., and Rignot, E.: Ice flow sensitivity to geothermal heat flux of Pine Island Glacier, Antarctica, J. Geophys. Res.-Earth, 117, F04023, https://doi.org/10.1029/2012jf002371, 2012a. a

Larour, E., Schiermeier, J., Rignot, E., Seroussi, H., Morlighem, M., and Paden, J.: Sensitivity Analysis of Pine Island Glacier ice flow using ISSM and DAKOTA, J. Geophys. Res.-Earth, 117, F02009, https://doi.org/10.1029/2011jf002146, 2012b. a, b, c, d, e, f, g, h, i, j

Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res.-Earth, 117, F01022, https://doi.org/10.1029/2011jf002140, 2012c. a, b, c

Larour, E., Caron, L., Morlighem, M., Adhikari, S., Frederikse, T., Schlegel, N.-J., Ivins, E., Hamlington, B., Kopp, R., and Nowicki, S.: ISSM-SLPS: geodetically compliant Sea-Level Projection System for the Ice-sheet and Sea-level System Model v4.17, Geosci. Model Dev., 13, 4925–4941, https://doi.org/10.5194/gmd-13-4925-2020, 2020. a

Lindgren, F., Rue, H., and Lindström, J.: An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach, J. R. Stat. Soc. B, 73, 423–498, https://doi.org/10.1111/j.1467-9868.2011.00777.x, 2011. a, b, c, d, e, f, g, h, i, j

MacAyeal, D. R.: Large-scale ice flow over a viscous basal sediment: Theory and application to ice stream B, Antarctica, J. Geophys. Res.-Sol. Ea., 94, 4071–4087, https://doi.org/10.1029/jb094ib04p04071, 1989. a

Minasny, B. and McBratney, A. B.: The Matérn function as a general model for soil variograms, Geoderma, 128, 192–207, https://doi.org/10.1016/j.geoderma.2005.04.003, 2005. 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, L14502, https://doi.org/10.1029/2010gl043853, 2010. a

Nitsche, F. O., Jacobs, S. S., Larter, R. D., and Gohl, K.: Bathymetry of the Amundsen Sea continental shelf: Implications for geology, oceanography, and glaciology, Geochem. Geophy. Geosy., 8, Q10009, https://doi.org/10.1029/2007gc001694, 2007. a

Petra, N., Martin, J., Stadler, G., and Ghattas, O.: A Computational Framework for Infinite-Dimensional Bayesian Inverse Problems, Part II: Stochastic Newton MCMC with Application to Ice Sheet Flow Inverse Problems, SIAM J. Sci. Comput., 36, A1525–A1555, https://doi.org/10.1137/130934805, 2014. a, b, c, d

Rasmussen, C. E. and Williams, C. K. I.: Gaussian Processes for Machine Learning, Adaptive Computation and Machine Learning, The MIT Press, Cambridge, MA, 2006. a

Rignot, E.: Changes in West Antarctic ice stream dynamics observed with ALOS PALSAR data, Geophys. Res. Lett., 35, L12505, https://doi.org/10.1029/2008GL033365, 2008. a

Ritz, C., Edwards, T. L., Durand, G., Payne, A. J., Peyaud, V., and Hindmarsh, R. C. A.: Potential sea-level rise from Antarctic ice-sheet instability constrained by observations, Nature, 528, 115–118, https://doi.org/10.1038/nature16147, 2015. a, b

Roininen, L., Huttunen, J. M. J., Lasanen, S., and and: Whittle–Matérn priors for Bayesian statistical inversion with applications in electrical impedance tomography, Inverse Probl. Imag., 8, 561–586, https://doi.org/10.3934/ipi.2014.8.561, 2014. a

Rue, H. and Held, L.: Gaussian Markov random fields: Theory and Applications, vol. 104 of Monographs on statistics and applied probability, Chapman & Hall/CRC, Boca Raton, FL, 2005. a

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., and Tarantola, S.: Global Sensitivity Analysis: The Primer, John Wiley & Sons, West Sussex, UK, 2008. a

Schabenger, O. and Gotway, C. A.: Statistical Methods for Spatial Data Analysis, Chapman & Hall, Boca Raton, FL, ISBN 1-58488-322-7, 2005. a

Schlegel, N.-J. and Larour, E. Y.: Quantification of Surface Forcing Requirements for a Greenland Ice Sheet Model Using Uncertainty Analyses, Geophys. Res. Lett., 46, 9700–9709, https://doi.org/10.1029/2019gl083532, 2019. a

Schlegel, N.-J., Larour, E., Seroussi, H., Morlighem, M., and Box, J. E.: Decadal-scale sensitivity of Northeast Greenland ice flow to errors in surface mass balance using ISSM, J. Geophys. Res.-Earth, 118, 667–680, https://doi.org/10.1002/jgrf.20062, 2013. a

Schlegel, N.-J., Larour, E., Seroussi, H., Morlighem, M., and Box, J. E.: Ice discharge uncertainties in Northeast Greenland from boundary conditions and climate forcing of an ice flow model, J. Geophys. Res.-Earth, 120, 29–54, https://doi.org/10.1002/2014jf003359, 2015. 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, https://doi.org/10.5194/tc-12-3511-2018, 2018. a, b, c

Scott, D. W.: Multivariate Density Estimation: Theory, Practice, and Visualization, 2nd edn., Wiley Series in Probability and Statistics, Wiley, Hoboken, NJ, https://doi.org/10.1002/9781118575574, 2015. a

Simpson, D., Lindgren, F., and Rue, H.: In order to make spatial statistics computationally feasible, we need to forget about the covariance function, Environmetrics, 23, 65–74, https://doi.org/10.1002/env.1137, 2011. a

Sobol', I.: Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Math. Comput. Simulat., 55, 271–280, https://doi.org/10.1016/s0378-4754(00)00270-6, 2001. a

Stein, M. L.: Interpolation of Spatial Data: Some theory for Kriging, Springer Series in Statistics, Springer, New York, NY, https://doi.org/10.1007/978-1-4612-1494-6, 1999. a, b

Stuart, A. M.: Inverse problems: A Bayesian perspective, Acta Numer., 19, 451–559, https://doi.org/10.1017/s0962492910000061, 2010. a

Vaughan, D. G., Bamber, J. L., Giovinetto, M., Russell, J., and Cooper, A. P. R.: Reassessment of Net Surfcace Mass Balance in Antarctica, J. Climate, 12, 933–946, https://doi.org/10.1175/1520-0442(1999)012<0933:ronsmb>2.0.co;2, 1999. a

Vio, R., Andreani, P., and Wamsteker, W.: Numerical Simulation of Non-Gaussian Random Fields with Prescribed Correlation Structure, Publ. Astron. Soc. Pac., 113, 1009–1020, https://doi.org/10.1086/322919, 2001. a

Western, L. M., Sha, Z., Rigby, M., Ganesan, A. L., Manning, A. J., Stanley, K. M., O'Doherty, S. J., Young, D., and Rougier, J.: Bayesian spatio-temporal inference of trace gas emissions using an integrated nested Laplacian approximation and Gaussian Markov random fields, Geosci. Model Dev., 13, 2095–2107, https://doi.org/10.5194/gmd-13-2095-2020, 2020. a, b, c

White, P. A., Reese, C. S., Christensen, W. F., and Rupper, S.: A model for Antarctic surface mass balance and ice core site selection, Environmetrics, 30, e2579, https://doi.org/10.1002/env.2579, 2019. a

Whittle, P.: On Stationary Processes in the Plane, Biometrika, 41, 434, https://doi.org/10.2307/2332724, 1954. a