2.1. Ice flow momentum balance forward model

Given a spatial domain with spatial coordinates specified by x, where for two dimensions x = [x,  y], our goal is to estimate the most likely spatial field of ice rigidity, B = B(x), conditional on observations of the flow of ice shelves and their geometry. To that end, we utilize a momentum balance method to estimate B that computes resistive stresses that optimally balance gravitational driving stresses. Within ice shelves, resistive (vertical) shear stresses at the base are negligible due to contact with relatively inviscid seawater. Ice is a thin film with small thicknesses relative to the horizontal dimensions (aerial extent). Thus, we are justified in employing the widely used shallow-shelf approximation (SSA), which assumes negligible vertical shearing in a thin film and vertically integrates viscosity and stresses in the ice column to obtain a simplified 2D framework for the governing equations of flow in ice shelves. We write the SSA momentum balance as:

(2)$$r_x = {\partial\over \partial x}\left(2 \eta h \left(2 {\partial u\over \partial x} + {\partial v\over \partial y}\right)\right) + {\partial\over \partial y}\left(\eta h \left({\partial u\over \partial y} + {\partial v\over \partial x}\right)\right)- \rho_{\rm i} g h {\partial s\over \partial x},\; $$

(3)$$r_y = {\partial\over \partial y}\left(2 \eta h \left(2 {\partial v\over \partial y} + {\partial u\over \partial x}\right)\right) + {\partial\over \partial x}\left(\eta h \left({\partial u\over \partial y} + {\partial v\over \partial x}\right)\right)- \rho_{\rm i} g h {\partial s\over \partial y},\; $$

where u and v are the horizontal velocity components of the velocity vector, u, along the x- and y-directions, respectively, and taken to be constant with depth; h is the ice thickness; s is the ice surface elevation; $2\eta = B \dot {\epsilon }_{\rm e}^{{1-n}/{n}}$ is the effective dynamic viscosity of ice (Eqn 1); ρ _i is the mass density of ice (assumed to be constant at ρ _i = 917 kg/m³); and g is the gravitational acceleration. In the above formulation, r _x and r _y on the left-hand side represent residual terms in the momentum balance. For the general SSA as applied to flowing ice where ice is grounded, such as at ice streams, these residual terms are non-zero and correspond to basal drag. However, since drag at the base of ice shelves is assumed to be negligible because of the negligible viscosity of seawater, r _x and r _y are nominally zero. Thus, we seek to construct the field B(x) that provides an optimal balance between the membrane stresses (first two terms on the right-hand side) and the driving stress (last term on the right-hand side) such that r _x and r _y are close to zero.

2.2. Reparameterization for ice rigidity

Before proceeding, we first introduce a commonly used reparameterization of the rigidity B:

(4)$$B( {\bf x}) = B_0( {\bf x}) \, {\rm e}^{\theta ( {\bf x}) },\; $$

where B ₀ (x) represents some reference field of the rigidity and θ(x) corresponds to logarithmic rigidity variations about B ₀ (x). Since B is a strictly positive quantity, this reparameterization allows for the transformation of an inequality-constrained inference problem on B to an unconstrained inference problem on θ while also compressing the dynamic range of rigidity variations to log-space (Shapero and others, Reference Shapero, Badgeley, Hoffman and Joughin2021; Brinkerhoff, Reference Brinkerhoff2022). Analysis of inferred variations in B can thus be performed by inserting any inferred variations in θ into the above equation. Proper choice of the reference rigidity B ₀ depends on the prior information available, which we discuss in detail in later sections.

2.3. Probabilistic inference of ice rigidity from remote-sensing data

Observations of horizontal ice velocity over ice-sheet margins have been widely available for the past decade thanks to the prevalence of remote-sensing platforms and efficient data processing methodologies (Joughin and others, Reference Joughin, Smith, Howat, Scambos and Moon2010; Mouginot and others, Reference Mouginot, Rignot, Scheuchl and Millan2017; Gardner and others, Reference Gardner, Fahnestock and Scambos2019). At the same time, improved integration of ice-penetrating radar and surface velocities using mass conservation techniques have allowed for more accurate and higher resolution maps of ice thickness and bathymetry (Morlighem and others, Reference Morlighem2017). Specifically, over ice shelves, it is common practice to convert observations of surface elevation, which are well constrained, to ice thickness by assuming hydrostatic equilibrium and applying corrections for firn layers derived from in situ thickness data (Morlighem and others, Reference Morlighem2020). Thus, we can estimate spatial gradients and compute the SSA momentum balance directly for a given rigidity field when velocity and thickness observations are spatially continuous over an ice shelf. However, observation noise and data gaps generally degrade estimates of observation gradients, resulting in non-physical estimates of SSA terms that require an additional gradient operation. Therefore, we seek to formulate a method for approximating the continuous functions underlying the surface velocity, ice thickness and ice rigidity fields that balance the reconstruction accuracy of the observed data while resulting in minimal SSA residuals, r. Furthermore, such a method should allow for rigorous quantification of the uncertainties associated with θ, which will be driven by observation noise, varying sensitivities of the SSA equations to spatial variables, and prior knowledge on the range of values for θ.

To that end, we utilize Bayes’ Theorem to construct the joint posterior probability distribution for rigidity and the reconstructed surface velocity and ice thickness, conditioned on a set of velocity and thickness observations assimilated into the SSA momentum balance. Let us first define the observation vector d = [u,  v,  h], reconstructed observation vector $\hat {{\bf d}} = [ \hat {u},\; \, \hat {v},\; \, \hat {h}]$, and SSA residual vector r = [r _x,  r _y]. In the following, we consider that all three vectors vary spatially over an ice shelf, i.e. d = d(x), $\hat {{\bf d}} = \hat {{\bf d}}( {\bf x})$ and r = r(x). We can then write the joint posterior distribution for θ and $\hat {{\bf d}}$ as (Appendix A):

(5)$$p( \theta,\; \hat{{\bf d}} \vert {\bf d},\; {\bf r}) \ \propto\ p( {\bf d} \vert \hat{{\bf d}}) p( {\bf r} \vert \theta,\; \hat{{\bf d}}) p( \theta) .$$

The first two terms on the right-hand side of the equation are the data likelihood distributions, where $p( {\bf d} \vert \hat {{\bf d}})$ measures the probability of observing d for a given set of predictions $\hat {{\bf d}}$ while $p( {\bf r} \vert \theta ,\; \, \hat {{\bf d}})$ measures the probability of computing r (which is nominally [0,  0] for ice shelves) through evaluation of the SSA momentum balance using $\hat {{\bf d}}$ and θ. The last term on the right is the prior distribution p(θ), which encodes our prior knowledge on θ values without having seen any observations.

2.4. Variational inference of ice rigidity

Due to non-linearities in the SSA momentum balance and potentially non-Gaussian data likelihoods, the posterior for θ and $\hat {{\bf d}}$ must be approximated, either by drawing random samples from $p( \theta ,\; \, \hat {{\bf d}} \vert {\bf d},\; \, {\bf r})$ or by constructing a suitable approximating distribution. The former strategy is based on the general class of Markov Chain Monte Carlo (MCMC) approaches and tends to be suitable for a low or moderate number of model dimensions. However, for spatial domains with sizes typical of Antarctic ice shelves, the number of model dimensions will be quite high, regardless of the model used to represent the spatial fields. Therefore, we instead employ a variational inference framework wherein we aim to construct an approximating distribution, $q( \theta ,\; \, \hat {{\bf d}})$, for θ and $\hat {{\bf d}}$ that is minimally divergent from the true posterior $p( \theta ,\; \, \hat {{\bf d}} \vert {\bf d},\; \, {\bf r})$. The choice of approximating distribution generally involves a trade-off between computational complexity (e.g. number of trainable parameters) and accuracy in capturing relevant statistics and correlations in the target posterior (Blei and others, Reference Blei, Kucukelbir and McAuliffe2017). We discuss the choice of approximating distribution in the next section.

