Prediction under uncertainty : from models for marine-terminating glaciers to Bayesian computation

Abstract

The polar ice sheets have enormous potential impact on future global mean sea level rise. Recent observations suggest they are losing mass to the ocean at an accelerated rate. Skillful prediction of the ice sheets' future mass loss remains difficult, however; observations of key variables are insufficient and physical processes are poorly understood. Even when a relatively accurate dynamical model is available, computational limitations make it difficult to characterize uncertainties associated with the model's predictions. To address this prediction challenge, this thesis presents complementary developments in glaciology and in Bayesian computation.

In particular, (i) we develop new models of marine-terminating glaciers whose dynamics are controlled by an extended set of physical processes and geometric constraints; and (ii) we develop new sampling algorithms to efficiently characterize selected marginals of a high-dimensional probability distribution describing uncertain parameters. The latter algorithms have broader utility in Bayesian modeling and inference with computationally intensive models. We begin by studying laterally confined ice streams that terminate in the ocean, where they may form floating ice shelves. Such marine-terminating outlet glaciers are the main conduits by which Greenland and Antarctica drain their ice mass into the ocean. Ice shelves play an important role in buttressing the grounded inland ice. The seaward ice flow is typically accompanied by acceleration and thinning. Increased thinning eventually leads to flotation of the ice supported by buoyant forces from the ocean.

The transition region from grounded to floating ice is referred to as the grounding line (or zone), and the mass transport across the grounding line as the output flux. Previous work by Weertman (1974) and Schoof (2007) considers laterally unconfined ice streams, showing that their output flux is a monotonically increasing function of the bedrock rock depth at the grounding line. This scenario leads to the marine ice sheet instability (MISI): retreating into deeper water increases the output flux, and retreat accelerates. Therefore, stable steady states cannot exist on downward sloping beds. We extend this analysis to laterally confined glaciers and investigate when side-wall drag is sufficient to stabilize glaciers on downward sloping beds. Additionally, we include a parameterization of sub-shelf melt. We find that, whereas lateral drag can stabilize glaciers that would otherwise be subject to the MISI, sub-shelf melt can destabilize them.

Our ultimate goal is to predict future ice sheet volume and to quantify its uncertainty. We do so in the Bayesian statistical setting, conditioning our prediction on available observations. Yet characterizing a posterior distribution-using, for example, Markov chain Monte Carlo (MCMC)-involves repeated evaluations of an ice stream model, which are prohibitively expensive. Furthermore, the model parameters that need to be inferred are high dimensional, even though we are primarily interested in a low dimensional quantity: the future ice volume. We address this computational challenge by developing new structure-exploiting Monte Carlo methods that combine marginalization with surrogate modeling. Given a high-dimensional (posterior) distribution on the model parameters, whose density evaluations are computationally intensive, we construct an MCMC chain that directly targets a particular low-dimensional marginal of interest. In general, the marginal density is not available analytically.

Instead, we can compute unbiased noisy estimates of this density. Our MCMC algorithm incrementally constructs a local regression approximation of the target marginal density using these estimates. Continual refinement of the approximation, as MCMC sampling proceeds, leads to an asymptotically exact characterization of the desired marginal distribution. Analysis of the bias-variance tradeoff guides an ideal refinement strategy that balances the decay rates of different components of the error. Our approach exploits regularity in the marginal density to significantly reduce computational expense relative to both full-dimensional and pseudo-marginal MCMC.