A new sampling capability for uncertainty quantification in the Ice-sheet and Sea-level System Model v4.19 using Gaussian Markov random fields

Abstract. 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 involves first the (probabilistic) characterization of the sources of uncertainty followed by the propagation and sensitivity analysis of these sources of uncertainty. Previous studies 5 concerned with UQ in ice-sheet models have generally focused on the last two steps but 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 in ice-sheet models, we propose in this paper to rep10 resent them as 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. Samples from a Gaussian random field with Matérn covariance function can be generated efficiently by solving a certain stochastic partial differential equation. Discretization of this stochastic partial differential equation by the finite element method results in a sparse approximation known as a Gaussian Markov random field. We solve this equation efficiently 15 using the finite element method within the Ice-sheet and Sea-level System Model (ISSM). 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 well the desired spatial and temporal correlations. Finally, we demonstrate the interest of this sampling capability in an illustration concerned with assessing the impact of various sources of uncertainties on the Pine Island Glacier, West Antarctica. We find that both larger spatial and temporal correlations lengths will likely result 20 in increased uncertainty in the projections.


Introduction
Despite large improvements in ice-sheet modeling 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 modeled 25 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. 30 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(Bulthuis et al., , 2020Edwards et al., 2019Edwards et al., , 2021Larour 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 35 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 unmodeled dynamics uncertainty (Bulthuis et al., 2019;Edwards et al., 2019Edwards et al., , 2021Ritz et al., 2015). These uncertain lumped parameters then appear in parameterizations and reduced-order models of complex physical processes such as basal friction or oceaninduced melting. However, such a characterization does not reflect the fact that most sources of uncertainty in numerical 40 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. (2013Schlegel et al. ( , 2015Schlegel et al. ( , 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, function (Rasmussen and Williams, 2006). In practice, the covariance function is often assumed to be a function of the space lag s − s (stationary covariance function) or even the spatial distance s − s (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 modeling 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 90 (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, Γ the Gamma function, σ 2 the variance of the random field, and κ a positive scaling parameter. The scaling parameter κ is more naturally 95 interpreted as a range parameter ρ, with the empirically derived relationship ρ = √ 8ν/κ corresponding to correlations near 0.1 at the distance ρ . 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 modeling, 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 100 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. Various methods exist to sample Gaussian random fields, see for instance Hristopulos (2020); yet these methods can be 105 computationally prohibitive for applications in glaciology. For instance, methods based on a factorization of the covariance matrix typically results 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 110 stochastic partial differential equation where κ and τ are positive parameters, α = ν + d/2 (where d is the spatial dimension), and W(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, for instance, Khristenko et al. (2019)): An intuitive interpretation for the rationale behind the SPDE (3) 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 120 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 (3) requires a choice of boundary conditions. Here, we consider either homogeneous Neumann boundary conditions 125 ∇x(s) · n = 0, or Robin boundary conditions βx(s) + ∇x(s) · n = 0, 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 130 the boundary; see, for instance, Daon et al. (2018) and Khristenko et al. (2019) for a discusion and approaches to reduce the impact of boundary effects. For α = 2 and Robin boundary conditions, the choice of β = (κτ )/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. (3) can be solved efficiently by 135 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 to impose the covariance structure explicitly; see Simpson et al. (2011) for a discussion. For instance, the parameters κ and τ in Eq.
(3) can vary spatially in order to represent non-stationary Gaussian random fields with spatially varying range parameter and marginal variance. Similarly, the Robin coefficient can vary spatially to mitigate boundary effects (Daon et al., 2018).

140
Other extensions of the SPDE approach include anisotropic Gaussian models  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 .

Numerical implementation
We assume at this stage that α = 2. The SPDE (3) can be solved efficiently by the finite element method; see, for instance, 145 Lindgren et al. (2011) and Chu and Guilleminot (2019). The finite element representation of the solution writes as where {ψ k } N k=1 is a finite element basis (here piecewise linear functions), N the total number of nodes, and {x k } N k=1 are the stochastic nodal values. The discretized weak form of the differential operator in Eq. (3) 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 (3) is a Gaussian random 155 vector that satisfies where the covariance matrix is given by where the superscript underpins the order 2. Let Σ (2) = L (2) L (2) T be a matrix factorization of Σ (2) (e.g., Cholesky or QR 160 decomposition). Then, 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 x can obtained by solving the following linear system of equations

165
More generally, the covariance matrix Σ (α) for any arbitrary integer order α is given by ) where Σ (1) = K −1 and Σ (2) = K −1 MK −1 . Also the matrix L (α) for any arbitrary integer order α is given by where L (1) = K 1/2 and L (2) = K −1 M 1/2 . Therefore, samples for any arbitrary integer order α > 2 can be generated in a re-170 cursive 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 (3) to non-integer orders.
Algorithm 1 SPDE-based sampling algorithm for Gaussian random fields with Matérn covariance function Set input data α, τ , κ and β; Create FEM discretization: Build stiffness matrix K and mass matrix M; Create random vector ξ with independent entries ∼ N (0, 1); if α%2 = 0 then Compute the square root of K; return Sample x (α) ; In order to speed up the computation, the mass matrix M can be approximated as a diagonal lump mass matrix 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 ; elements in the precision matrix are non zero only for diagonal elements and neighboring 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 180 other nodes only depends on the neighboring nodes.

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

185
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 190 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, obtained as the solution of the SPDE (3). If x 0 is a stationary Gaussian random field with variance σ 2 and t is a stationary Gaussian random field with variance σ 2 = σ 2 (1 − φ 2 ), then x t is a stationary process in time with zero mean and marginal variance σ 2 if |φ| < 1. The process is non stationary if |φ| ≥ 1, with the case φ = 1 corresponding to a discrete random walk. The temporal 195 correlation function, or autocorrelation function, is given by and is therefore stationary 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. (3) and (17) defines a spatio-temporal Gaussian random field discretized in time.

200
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 beyond the scope of this work, non-separable models can however be formulated based on the SPDE approach (Krainski, 2018;Bakka et al., 2020).

Sampling capability: Verification
As a means of verifying and testing the implementation of our new sampling capability, we carry out a set of numerical experiments on a synthetic ice sheet/ice shelf of size 100 km × 100 km. 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 (3) for the orders α = 2 and 4, then discuss the impact of the boundary conditions on the simulations in a second 210 experiment. We finally perform transient simulations to verify the temporal covariance structure.

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

Experiment 2: Influence of the boundary conditions
The second experiment investigates the influence of the boundary conditions on the marginal variance of the samples gener-225 ated following the SPDE approach. Figure  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 245 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 thirteen fluxgates positioned along the main tributaries of the PIG (see Figure 5). The mass flux M i through gate i is given by where the integral is a line integral along the fluxgate contour C i , v the depth-averaged horizontal velocity computed using the Shallow Shelf-Stream Approximation (SSA) (MacAyeal, 1989), n the downstream unit normal vector to the fluxgate contour, and ρ ice the ice density. All thirteen mass fluxes are computed simultaneously for each forward run of ISSM and each sample of the input quantities.

Characterization of uncertainty
We represent each source of uncertainty in the following way where q ref is a reference value for each source of uncertainty, σ q 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 260 topography (Nitsche et al., 2007) and surface elevation Griggs and Bamber, 2009) data. The temperaturedependent 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 ( For the sensitivity analysis in Sect. 4.5, we specify a spatially uniform 5% error margin on 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 we represent samples of the random perturbation for ρ = 10 km in Fig. 6, ρ = 30 km in Fig. 7, and ρ = 100 km in Fig. 8. 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 our simulations, we found that n = 2500 samples were sufficient to ensure reasonably converged statistical estimates.  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

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 300 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 as sensitivity measure the coefficient of correlation between the random field and the mass flux M i , i.e., where Cov denotes the covariance between two random variables and V the variance. We mention that other sensitivity mea-     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 320 flux estimates.

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 325 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 {x i (s), s ∈ D} 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 330 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 335 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 on 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.
where smb ref is a reference value for the smb taken from Vaughan et al. (1999), σ smb 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. (17). 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 one-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 350 temporal correlation of 0.5 and 0.95, respectively.

Discussion
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 modeling and 365 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 new sampling capability 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 new sampling capability provides an alternative to the pre-existing UQ mesh partitioning approach implemented within the ISSM-DAKOTA UQ framework. The 370 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 new sampling capability 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 375 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 centered around the 380 reference value. The overall uncertainty is relatively small and does not exceed at most a few percents of the mean value.
This suggests that the response of the PIG to changes in ice thickness is relatively smooth even if the ice flow dynamics is nonlinear. 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 385 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 a local sensitivity analysis as performed in 390 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 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.

395
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, heavy tails, or simply be positive or bounded. While this may sound restrictive, non-Gaussian random fields can be obtained through nonlinear transformations of Gaussian random fields (Vio et al., 2001;400 Hristopulos, 2020) or by solving the SPDE with a non-Gaussian white noise (Bolin, 2014). Second, representing spatiotemporal 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 405 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, modeling 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 410 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 new sampling capability 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 in-415 put 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 420 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 to control intrinsically the spatial and temporal correlations of the realizations as well as their variance and smoothness.
We applied this new capability 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 425 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 430 the need to better constrain the spatial and temporal variability of physical processes impacting ice-sheet dynamics through data assimilation and modeling efforts.
Code and data availability. The ISSM code can be downloaded, compiled, and executed following the instructions available on the ISSM