At this stage, let us also enforce that $q( \theta ,\; \, \hat {{\bf d}})$ is the product of two independent variational distributions, $q( \theta ,\; \, \hat {{\bf d}}) = q( \theta ) q( \hat {{\bf d}})$. This independence will allow us to use separate machine learning models to reconstruct observations and predict rigidity while using a joint objective function to tune their parameters simultaneously. Moreover, we restrict $q( \hat {{\bf d}})$ to produce only the mean reconstructed velocity and thickness, assuming that information about formal observation uncertainties can be obtained. Later, we will discuss how this assumption on $q( \hat {{\bf d}})$ will affect the posterior variance estimates for θ.

As is commonly done in variational inference, we utilize the Kullback–Leibler (KL) divergence for measuring the difference between two probability distributions:

(6)$${\cal KL}\left[q( \theta,\; \hat{{\bf d}}) \Vert p( \theta,\; \hat{{\bf d}} \vert {\bf d},\; {\bf r}) \right] = \int q( \theta,\; \hat{{\bf d}}) \log {q( \theta,\; \hat{{\bf d}}) \over p( \theta,\; \hat{{\bf d}} \vert {\bf d},\; {\bf r}) }\, d\theta\, d\hat{{\bf d}}.$$

By minimizing the KL-divergence, we are tuning the variational distribution to be close to the target posterior distribution from an informational perspective. However, the existence of the posterior distribution in the denominator of the log term generally makes computing the KL-divergence intractable. As shown in Appendix A, we can instead maximize a stochastic variational lower bound, referred to as the Evidence Lower Bound (ELBO):

(7)$${\text ELBO} = E_{\theta \sim q( \theta) } \left[\log p( {\bf r} \vert \theta,\; \hat{{\bf d}}) \right] + \log p( {\bf d} \vert \hat{{\bf d}}) - {\cal KL}\left[q( \theta) \Vert p( \theta) \right].$$

The first term for the ELBO is the expected value of the log-likelihood for the SSA residual vector, r, for a given set of predictions $\hat {{\bf d}}$ and integrated over the variational distribution for θ (see Appendix A for estimation of this integral over ice shelves). The second t erm for the ELBO is the data log-likelihood for the observations d given predictions $\hat {{\bf d}}$. The last term is the KL-divergence between the variational distribution q(θ) and prior p(θ) for the rigidity, which can be computed analytically for certain pairs of distributions, e.g. two multivariate normal distributions. Overall, the ELBO provides an objective function for learning an approximating distribution for the spatial field of rigidity, as well as the mean ice velocity and thickness fields. In the next section, we discuss the choice of approximating distribution q(θ) and the parameterization of spatial fields such that the ELBO can be computed efficiently for large datasets and spatial domains.

2.5. Parameterization of spatial fields and physics-informed machine learning

We seek representations of the mean spatial fields $\hat {{\bf d}}( {\bf x})$ and an approximating distribution q(θ) that permit efficient evaluation of the ELBO in Eqn (7). For ice-sheet modeling in general, spatial fields are typically discretized into an irregular mesh of triangular finite elements. However, this approach can lead to a large number of nodal parameters for a given spatial field, particularly for fields with high spatial variability with length scales on the order of a few ice thicknesses. Large numbers of nodes can severely restrict the choice of approximating distributions that are computationally feasible. For example, for a mesh with N nodes, an approximating multivariate normal distribution with a mean vector of N elements and a covariance matrix of N ² elements would require excessive computer memory when N is greater than a few thousand nodes.

Our strategy in this work is to adopt mesh-free parameterizations of spatial fields using two different machine learning models: (1) a neural network for the mean fields $\hat {{\bf d}}( {\bf x})$; and (2) a variational Gaussian Process (VGP) for the approximating distribution for the rigidity field θ(x). Both models fall under the broad classification of function approximators for a given set of data, while the VGP (and Gaussian Processes in general) also approximates the probability distribution of functions that generate the given dataset. For the mean fields $\hat {{\bf d}}( {\bf x})$, a dense, feedforward neural network, f _ψ, is used to represent the surface observations on a point-by-point basis:

(8)$$\hat{{\bf d}}_i = f_{{\boldsymbol \psi}}\left({\bf x}_i\right),\; $$

where $\hat {{\bf d}}_i$ is the vector of neural network predictions at the i-th coordinate x_i, and ψ represents the total set of weights and biases of the hidden layers. The choice of feedforward networks in this work is motivated by the ability to generate predictions of $\hat {{\bf d}}_i$ at arbitrary coordinates, which facilitates mini-batch training by stochastic gradient descent. For the ice shelf models investigated here, we found that four-layer feedforward networks with 50–100 nodes per layer provide a suitable balance between expressivity and computational efficiency. Additionally, we use tanh activation functions as they are continuously differentiable and result in spatial gradients that are spatially smooth (Riel and others, Reference Riel, Minchew and Bischoff2021).

For the rigidity variational distribution, we exploit the ability of Gaussian Processes (GPs) to construct an approximating multivariate normal distribution through computation of a mean vector for any set of coordinates and a covariance matrix for any pairwise set of coordinates (Rasmussen, Reference Rasmussen2003). Multivariate normal distributions permit modeling of correlations between variables and have been shown to agree well with MCMC-sampled posteriors for ice dynamics inverse problems (Petra and others, Reference Petra, Martin, Stadler and Ghattas2014). The mean and covariance both require the evaluation of a kernel function that describes the expected similarity between data points for a given coordinate pair. Here, we use a squared exponential kernel:

(9)$$k( {\bf x},\; {\bf x}') = \sigma^2 \exp \left(-{( {\bf x} - {\bf x}') ^T ( {\bf x} - {\bf x}') \over 2L^2}\right),\; $$

where σ ² is the kernel variance and L is a covariance length scale that describes the characteristic distance up to which function values at x and x′ are correlated. However, as described in Appendix B, the use of standard GPs to infer an approximating distribution for a given dataset will be computationally prohibitive for large datasets, which stems from the need to evaluate the kernel at N coordinates in order to construct and invert a covariance matrix of size N × N. For our study areas, N is usually on the order of tens of thousands of points. Variational GPs are a modification of standard GPs that utilize the concept of sparse inducing index points for overcoming the prohibitive ${\cal O}( N^2)$ memory and ${\cal O}( N^3)$ time requirements of large sets of index points for standard GPs (Titsias, Reference Titsias2009). Rather than evaluating the kernel at N coordinates, we evaluate the kernel at M inducing index coordinates, z, where M ≪ N. The inducing index points are sometimes referred to as support points and act as a low-dimensional latent representation of the data (Quinonero-Candela and Rasmussen, Reference Quinonero-Candela and Rasmussen2005). As explained in Appendix B, computation of a posterior mean and covariance matrix at an arbitrary number of prediction points involves evaluation of the kernel at the prediction and inducing index points and combining the kernel evaluations with a full-rank covariance matrix at the inducing points. The number of inducing points, M, is a prescribed parameter that provides additional control on the complexity of the predicted variable field. For the ice shelves studied here, we find that M between 300 and 1000 provides sufficient expressiveness for the rigidity fields. Overall, the total set of trainable variables, ϕ, for the VGP consists of the inducing point locations z, the inducing point mean values, the covariance matrix values at the inducing points and the hyperparameters of the kernel (σ ² and L). In the following, we denote the approximating multivariate normal distribution as q _ϕ(θ).

While GPs and VGPs are usually trained in a supervised fashion for a training dataset of N examples, we do not have a set of ground truth rigidity values for θ. We instead insert q _ϕ(θ) into the ELBO in Eqn ( 7) in order to train f _ψ and q _ϕ(θ) simultaneously. This style of training is consistent with recent developments in physics-informed machine learning where available physics knowledge (in our case, the SSA momentum balance) is used to augment observables with pseudo-observables, which can be useful in domains where direct observations are impossible or difficult to obtain (Raissi and others, Reference Raissi, Perdikaris and Karniadakis2019; Riel and others, Reference Riel, Minchew and Bischoff2021). Here, ice surface velocity and thickness are observable, and their reconstructed values (generated by f _ψ) are used in conjunction with rigidity samples drawn from q _ϕ(θ) in order to compute the SSA residual pseudo-observables, which should nominally be zero. We evaluate the SSA residuals at a random set of C ‘collocation’ coordinates, x_c, not included in the training data for f _ψ in order to maximize the information extracted from the pseudo-observables (Raissi and others, Reference Raissi, Perdikaris and Karniadakis2019).

The ELBO in Eqn (7) can now be written in terms of the neural network parameters ψ and the VGP parameters ϕ:

(10)$$\eqalign{{\text ELBO}( {\boldsymbol \psi},\; {\boldsymbol \varphi}) = & E_{\theta_c \sim q_{{\boldsymbol \varphi}}( \theta_c) } \left[\log p( {\bf r} \vert \theta_c,\; f_{{\boldsymbol \psi}}( {\bf x}_c) ) \right] + \log p( {\bf d} \vert f_{{\boldsymbol \psi}}( {\bf x}) ) \cr & - {\cal KL}\left[q_{{\boldsymbol \varphi}}( \theta_z) \Vert p( \theta_z) \right],\; }$$

where we use the shorthand notation θ(x_c) = θ _c and θ(z) = θ _z to indicate the rigidity field evaluated at the collocation points x_c and the inducing index points z, respectively. The neural network f _ψ is evaluated at the surface observation coordinates x. In this work, we minimize the negative of the ELBO using the Adam optimizer (Kingma and Ba, Reference Kingma and Ba2014) with a learning rate between 1e⁻⁴ and 5e⁻⁴ and batch sizes varying from 256 to 1024 examples. Additionally, we often found it beneficial to pre-train the network f _ψ using only the middle term of the ELBO (i.e. fitting the observations without the rigidity field). The pre-training results in a f _ψ that can then be fine-tuned with the full ELBO in order to improve overall training convergence with the VGP. Typically, we make a qualitative assessment of convergence by monitoring the loss values for training and testing datasets (using a train-test split of 85 and 15% of the data, respectively) and stop training when reductions in the loss value are negligible (Fig. S1). A summary of the neural network and variational inference training framework is described in Figure 1.

Fig. 1. Illustration of physics-informed variational inference framework. The two main components of the framework are the neural network f _ψ tasked with reconstructing surface observations (grey box) and the variational Gaussian Process (VGP) tasked with generating the distribution q(θ) (white box), which is a variational approximation to the posterior distribution of the normalized ice rigidity parameter θ. The training loss function is the sum of: (1) a data likelihood loss measuring the consistency between the reconstructed (predicted) and observed observations; (2) the expected value of the SSA residual likelihood given predicted rigidity values and observations; and (3) the KL-divergence between the variational and prior distributions for rigidity. For the data likelihood, training data and corresponding uncertainties and spatial coordinates are sampled from remote-sensing observations (red dots). For the SSA residual likelihood, an independent set of collocation coordinates (x _c and y _c, blue dots) are sampled from the model domain, which are input to f _ψ and the VGP in order to evaluate the SSA momentum balance at those coordinates.

2.6. Selection of data likelihood and prior distributions

To concretize the ELBO training objective for the neural network f _ψ and VGP q _ϕ(θ), we now discuss specification of the data likelihoods for the surface observations and SSA residuals and the prior distribution for the rigidity.

2.6.1. Data likelihood specification

For the likelihood $p( {\bf d} \vert \hat {{\bf d}})$, we use independent normal distributions for the velocity components and ice thickness:

(11)$$\eqalign{\,p( {\bf d} \vert \hat{{\bf d}}) & = p( u \vert \hat{u}) p( v \vert \hat{v}) p( h \vert \hat{h}) ,\; \cr p( u \vert \hat{u}) & = {\cal N}( u;\; \hat{u},\; \sigma_{u}^2) ,\; \cr p( v \vert \hat{v}) & = {\cal N}( v;\; \hat{v},\; \sigma_{v}^2) ,\; \cr p( h \vert \hat{h}) & = {\cal N}( h;\; \hat{h},\; \sigma_{h}^2) ,\; }$$

where the different σ ² variables correspond to the variances of each observable. The observation variances may be prescribed or learned, and in this work, we prescribe the variances to be scaled values of formal observation uncertainties provided with the datasets. We found that a scaling factor between 10 and 30 prevented the neural network f _ψ from overfitting the velocity data in order to generate velocities that are more consistent with the SSA momentum balance (e.g. mitigating high spatial frequency velocity variations that are not well-modeled by the SSA approximation). We further note that we apply a limited amount of spatial smoothing to the velocity and thickness data (using a Gaussian filter with a window size of roughly half of the mean ice thickness of the shelf) in order to mitigate high-frequency observation noise and improve training convergence.

For the likelihood $p( {\bf r} \vert \theta ,\; \, \hat {{\bf d}})$, we again use a normal distribution:

(12)$$p( {\bf r} \vert \theta,\; \hat{{\bf d}}) = {\cal N}( {\bf r}( \theta,\; \hat{{\bf d}}) ;\; {\bf 0},\; \sigma_{{\bf r}}^2) ,\; $$

where the mean of zero encodes that we expect ${\bf r} = {\bf r}( \theta ,\; \, \hat {{\bf d}})$ to be zero on ice shelves, and $\sigma _{{\bf r}}^2$ is the expected variance of the residuals. Prescribing a proper value of $\sigma _{{\bf r}}^2$ for the likelihood distribution involves a careful consideration of the observation uncertainties and the expected uncertainties in θ. In the general case, one could perform a Monte Carlo estimate of $\sigma _{{\bf r}}^2$ by inputting random realizations of the velocity data and samples of θ from the prior into the SSA equations to get an expected range of values for r. Using this approach, we find that σ _r ≈ 1 kPa works well for most cases and can be decreased for lower observation noise. Note that this approach does not explicitly consider situations where the SSA (with zero basal drag for ice shelves) is expected to perform poorly, such as pinning points where ice shelves become locally grounded over bathymetric highs. These epistemic uncertainties should correspond to an increase in $\sigma _{{\bf r}}^2$. An additional improvement to the methods presented here is to use a full variational distribution for $q( \hat {{\bf d}})$, rather than using only the mean values $\hat {{\bf d}}$. Our restriction to the mean values for $\hat {{\bf d}}$ likely reduces the estimated posterior variance for θ, which can be partially compensated by inflating $\sigma _{{\bf r}}^2$. Overall, the likelihood formulation can be improved by modeling the full probability distribution for the surface observations, explicit handling of epistemic uncertainties, incorporation of spatially correlated uncertainties (e.g. Brinkerhoff, Reference Brinkerhoff2022) and using a non-Gaussian probability distribution. We leave these improvements for future exploration.

2.6.2. Prior distribution

As discussed earlier, the VGP uses a squared exponential kernel function for posterior inference of rigidity (Eqn 9). Likewise, we use the squared exponential kernel to construct the GP prior p(θ) in order to encourage spatial coherence between predictions of θ at different coordinates:

(13)$$k_{\text{prior}}( {\bf x},\; {\bf x}') = \sigma_\theta^2 \exp \left(-{( {\bf x} - {\bf x}') ^T ( {\bf x} - {\bf x}') \over 2L_\theta^2}\right),\; $$

where the prior variance and length scale, $\sigma _\theta ^2$ and $L_\theta$, are prescribed and not tuned during minimization of the negative ELBO. For ice dynamics well-described by the SSA approximation, L _θ can usually be chosen to be some finite multiple of the mean ice shelf thickness, which is commensurate with the longitudinal stress coupling length scale (Cuffey and Paterson, Reference Cuffey and Paterson2010). The complete prior distribution is written as:

(14)$$p( \theta( {\bf x}) ) = {\cal N}( {\bf 0},\; k_{\text{prior}}( {\bf x},\; {\bf x}) ) .$$

Since we enforce a prior mean of zero, the prior variance $\sigma _\theta ^2$ will depend on the value of B ₀ used for the reparameterization of rigidity. In this work, we investigate two different strategies: (i) assume that B ₀ is uniform; or (ii) estimating a spatial field B ₀(x) independently through an inversion using traditional control methods. For the second strategy, recall that traditional control methods generally use an optimization objective based on the misfit between observed and predicted ice velocities, which is different from the momentum balance optimization objective used here. Therefore, the velocity misfit-based objective will implicitly have different spatially varying sensitives to the parameter field than the momentum balance objective, which provides an opportunity to combine the two objectives in a complementary manner. For the two different strategies for selecting B ₀, we also correspondingly adjust the prior variance. For a uniform B ₀, we set $\sigma ^2_\theta = 1$ to allow for a relatively large variation in θ over the ice shelf. For a B ₀ obtained from a control method inversion, we reduce $\sigma _\theta ^2$ to 0.2, which encodes our belief that the values of B from the inversion are relatively well-constrained, and θ thus represents smaller deviations of B dictated by the momentum balance optimization objective.

2.7. Generating shelf-wide samples of the ice rigidity

While the VGP utilizes inducing index points to allow for per-batch prediction of θ and the corresponding posterior covariance matrices, we require an algorithm for generating a random sample of θ with a given spatial resolution over the entire modeling domain. Even with the use of inducing points, assembling global posterior mean and covariance matrices for the entire domain would require a very large amount of memory for uniform grids with sizes exceeding tens of thousands of grid points (${\cal O}( NM)$ memory requirements for N grid points and M inducing points). We therefore utilize an MCMC-based approach where a random sample of θ is generated at a limited set of coordinates within the ice shelf and used as a seed to grow a full chain over the entire shelf, which is equivalent to block-sampling approaches for sampling Gaussian Markov Random Fields (e.g. Rue, Reference Rue2001). Note that this form of MCMC utilizes the posterior statistics learned from variational inference and does not involve evaluation of a physics-based forward model. Specifically, we apply Gibbs sampling on a block-by-block basis where a block is defined as a small subset of the uniform grid. For each block, the mean ${\boldsymbol \mu }_\theta$ and covariance matrix ${\bf \Sigma }_\theta$ are computed using the trained VGP. In this way, we can grow a random chain through successive evaluation of smaller multivariate normal distributions rather using a multivariate normal with a very large global covariance matrix. As an example, assume that the first two blocks are combined in the following partition:

(15)$${\bf x} = \left[\matrix{{\bf x}_1 \cr {\bf x}_2 }\right],\; \quad {\boldsymbol \theta} = \left[\matrix{{\boldsymbol \theta}_1 \cr {\boldsymbol \theta}_2 }\right],\; \quad {\boldsymbol \mu}_\theta = \left[\matrix{{\boldsymbol \mu}_1 \cr {\boldsymbol \mu}_2 }\right],\; \quad {\bf \Sigma}_\theta = \left[\matrix{{\bf \Sigma}_{11} & {\bf \Sigma}_{12} \cr {\bf \Sigma}_{21} & {\bf \Sigma}_{22} }\right]$$

Then, initialization of a Gibbs chain would use the following identities:

(16a)$${\boldsymbol \mu}_\theta,\; \enspace {\bf \Sigma}_\theta \leftarrow q_{{\boldsymbol \varphi}}( \theta( {\bf x}) ) ,\; $$

(16b)$${\boldsymbol \theta}_1 \sim {\cal N}( {\boldsymbol \mu}_1,\; {\bf \Sigma}_{11}) ,\; $$

(16c)$${\boldsymbol \theta}_2 \sim {\cal N}\left({\boldsymbol \mu}_2 + {\bf \Sigma}_{21} {\bf \Sigma}_{11}^{-1} ( {\boldsymbol \theta}_1 - {\boldsymbol \mu}_1) ,\; {\bf \Sigma}_{22} - {\bf \Sigma}_{21}{\bf \Sigma}_{11}^{-1} {\bf \Sigma}_{12}\right).$$

In words, we use the trained VGP to generate a mean vector and covariance matrix for two blocks of coordinates, which are then partitioned according to Eqn (15). A sample θ ₁ = θ(x₁) is drawn directly from the mean and covariance of the first block. Then, a sample θ ₂ = θ(x₂) is drawn from a multivariate normal with a mean and covariance that are computed using the Schur complement. We then set μ₁ ← μ₂, θ₁ ← θ₂, Σ₁₁ ← Σ₂₂, and use the VGP to predict the mean and covariance for the next set of blocks. The process is repeated for all blocks in the modeling domain. With this approach, each sample is independent of the previous sample and is accepted with a probability of 1. In our experiments, we found that if the block size was chosen to be too small, the variance of the final sample was artificially large, likely due to excessive truncation of the covariance matrix (depending on the resolution of the uniform grid relative to the length scale of the posterior variations in rigidity). Here, we found a block size of 1000 grid points works well for grids with cell sizes equal to roughly half or a quarter of the prior covariance lengthscale. Overall, we can validate the samples derived from the Gibbs sampler by viewing traceplots (plots of samples at independent coordinates) and cross-comparing the standard deviation fields between the sample standard deviation and the standard deviation predicted directly by the VGP (Fig. S6).

2.8. Related work

Bayesian inference has long been applied to geophysical inverse problems, and as computational resources and inference algorithms improve, the complexity and size of the physical models investigated has increased. Within glaciological inverse problems, Bayesian formulations of the posterior distributions have been used as cost functions for obtaining point estimates of basal topography and friction for grounded ice streams (Pralong and Gudmundsson, Reference Pralong and Gudmundsson2011). For fully Bayesian inference, Petra and others (Reference Petra, Martin, Stadler and Ghattas2014) developed an MCMC method for estimating the posterior distribution for ice-sheet models with a large number of parameters, utilizing low-rank approximations of data likelihood Hessian matrices in order to reduce computational complexity while improving sample efficiency. Similarly, Gopalan and others (Reference Gopalan, Hrafnkelsson, Aḥhalgeirsdóttir and Pálsson2021) used a Gibbs sampler in order to sample for ice stream model parameters for a simpler model applicable to slower-flowing ice. While MCMC methods generally serve as ‘gold standards’ for Bayesian inference, they do not scale well to large problem sizes. MCMC methods that invoke simpler proposal distributions usually require many more samples in order to sufficiently sample the posterior, whereas methods that can utilize the problem structure to improve sample efficiency require more computational resources (Petra and others, Reference Petra, Martin, Stadler and Ghattas2014).

Methods that approximate the posterior distribution, rather than sample from it, provide appealing alternatives to MCMC. Both Isaac and others (Reference Isaac, Petra, Stadler and Ghattas2015) and Babaniyi and others (Reference Babaniyi, Nicholson, Villa and Petra2021) utilize a Gaussian approximation of the posterior centered on the maximum a posteriori (MAP) point (i.e. a Laplace approximation) in order to infer basal drag parameters for ice sheets. While Laplace approximations subvert the need for generating posterior samples (and the forward model evaluations associated with each sample), they can lead to posterior approximations that fail to capture much of the probability mass when the posterior is sufficiently non-Gaussian or multi-modal (Penny and others, Reference Penny, Kiebel and Friston2007). In contrast, variational methods that utilize the KL-divergence as an optimization criterion (as done here) tend to favor approximating distributions that match the moments of the target distribution (e.g. mean and variance), which tends to capture more probability mass. Sufficient capturing of probability mass can be especially important for posterior predictive modeling where non-linearities can lead to a large spread of predictions (e.g. see the section on ice shelf buttressing).

To that end, Brinkerhoff (Reference Brinkerhoff2022) introduced a variational inference method to jointly infer basal drag and ice rheology at a catchment-scale for glaciers. Importantly, the KL-divergence was used to estimate an optimal approximating distribution that also uses a Gaussian process prior, similar to the approximating distributions used in our work. A finite number of eigenvectors of the prior covariance are used to construct a linear model that permits inference at a lower dimension, similar to the sparse inducing points used in the VGP. The construction of the eigenvectors utilizes a coarse grid in the Fourier domain to model signals with a discrete set of spatial frequencies. Thus, the method of Brinkerhoff (Reference Brinkerhoff2022) shares many of the same features proposed here, with two main differences. Firstly, we take the ice surface velocity and thickness as given (subject to smoothing) and use the momentum balance based on the SSA as our forward model in order to compute r. On the other hand, all of the previous approaches strictly enforce r = 0 (i.e. they take satisfaction of the SSA momentum balance as a constraint) and use predicted velocities as the forward model. The difference between the two approaches may be interpreted as taking two different pathways to the same solution, where the difference in forward models will lead to different sensitivities to the ice rheology parameters (see Discussion section). Furthermore, our parameterization of the spatial fields with neural networks and VGPs permits independent evaluation of the spatial gradients in the SSA momentum balance using automatic differentiation. This independence allows for batch-based stochastic gradient descent, which is advantageous for large datasets.

3.1. 1D ice shelf

We first simulate a laterally confined ice shelf using 1D SSA equations where lateral drag is parameterized assuming a rectangular bed with width w in the across-flow direction (Nick and others, Reference Nick, Van der Veen, Vieli and Benn2010). The momentum equation in the along-flow direction reduces to:

(17)$${\partial\over \partial x}\left(4 \eta h {\partial u\over \partial x}\right)- {2Bh\over w}\left({5u\over w}\right)^{1/n} = \rho_igh{\partial s\over \partial x}.$$

As there are no time derivatives in the momentum balance for ice flow, the above equation is solved for u for a given thickness profile h. The thickness is then evolved using mass conservation:

(18)$${\partial h\over \partial t} = a - {\partial q\over \partial x},\; $$

where a is the accumulation rate (set to 0 for this example) and q = uh is the horizontal ice flux through a vertical column of ice.

3.1.2. Evaluation of variational inference

As is commonly done in studies investigating variational inference techniques to approximate a target posterior distribution, we compare the estimated variational distribution with direct samples from the posterior using MCMC. Here, we utilize a No U-Turn Sampler scheme implemented in the NumPyro Python package (Phan and others, Reference Phan, Pradhan and Jankowiak2019), which uses automatic differentiation to efficiently generate sample trajectories for moderately high numbers of model parameters. We use the velocity and thickness predictions from the neural network as the observations such that MCMC only samples the rigidity posterior distribution, which allows for a more direct comparison between the variational and sampled posterior.

We find that both MCMC and variational inference recover a posterior mean profile for θ that is close to the true values for areas >20 km upstream from the ice front (Fig. 2). Close to the ice front, both methods predict uncertainties that are substantially larger due to thinner ice, which reduces the sensitivity of the longitudinal and lateral membrane stresses (left-hand terms in Eqn (17)) to rigidity variations. These uncertainties near the ice front can likely be reduced by augmenting the SSA residual loss with a dynamic boundary condition at the ice front that balances longitudinal stresses normal to the ice front with the difference in hydrostatic pressure between the ice and ocean water (Larour and others, Reference Larour, Rignot, Joughin and Aubry2005). Pair plots of marginal distributions of θ at different locations along the ice shelf show that the variational approach is able to recover strong covariances between θ samples for locations that are relatively close to each other while ensuring samples are uncorrelated for larger pair-wise distances. In general, the strength of the posterior covariance will be modulated by the physical model as well as the prior correlation lengthscale. Closer to the ice front, the marginal distributions derived from MCMC indicate a slight deviation from Gaussian behavior, which is again likely due to the lower ice thicknesses limiting ice stress sensitivity to rigidity variations. Since the variational distribution is constrained to be a multivariate normal, it is unable to recover the skewed, non-Gaussian behavior in the marginals in these areas. For applications where it is desirable to place more probability mass in the longer-tails of the posterior distribution, one could simply increase relevant variances ($\sigma _\theta$, σ _r) to inflate the variational posterior variances or use a more flexible variational approximation not restricted to multivariate normals (Rezende and Mohamed, Reference Rezende and Mohamed2015).

Fig. 2. Posterior inference of ice rigidity for simulated 1D ice shelf. (a) Posterior marginal distributions for θ at different locations along the ice shelf. The grounding line is at X = 0 km and the ice front is at X = 100 km. Diagonal plots show the 1D marginals computed from posterior samples generated with MCMC (blue) and the variational Gaussian process (VGP; orange). Off-diagonal plots show 2D covariance plots for the same sample set. All marginals have been smoothed using a Gaussian kernel density estimator. (b) Velocity (blue) and ice thickness (orange) of ice shelf used for posterior inference. (c) Comparison of the true θ(X) against the mean θ(X) computed from the MCMC samples (blue) and the VGP (orange). The shaded regions correspond to 2σ posterior uncertainties. Overall, the posterior distributions for MCMC and VGP are very similar. The largest deviations occur near the ice front where the marginals exhibit stronger non-Gaussian behavior, which cannot be modeled by the VGP.

3.2. 2D ice shelf

We now simulate a 2D ice shelf using the icepack ice flow modeling software (Shapero and others, Reference Shapero, Badgeley, Hoffman and Joughin2021). Similar to the synthetic ice shelf presented in Shapero and others (Reference Shapero, Badgeley, Hoffman and Joughin2021), we prescribe a semi-circular shelf geometry with four inlet glaciers of varying widths (Fig. 3). Additionally, we prescribe a bed topography that results in a few pinning points where the flotation height is positive, i.e. the ice is actually grounded at these locations. Under the shallow-stream approximation, we prescribe a basal drag friction coefficient proportional to the flotation height such that friction is only non-zero for grounded ice. Such pinning points in the form of ice rumples are common in ice shelves in Antarctica. However, assuming a fully floating ice shelf during inversion for rheological parameters will introduce errors into the inferred parameter field due to model mismatch. Therefore, by purposefully injecting modeling errors into the estimation procedure, we can assess how the two different cost functions and the estimated parameter uncertainties respond to such errors.

Fig. 3. 2D ice shelf with simulated damage evolution and pinning points. The simulation outputs shown are computed from 700 years of spinup in order to achieve steady state. The ice shelf is fed by four inlet ice streams, as evidenced by the flow speed (a) and ice thickness (b). Height-above-flotation (HAF) in (c) shows the location and orientation of the prescribed pinning points. The steady-state ice rigidity B (d) reflects damage accumulation due to shear margin weakening and ice thinning due to large strain-rates over the pinning points. The effective strain rate (e) and effective dynamic viscosity (f) are approximately inversely related and show strong shearing in the ice between the inlet flow, as well as over the pinning points. Strain rates are lower closer to the ice front. The effective viscosity exhibits a mix of long-wavelength variations within flow units and short-wavelength variations near the shear margins.

We use icepack to first simulate the evolution of shelf velocity and thickness for roughly 500 years with a constant ice rigidity, B ₀, corresponding to an ice temperature of $-5\, ^\circ$C, a net surface mass balance of 0, and a stress exponent of n = 4. Here, we choose n = 4 in order to evaluate the sensitivity of the rigidity inference to 2D ice stress variations that are more likely to be found in natural environments of fast-flowing ice (Bons and others, Reference Bons2018; Millstein and others, Reference Millstein, Minchew and Pegler2022). After the first simulation stage, we apply a continuum damage mechanics model that modulates the rigidity field with an evolving damage factor, D, such that B _D = (1 − D)B ₀ (Borstad and others, Reference Borstad, Rignot, Mouginot and Schodlok2013; Albrecht and Levermann, Reference Albrecht and Levermann2014). This approach provides a physically realistic means to obtain a spatially varying prefactor field with rheology-modifying processes such as shear weakening. We run the damage-enhanced model for an additional 200 years to achieve approximate steady-state. At the end of the simulation, we can observe substantial spatial variation in damage, where ice is nearly undamaged at the grounding line (due to a zero-damage boundary condition) and highly damaged near the ice front, at shear margins and downstream of the pinning points. The dynamic effective viscosity field shows concentrated low viscosities near the pinning points and higher viscosities between the inlet ice streams where deformation rates are lower. Overall, the viscosities exhibit a mix of short- and long-wavelength features, which are mirrored in the effective strain rate field.

3.2.2. Evaluation of rigidity reconstruction errors

A more detailed comparison of the recovery error for B (recovered using $B = B_0\, {\rm e}^\theta$) between the control method inversion and variational inference reveals that the two methods are complementary. The control method inversion has the lowest overall error bias, but the areas where the errors are largest are systematically upstream of the pinning points (Fig. 4). Since we assume all ice is floating for the forward model, the missing resistive stress provided by drag at grounded ice is compensated by artificially making the ice stiffer upstream of the pinning points, which acts to slow the ice down in a manner that allows the predicted velocities to match the observed velocities. In contrast, the B recovered by variational inference (which uses the SSA momentum balance as a forward model) shows larger errors directly over the pinning points, as well as in areas where ice is stagnant (low strain rate). Over the pinning points, the true rheology sharply varies from about 250 to 200 prior yr^1/4 kPa (Fig. 3d). However, the prior lengthscale of 15 km encourages spatially smoother fields of ice rigidity, which limits the dynamic range of ice stresses that can be modeled in order to satisfy the SSA momentum balance. Since the driving stress variations over the pinning points are overestimated due to the buoyancy assumption (Fig. S2), the preferred solution is to smooth out all stress variations over the grounded ice in order to minimize the residual SSA terms. Upstream of the pinning points and closer to the grounding line, the recovery errors are actually lower using variational inference as compared to the control method inversion. The spatial patterns in the recovery errors are similar to the patterns of residual SSA components (Fig. S4). Finally, by using the control method inversion as the reference B ₀ for variational inference, we can minimize much of the recovery errors closer to the ice front and in areas where strain rates are lower but flow speeds are still high, i.e. areas where the inversion has greater sensitivity and where the dynamic boundary condition at the ice-ocean interface provides additional constraints on the rigidity (Fig. 4c).

Fig. 4. Comparison of reconstruction errors for mean inferred ice rigidity between control method inversion and the proposed variational inference method. (a) Error (B − B _true) for control method inversion using icepack. (b) Error for variational inference with a uniform B ₀ field. (c) Error for variational inference using the control method inversion for B ₀. (d) Inferred uncertainty for normalized rigidity parameter θ using the control method inversion for B ₀. Black dashed lines correspond to a thickness contour of 150 m while the white dashed lines correspond to an effective strain rate contour of 10^−2.6 a⁻¹. (e) Histograms of errors for different methods. Higher reconstruction errors and uncertainties are mostly concentrated in thinner ice and areas with lower effective strain rates.

3.2.3. Evaluation of rigidity uncertainties

The predicted uncertainties for θ are consistent with the reconstruction errors: uncertainties are higher closer to the ice front where ice thicknesses are lower (as observed in the 1D case), as well as in more stagnant ice where strain rates are lower (Figs 4d, 5). Uncertainties near the ice front are reduced when using B ₀ from the control method inversion (Fig. S3), which is again likely due to the incorporation of the dynamic boundary condition at the ice-ocean interface which also reduced reconstruction errors. In areas where ice is thinner but strain rates are higher (e.g. higher shear strain rate in the areas between the fast-flowing ice), the balance between extensional stresses and lateral drag also provides sufficient signal for reducing uncertainties. In a few isolated patches, even when effective strain rates are low and ice is relatively thin, slightly positive lateral forces that act as a ‘pull’ on the ice can also reduce uncertainties (Fig. 5). Directly over the pinning points, uncertainties are low due to the high strain rates there. This mismatch between low uncertainties and high SSA residuals over the pinning points may be compensated by inflating the variance $\sigma _{\bf r}^2$ over areas where measured flotation height is near zero (see Discussion section). Downstream of the pinning points, we observe higher uncertainties due to downstream thinning of the ice. Overall, the uncertainty maps for θ indicate that areas with thicker ice and higher strain rates are better constrained, and targeted inverse modeling (e.g. estimating B ₀ from a control method inversion) can be an effective tool for further reducing uncertainties (Fig. S3).

Fig. 5. Uncertainty for normalized rigidity θ vs ice thickness, along-flow lateral drag and effective strain rate for simulated 2D ice shelf. (a) Effective strain rate vs ice thickness with colors corresponding to uncertainty in θ. (b) Same as (a) but for effective strain rate vs along-flow lateral drag. Here, lateral drag is computed as the transverse gradients of the SSA momentum balance projected to the along-flow direction, where negative values denote flow resistance. While ice thickness is the first-order control on rigidity uncertainties, higher strain rates can reduce uncertainties in thinner ice. Positive lateral forces can also reduce uncertainties where effective strain rates are low.

4.1. Remote-sensing data and pre-processing strategy

For RIS and B-SW-RL, we use the MEaSUREs velocity mosaic (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011; Mouginot and others, Reference Mouginot, Scheuchl and Rignot2012), which combines speckle tracking of SAR images from various satellite platforms with feature tracking of Landsat 8 images and has a nominal temporal coverage between 2009 and 2016. For LCIS and FRIS, we use a 2020 annual velocity mosaic provided by ITS_LIVE, which is derived from feature tracking of Landsat 7 and 8 images over Antarctica (Gardner and others, Reference Gardner, Fahnestock and Scambos2019). From a visual inspection, we found that the ITS_LIVE mosaic exhibited fewer velocity artifacts for LCIS and FRIS, whereas the MEaSUREs mosaic exhibited fewer artifacts over B-SW-RL and provided full coverage over RIS. Ice thickness and elevation data are derived from BedMachine V2 (Morlighem and others, Reference Morlighem2020), which combines radar-estimated thickness profiles with mass conservation constraints and firn corrections in order to obtain continuous thickness maps. While the nominal year for the thickness data is 2015, the correspondence between the velocity and thickness data are sufficient for the spatial resolution of our analysis (assuming an upper bound of ≈5 km of motion for feature advection). For all velocity and thickness rasters, we first perform a void-filling operation that uses a spring-based PDE constraint to fill in missing data (D'Errico, Reference D'Errico2012). The rasters are then filtered to ≈10−15 times the average ice thickness using a Savitzky–Golay filter in order to remove high-frequency components not resolvable by the SSA momentum balance.

4.2. Variational inference setup

As with the simulated ice shelf, we first invert for B using icepack in order to optimize a cost function combining a velocity misfit term (weighted by the formal uncertainties for the velocity estimates) and a regularization term based on first-order spatial gradients to encourage smoother solutions. All ice shelves are discretized into 2D triangular finite elements with an average size of 5 km. A penalty parameter controlling the relative contribution of the regularization term is selected with a standard L-curve analysis, independently for each ice shelf. For each ice shelf, we use feedforward neural networks with four layers of 100 nodes each and VGPs with 600–900 inducing index points. The training data for the feedforward neural networks and the collocation points for the VGP are independently constructed with a point density of 2 points per 10 square kilometers, which is commensurate with the size of the mesh elements used for the icepack inversion. This point density for the VGP is sufficient to model a moderate-resolution (~20 km) θ field, which when combined with the higher-resolution B ₀ from the control method inversion results in rigidity samples that resolve important variations near shear margins and rifts. The neural network training observations are randomly sampled directly from the raster grids, whereas the collocation points are randomly generated using a latin hypercube sampling scheme. Finally, we use the estimated B field as the reference field B ₀ for parameterization of θ, setting the prior variance for θ to 0.2².

4.3. Analysis of rigidity mean and uncertainty fields

To a first-order approximation, ice is inferred to be stiffer for FRIS and RIS than for LCIS and BWSRL, and average rigidity values for LCIS are the lowest of the four (Fig. 6). These first-order trends are well-matched by modeled ice shelf surface temperatures where temperatures for FRIS are generally around −25 to $-30\, ^\circ$C, whereas for LCIS they range from −15 to $-10\, ^\circ$C (Fig. S8). However, all ice shelves exhibit significant spatial variability in inferred ice rigidity beyond surface temperature variations in order to minimize SSA residuals (Figs S9, S10). For FRIS, the estimated mean B field is broadly consistent with results from prior studies (e.g. MacAyeal and others, Reference MacAyeal, Rignot and Hulbe1998; Larour and others, Reference Larour, Rignot, Joughin and Aubry2005). Ice is inferred to be substantially softer in the shear margins where strong lateral shearing leads to viscous dissipation and elevated ice temperatures. These shear margins are prominent in the Ronne Ice Shelf where fast-flowing floating ice is in contact with rock (along the Orville coast and Berkner Island) or stagnant ice, as is the case downstream of the Korff Ice Rise. As discussed in Larour and others (Reference Larour, Rignot, Joughin and Aubry2005), larger basal melt rates on the northern tip of the Henry Ice Rise are coincident with softer ice. Within the Filchner Ice Shelf, lower overall values of B indicate softer ice, again in the shear margins where ice streams flow onto the shelf and are in contact with stagnant ice. A large lateral surface crevasse close to the ice front is also associated with higher strain rates and softer ice. We can also observe localized regions of substantially stiffer ice, such as downstream of the Foundation Ice Stream and upstream of the Korff and Henry Ice Rises. These regions are associated with larger driving stresses (Fig. S7) such that ice is inferred to be stiffer in order to provide enough resistive stresses to balance those driving stresses. Ice is also inferred to be stiffer closer to the grounding line where colder ice is advected by the ice streams. For all ice shelves, rigidity uncertainties are mostly lower where ice is thicker and strain rates are larger, similar to what was observed for the simulated 1D and 2D ice shelves (Fig. 7).

Fig. 7. Estimated 1-σ B uncertainties for West Antarctic ice shelves. Uncertainties are generally larger for higher B values (scale-dependence) and for areas with thinner ice and lower driving stresses. Uncertainties tend to be lower closer to the grounding line.

Similar to FRIS, the ice in the central portions of RIS are inferred to be more rigid, likely due to relatively cold surface temperatures of −20 $^\circ$C. However, we can also observe zones of softer ice near shear margins and localized areas of grounding. At the inlet of the Byrd Glacier to the west, prominent shear margins separating the fast-flowing inlet ice from more stagnant shelf ice are coherent for more than 300 km downstream of the grounding line (Fig. S7), which results in substantial shear weakening. In the central trunk of the Byrd Glacier inlet, the reduction in flow speed as the ice flows onto RIS leads to enhanced compressional stress and thickening of the ice, leading to inferred higher B values. On the east side of RIS, the Shirase Coast Ice Rumples (SCIR) at the outlet of the MacAyeal and Bindschadler Ice Streams significantly modify the flow field and ice thickness due to grounding of the ice, consistent with the simulated pinning points for the synthetic ice shelf. Thinning of ice downstream of SCIR and diversion of the shear margins toward Roosevelt Island (RI) are both dynamical effects that modify the buttressing capability of ice in this region (Still and others, Reference Still, Campbell and Hulbe2019; Still and Hulbe, Reference Still and Hulbe2021). In our inferred B field, the ice covering the rumples is inferred to be softer while the downstream ice connected to RI is inferred to be stiffer. Alternatively, the ice upstream of Steershead Ice Rise (SIR) is near-stagnant, leading to very high inferred values for B. Downstream of SIR is a streakline of thin ice coincident with the shear margin of the inlet of MacAyeal and Bindschadler Ice Streams, leading to a narrow zone of soft ice that persists nearly all the way to the ice front.

At LCIS, the softest ice is inferred within highly localized areas corresponding to surface crevassing, including the large rift originating from the Gipps Ice Rise (Khazendar and others, Reference Khazendar, Rignot and Larour2011; Larour and others, Reference Larour, Rignot, Poinelli and Scheuchl2021). It is likely that some fraction of the inferred softness is due to not explicitly including rifts (geometrically and dynamically) within the ice flow model, which can reproduce a significant proportion of the observed strain rates with active opening/closing of rifts (Larour and others, Reference Larour, Rignot, Poinelli and Scheuchl2021). As is the case with FRIS, stiffer ice is inferred near the grounding line where colder and thicker ice is advected downstream by the inlet ice streams. Within the ice shelf, areas in between faster flowing ice correspond to thinner ice and higher strain rates, resulting in softer ice. Unlike FRIS, the proximity of the fast flowing inlet ice streams with one another limits the areal extent of stagnant ice over Larsen C. High effective strain rates between ice streams are aligned with the initiation of suture zones where mechanically weak marine ice (sourced from warmer ocean water) has been observed to accumulate at the base of LCIS (Kulessa and others, Reference Kulessa, Jansen, Luckman, King and Sammonds2014). The initial portion of the suture zones within approximately 20–30 km downstream of promontories and peninsulas are associated with inferred softer ice. Upstream of the Bawden Ice Rise (BIR), strain rates are substantially lower and correspond to larger inferred B values. Here, the correspondence between large fractures and a simulated confluence of meltwater plumes is hypothesized to stimulate abundant accretion of marine ice, which can actually lead to ice stiffening (Khazendar and others, Reference Khazendar, Rignot and Larour2011).

Finally, for B-SW-RL, ice is inferred to be substantially softer in the mélange area that separates the Brunt Ice Shelf from the Stancomb-Wills Ice Tongue, as well as in the mélange that separates the latter from the Riiser-Larsen Ice Shelf. These areas, which contain a heterogeneous mixture of marine ice, sea ice and ice shelf debris, have previously been inferred to exhibit lower rigidity values (within a continuum mechanics model) and act to bind large ice fragments to the coast (Khazendar and others, Reference Khazendar, Rignot and Larour2009). Since the mélange is less coherent than meteoric ice advected from the ice streams, it deforms readily and corresponds to high strain rates. Additionally, prominent surface crevasses throughout B-SW are also associated with softer ice, including several transverse rifts close to the grounding line of Brunt Ice Shelf and a frontal rift separating the northeastern corner of Brunt Ice Shelf from the Stancomb-Wills Ice Tongue. Since the nominal temporal coverage of the MEaSUREs velocity data is 2009–2016, the Halloween Crack has not yet initiated (De Rydt and others, Reference De Rydt, Gudmundsson, Nagler and Wuite2019). At the southern edge of Brunt Ice Shelf at the base of Chasm 1, ice is actually inferred to have high mean B, but since uncertainties are large here (Fig. 7), we consider this to be a smoothing artifact stemming from larger thickness errors near the large rifts. Upstream of the prominent pinning point on the Riiser-Larsen Ice Shelf (PP in Fig. 6), ice is inferred to be stiffer, similar to what we observed with the pinning points for the simulated ice shelf as a compensation for unmodeled basal drag. The thinner ice downstream of the pinning point is correspondingly inferred to be softer. We do note that the orientation of the flow field relative to the pinning point is more oblique than that of our simulated shelf, which likely is the source of the more complex strain rate pattern adjacent to the pinning point (Figs S7, S11). Finally, upstream of Lyddan Island in the mélange at the eastern edge of Stancomb-Wills, ice is inferred to have high rigidity, but as this area corresponds to both low strain rates and low driving stress, the uncertainty in rigidity is very large.

4.4. Posterior predictive distributions and ice shelf buttressing

After obtaining the variational distribution that best approximates the posterior distribution for the ice rigidity, we can compute a posterior predictive distribution for any quantity or forward model that depends on the rigidity. The most straightforward way to accomplish this is to generate random samples from the variational distribution and pass each sample through the forward model of interest, i.e. Monte Carlo approximation. For example, one could perform a dynamic perturbation analysis on specific ice shelves by applying some form of stress perturbation at the ice front (calving event, gain/loss of buttressing sea ice, etc.) and running prognostic simulations for different realizations of the rigidity, sampled from the posterior distribution. This type of analysis has been performed in many studies to assess sensitivity of ice shelves to changing climate conditions (e.g. Schlegel and others, Reference Schlegel2018; Nias and others, Reference Nias, Cornford, Edwards, Gourmelen and Payne2019), but usually the rigidity field is varied by choosing some uniform upper and lower bound guided by expected temperature variations or other a priori knowledge on creep mechanisms. By instead using the posterior distribution to draw samples of the rigidity, we automatically incorporate information derived from surface observations while also allowing known physical laws (e.g. SSA equations) to induce realistic covariances between values of the rigidity over finite length scales. In other words, the combined information from data and flow equations results in more realistic samples of physical parameters consistent with all available knowledge.

Since one of the most important physical implications of ice shelf rheology is the amount of buttressing applied to inland grounded ice, we use the variational distribution for B to compute the distribution of maximum buttressing factors following Fürst and others (Reference Fürst2016). The normal buttressing number is defined as (Gudmundsson, Reference Gudmundsson2013):

(19)$$K_n = 1 - {\hat{{\bf n}} \cdot {\boldsymbol \tau} \hat{{\bf n}}\over N_0},\; $$

where τ is the deviatoric stress tensor, the quantity $N_0 = {1\over 2} \rho _i \left(1 - {\rho _i}/{\rho _w}\right)gh$ is the vertically integrated pressure exerted by the ocean (with density ρ _w) on the ice shelf, and $\hat {{\bf n}}$ is a normal vector selected to be aligned with the second principal stress, following Fürst and others (Reference Fürst2016). By performing systematic calving simulations where ice is removed from an ice shelf up to different buttressing factor isolines, the increase in ice flux across the ice front or grounding line can be predicted for various buttressing factors (Fürst and others, Reference Fürst2016). The buttressing factor above which ice flux is projected to rapidly increase then serves as a buttressing threshold for a given ice shelf. The isoline corresponding to the threshold can then delineate regions of ‘passive’ shelf ice (PSI), defined as ice that can be removed without significantly altering the flow dynamics of the adjacent ice. As the normal force in the buttressing factor is computed from the ice stress tensor, which itself depends on the rigidity B to estimate the stress components, the buttressing factor will be subject to random variations consistent with the posterior samples of B. We can therefore estimate the expected variation in PSI consistent with the surface observations. To estimate a more realistic estimate of PSI area specific to calving, we only include buttressing factor isolines that form polygons that intersect the ice front, meaning we exclude areas of isolated PSI closer to the grounding line.

The buttressing thresholds originally presented by Fürst and others (Reference Fürst2016) corresponded to flux increases across the ice front, leading to threshold values of 0.3–0.4 for the ice shelves investigated here. Alternatively, thresholds defined for increased ice flux across the grounding line are found to be a better predictor for ice shelf stability in response to instantaneous calving events (Reese and others, Reference Reese, Gudmundsson, Levermann and Winkelmann2018; Mitcham and others, Reference Mitcham, Gudmundsson and Bamber2022). These buttressing values tend to range from 0.8 to 0.9. For the purposes of comparison with the result of Fürst and others (Reference Fürst2016), we use a lower threshold of 0.4 roughly corresponding to a step increase in flux across the ice front. Due to slight biases between our inferred mean B fields and the fields estimated by Fürst and others (Reference Fürst2016), our threshold value of 0.4 is slightly higher than that used by Fürst and others (Reference Fürst2016) in order to roughly match the PSI regions in that study.

We observe variations in PSI that lie roughly within the bounds computed from ±10% variation of the mean B, following Fürst and others (Reference Fürst2016) (Fig. 8). However, we can observe additional spatial and statistical patterns beyond the simple ±10% variations. For the ice shelves that are laterally confined by embayments, there are a significant number of samples of the PSI boundary that exceed the upper and lower bounds. Over Larsen C, the PSI boundary samples are slightly skewed toward lower PSI areas. However, several posterior samples of B actually connect passive ice centered on the rift originating from GIR to passive ice at the ice front, which increases total PSI area and slightly reduces the vulnerability of Larsen C to ice loss. Over FRIS and Ross, the PSI distribution is more symmetrical, although the former has a long tail of lower PSI areas, which correspond to a slight increase in vulnerability of those shelves to ice loss. Finally, over B-SW-R, the distribution of PSI is near-symmetric and lies well within the ±10% bounds. However, the difference in spatial extent between the ±10% bounds is larger than for the other ice shelves, particularly for the Stancomb-Wills ice tongue, which indicates a greater sensitivity to variations and uncertainties in inferred ice rigidity. This sensitivity is likely reflective of the lack of lateral confinement and drag and highlights the importance of embayment geometry on ice shelf buttressing force. Overall, these results demonstrate that calibration of ice shelf rigidity and associated uncertainties using surface data can both inflate/deflate predictive uncertainties and needs to be performed on a shelf-by-shelf basis.

Fig. 8. Stochastic analysis of maximum buttressing factor for West Antarctic ice shelves, following Fürst and others (Reference Fürst2016). Background 2D buttressing fields are computed from the mean B inferred from variational inference for each ice shelf. The colormap is constructed to highlight a threshold buttressing value of 0.4, which roughly corresponds to a step increase in ice flux across the ice front for removal of ice up to the 0.4 buttressing isoline. Thus, blue areas correspond to ‘passive’ ice. The thick solid and dashed dark blue lines correspond to the 0.4 isoline for a ±10% variation of B about the mean, respectively. Thin gray lines correspond to the 0.4 isoline for B samples from the variational posterior distribution. For each ice shelf, a histogram is shown of the passive ice shelf area estimated from samples from the posterior, along with the same ±10% lower and upper bounds shown in the maps.