Edited by: Alun Hubbard, University of Tromsø, Norway and Aberystwyth University, UK
Reviewed by: Nathaniel K. Newlands, Federal Government of Canada, Canada; Stephen John Livingstone, University of Sheffield, UK; Ninglian Wang, Chinese Academy of Sciences, China; Mathieu Morlighem, University of California, Irvine, USA
*Correspondence: Douglas J. Brinkerhoff
This article was submitted to Cryospheric Sciences, a section of the journal Frontiers in Earth Science
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
We develop a Bayesian model for estimating ice thickness given sparse observations coupled with estimates of surface mass balance, surface elevation change, and surface velocity. These fields are related through mass conservation. We use the Metropolis-Hastings algorithm to sample from the posterior probability distribution of ice thickness for three cases: a synthetic mountain glacier, Storglaciären, and Jakobshavn Isbræ. Use of continuity in interpolation improves thickness estimates where relative velocity and surface mass balance errors are small, a condition difficult to maintain in regions of slow flow and surface mass balance near zero. Estimates of thickness uncertainty depend sensitively on spatial correlation. When this structure is known, we suggest a thickness measurement spacing of one to two times the correlation length to take best advantage of continuity based interpolation techniques. To determine ideal measurement spacing, the structure of spatial correlation must be better quantified.
Bed elevation is required to model glacier dynamics. Measurements are typically performed with spatially localized radar soundings (e.g., Allen et al.,
Many widely used digital elevation models (DEMs) of subglacial topography are based upon classical geostatistical techniques such as Kriging (Bamber et al.,
Inverse methods are widely used in glaciology. MacAyeal (
An important limitation exists in all of the above studies: none rigorously quantify the error bounds of their solutions in the sense that uncertainty estimates are linearized around the optimal solution or estimated empirically by comparison with independent data. Viewed in a Bayesian context, these methods report the maximum
Monte Carlo sampling techniques can provide distributions of model parameters. Several examples of the application of Monte Carlo techniques to glaciological problems exist. Petra et al. (
In this paper, we use the Metropolis-Hastings (MH) algorithm (Hastings,
We apply our method to three cases. We begin by considering a synthetically generated glacier, where the simulated velocity, mass balance, and thickness measurements (or prior estimates) are corrupted with noise, and the thickness field is recovered under several assumptions regarding covariance and measurement spacing. We then examine the degree of uncertainty induced by the choice of error structure and physical properties of the glacier. Next, we apply the method to Storglaciären, a ≈ 3 km long alpine glacier where dense measurements of thickness, velocity, and mass balance are available. This provides an interesting test case with which to gauge the uncertainties induced in topographic estimates due to real-world uncertainties. Finally, we apply the method to Jakobshavn Isbræ, the largest outlet glacier on the Greenland Ice Sheet, and assess resulting uncertainty estimates in the context of making a recommendation for flightline spacing during future airborne radar campaigns.
Bayes' theorem states that
Construction of the posterior requires two components. First, the likelihood
All inverse problems incorporate prior information in one form or another (smoothness for example), which is necessary because of the ill-posedness of such problems. One particular advantage of Bayesian methods is that these prior assumptions, which are often vacuously defined in other inverse methodologies, are defined precisely here and are subject to scrutiny. Another is that changes in assumptions and the addition of new information are easily incorporated by changing the definitions of likelihood and prior without any structural changes in the inference procedure.
We define data to be observed quantities that have not been considered when forming the prior distribution. We denote data with a hat. Data may include observations of surface speed
While we assume that the uncertainties Σ
Equation (4) is valid only at an instant in time. On the contrary, data are observed over finite time intervals. The length of these intervals may vary between observables, and are often not aligned with one another (e.g., velocity may have been observed over the entire winter of 2008, while thickness observations have been taken almost instantaneously during the summer of 2005). Observations may thus be far from the cotemporal and instantaneous values required to make Equation (4) valid. This incongruity in and between the characteristic time scales of the observation and model processes is responsible for an additional source of error that must be accounted for.
One way of constructing a valid mass conservation relationship over the entire observation period is to time integrate the forward model over the range
Monte Carlo methods such as the MH algorithm are computationally expensive because the likelihood (and thus the forward model) must be evaluated many times. Fortunately, as long as
Non-dimensionalization of Equation (10) produces an interesting result: the model depends on only a single non-dimensional parameter,
The specification of prior distributions on model parameters depends on the variable and the problem being considered. However, a universal requirement is that priors be chosen before considering any of the observations contained in
We assume no prior knowledge of the ice thickness
It is often preferable to specify a mass balance distribution in which the data have already been reanalyzed with a separate climate model or interpolated with some other method since these are more advanced than one we could include. Here we consider the estimated posterior distribution of a climate model to be the prior on mass balance, since observations have already been included. We assume that the mass balance prior is also a Gaussian process
Once again, surface velocities and depth-averaged velocities are only equal when sliding accounts for all glacier movement. The other end-member, pure deformation, bounds the value of surface velocity at
The posterior distribution has no closed form and must be characterized by sampling. With samples in hand we can evaluate their statistical properties as a proxy for the posterior distribution. Our method of choice for this procedure is the Adaptive Step Metropolis-Hastings algorithm (Hastings,
Pseudocode describing the Metropolis-Hastings algorithm. Note that
Choose initial parameter values |
|
Draw proposal sample |
|
Sample |
|
Set |
|
Set |
|
|
The MH algorithm operates by traveling through parameter space according to steps drawn from a proposal distribution, in this simple case an independent multivariate normal centered around the current point. If the posterior probability is greater at the proposed point than at the current point, then the proposed step is accepted, and the algorithm continues from the new point. If the likelihood is lower at the proposed point, then a step is taken with probability equal to the ratio of the current and proposed points. In the adaptive step variant of the MH algorithm, the width of the proposal distribution is adjusted such that an optimal proportion of proposals are accepted. The algorithm is ergodic (Hastings,
We used the finite element ice sheet model VarGlaS (Brinkerhoff and Johnson,
15 | km | Domain length | |
0 | m | Minimum elevation | |
1500 | m | Maximum linear elevation | |
200 | m | Sinusoidal topography amplitude | |
500 | m | Headwall amplitude | |
700 | m | Headwall decay length | |
σ |
50 | m | Random variability amplitude |
250 | m | Topographic correlation length | |
−80 | ma−1 | Minimum specific balance | |
10 | ma−1 | Maximum specific balance | |
3 | Shape factor | ||
250 | m | Thickness correlation length | |
1500 | m | Mass balance correlation length | |
0.01 | Relative thickness uncertainty factor | ||
0.1 | Relative speed uncertainty factor | ||
σ |
2 | ma−1 | Specific balance uncertainty |
σ |
10 | ma−1 | Minimum velocity uncertainty |
{11,6,3} | Number of data points | ||
100 | Number of grid cells |
As a first experiment, we used the MH algorithm to sample
We assumed that the error structure was known, and that the likelihood model was given by
We assumed no prior knowledge of the actual thickness values but that the thickness length scale was known, or that
Using the MH algorithm, we drew
Figure
Nonetheless, the maximum
An interesting feature of the posterior distributions of both mass balance and surface velocity is that the posterior distribution is more specific than the prior. This implies that not only are these fields contributing information to the estimation of thickness, but data and smoothness constraints on the thickness also feed back.
There are several mechanisms for assessing whether a distribution has become stationary, some heuristic and some quantitative. For a good approximation to the posterior distribution to be achieved, (a) the sampler must traverse the support of the sampled function many times, (b) the sampler must visit the entire support, and (c) the region traversed by the sampler should be insensitive to initial conditions.
Examining traces (i.e., the history of parameter values at each sampler iteration) gives a heuristic means to assess convergence. Figure
Figure
The Gelman-Rubin statistic (Gelman and Rubin,
Figure
The relative uncertainty with which the ice thickness can be recovered using mass conservation methods is a function of the relative uncertainties and covariance structures of surface velocity measurements, mass balance measurements (or assimilated model output), ice thickness measurements, and ice thickness measurement density. It is reasonable to suspect that an improvement in any of these factors would lead to a commensurate improvement in posterior predictive capabilities with respect to thickness. Nonetheless, it is not clear which of these factors contributes the lion's share of posterior variance. This information is key in forming plans of additional data acquisition; we would like to know which data improvements (and at what densities) will most improve thickness estimates.
Here we examine the influence that surface velocity uncertainty, mass balance uncertainty, and measurement density have on the relative thickness uncertainty (quantified here as
In order to assess the functional relationship between these uncertainties, we performed simulations over an array of 10 equally spaced relative velocity uncertainties between 0 and 20% of the true maximum surface velocity, 10 relative mass balance uncertainties between 0 and 50% of the true maximum mass balance, and the number of data points
Figure
Thickness precision as a function of mass balance uncertainty (Figure
Storglaciären is a ≈ 3 km long polythermal glacier in northwestern Sweden. It possesses the longest spatially distributed surface mass balance record of any glacier (Holmlund et al.,
Thickness and surface observations were adapted from Herzfeld et al. (
We estimated the flowline width by using a first order glacier flow model (Brinkerhoff and Johnson,
A prior distribution on mass balance was generated from the measured annual net balance for Storglaciären between years 2000 and 2013 (Jansson,
Thickness measurement uncertainty was assumed to be σ
We ran the Metropolis-Hastings algorithm three times for 106 iterations each, discarding the first 105 samples. Each instance was given different initial values drawn from the prior distribution. We assessed the convergence of the samples using the same methods as discussed in Section 3.3. This analysis produced similar results to those of the synthetic case (not shown).
Figure
Jakobshavn Isbræ is the largest basin (as ranked by discharge) on the Greenland ice sheet and is the fastest glacier on earth (Rignot and Kanagaratnam,
Absolute mass balance rates over the ice sheet are small compared to a mountain glacier, with a maximum accumulation of around half a meter per year, an order of magnitude lower than a typical mountain glacier. Relative measurement errors are large given identical methods, and the distances over which measurements must be extrapolated are longer. Surface velocities range over four orders of magnitude and relative errors are high in the low velocity interior of the ice sheet and low in the fast flowing outlet regions. Soundings of ice thickness are made by airborne radar in the form of discrete radar flightlines, mostly running normal to the dominant ice flow direction. These can have high observational uncertainty in topographically complex regions (Gogineni et al.,
We consider the time domain
We extracted a flowband from Jakobshavn Isbræ flowline by generating two streamlines from the InSAR derived velocity dataset of Rignot and Mouginot (
Velocity magnitudes were also derived from Rignot and Mouginot (
A simple uncertainty estimate that can account for the factors given above assumes that observations have an uncertainty of σ
Mass balance was drawn from the regional climate model RACMO2/GR (Ettema et al.,
The Greenland ice sheet is not in steady state (Motyka et al.,
We utilize ice thickness data obtained by aerial radar using the Multichannel Coherent Radar Depth Sounder (MCoRDS) (Allen et al.,
Variograms computed for IceBridge and pre-IceBridge data have placed the average range (the distance at which samples become uncorrelated) for Greenland between 58 km (Morlighem et al.,
We again ran the Metropolis-Hastings algorithm three times for 106 iterations each, discarding the first 105 samples, and each instance was given different initial values drawn from the prior distribution. We assessed the convergence of the samples using the same methods as discussed in Section 3.3, and consideration of traces, histograms drawn from different samplers, and the Gelman-Rubin statistic all indicated sampler convergence.
Figure
The algorithm makes significant adjustments to the mass balance and width functions in order to accomodate velocity observations, particularly in the ice sheet interior. The requirement that the algorithm finds a posterior mass balance distribution that is improbable with respect to the prior distribution implies that at least one of the data sets considered herein may have a mean that is far from the true value. One hypothesis is that the surface mass balance model is underestimating aeolian snow redistribution.
Figure
The method presented herein was limited to the flowline case, and we argue that this context is useful for assessing the uncertainty that we expect from using mass conservation methods. However, many modeling applications require fields over the map plane. Because Monte Carlo methods require the solution of a forward model many times, the transition from solving the balance velocity equation on a flowline, which requires only quadrature, to the map plane, which requires the solution of a linear system of equations (not to mention an inherently greater number of degrees of freedom) is inherently expensive. However, this is not to say that the problem is intractable; because the Metropolis-Hastings algorithm is subject to the Ergodic Theorem, it is efficiently parallelizable in the sense that rather than running a single instance of the algorithm for many iterations, we can run many (independent) instances of the algorithm, and concatenate the resulting samples (Murray,
The application of a geostatistical interpolation technique requires the selection of a model for the spatial covariance of the field in question, in this case ice thickness. In Bayesian methods this involves the
These two processes (covariance function definition and regularization) are equivalent: regularizing on the square of the gradient is the same as using a Gaussian covariance function. The regularization parameter thus has physical relevance in that it is proportional to correlation length, and care must be taken in the selection of a smoothness parameter. It is not appropriate to regularize away the non-uniqueness of the problem with L-curve analysis because this selects the smoothest solution for which the data retains a good fit. However, this level of smoothness may not be physically mandated, and would tend to reduce the feasible range of solutions. Furthermore, the selection of a particular norm on the thickness gradient to minimize has the effect of choosing a covariance function. Experimental variograms suggest that the appropriate model for Greenland is usually exponential, yet the 2-norm regularization commonly used in basal topography inversions implies a Gaussian covariance. This may be desirable if a large scale volume estimate is the goal (Bahr et al.,
Nonetheless, deducing the appropriate covariance model and associated parameters is a difficult task, due to the sparsity of bedrock elevation measurements. This is compounded by the fact that glacier flow tends to produce landscapes with anisotropic topographic variability. This variability is spatially heterogeneous, and ice sheet interiors may have different geomorphic properties than the ice sheet margin or heavily glaciated mountain regions. The problem should be addressed by detailed radar soundings of ice thicknesses over regions deemed geomorphically representative of large scale conditions, as well as through analysis of recently deglaciated terrain. For ice sheets, this includes understanding the topographic variability in the heretofore less observed interior regions.
The selection of bed elevation measurement spacing is of great importance for future data acquisition campaigns, both in mountain and ice sheet environments. Considering Figure
For Jakobshavn Isbræ, if we take the available ice-sheet wide empirical variograms computed from collected flightlines as guidance, then an isotropically oriented measurement spacing of approximately 10–20 km is appropriate. A more detailed look at existing measurements, particularly with an eye for discerning anisotropy in topographic covariance could improve estimates of uncertainty. This analysis could be improved further by collection of high density radar measurements over patches in order to better understand the short range topographic correlation structure that governs smoothness. For accurate assessment of the covariance structure, these patches would need to be at least 3 times the correlation length to a side. Regardless of topographic correlation structure, the admissible measurement spacing for a desired accuracy increases given more precise mass balance measurements (above a certain measurement spacing) and more precise velocities (to a point).
We have developed a Bayesian statistical model for inferring the posterior probability distribution of ice thickess given sparse observations thereof, coupled with observations of mass balance and surface velocity. Model variables are represented with Gaussian processes, which allow the specification of a covariance structure and prior information. These three parameters are linked through mass continuity. Our work advances upon previous methods in a few primary ways.
The continuity equation must be reformulated from acting instantaneously to over a finite time period. This point is subtle, because the resulting equation is similar to the original. However, the interpretation of the variables involved as time averages induces an additional step in the observation process, and we hence include an additional source of uncertainty (on top of measurement uncertainty) representing the deviation of observed quantities from their averages over the chosen averaging period of the continuity equation. Further constraining the form and magnitude of this uncertainty will be an important advance toward improving continuity-based interpolation schemes.
The Bayesian perspective allows the critical assumptions made in this model process to be elucidated. When assumptions of normality are used, they must be used explicitly. However, the method is general and assumptions of normality are not required. For example, this generality allows us to use a uniform distribution to model the relationship between surface and depth-averaged velocity.
Using Gaussian processes to represent fields such as thickness and mass balance allows for a straightforward and general way to impose smoothness constraints. We reiterate that there is a strong relationship between spatial covariance structure and regularization, and an accurate assessment of this value is critical in accurate modeling of basal topography as well as in determining the precision of those estimates.
For uncertainty estimates typical of velocities, mass balance, thickness, and data density in a mountain glacier, the algorithm reconstructed the basal topography, the distribution of which showed a relative uncertainty (as quantified by the normalized standard deviation) of around 10% at locations far from thickness measurement locations. Thickness uncertainty varies almost linearly with velocity uncertainty, due to the lack of an imposed covariance structure on velocity. On the contrary, the propogation of mass balance uncertainties is also influenced by the length scale of permissible variability in thickness and mass balance itself, as well as by data densities.
Application of the method to Storglaciären, a 3.1 km mountain glacier in the mountains of Sweden was successful, though the resulting uncertainties were relatively large. This is a result of the sparse nature of velocity measurements there, as well as the relatively large uncertainties in long-term average mass balance due to the high degree of interannual variability evident in the mass balance record.
At Jakobshavn Isbræ, the algorithm produced bed estimates with uncertainties varying greatly as a result of data density and relative uncertainty in different regions of the ice sheet. Consideration of mass conservation can improve thickness estimates in regions of low relative velocity and mass balance error, but tends to produce large uncertainties in regions of slow flow and small mass balance. We found that estimates of thickness uncertainty also depended strongly on correlation structure, which is equivalent to regularization in the PDE-constrained optimization context. Obtaining a better estimate of covariance structure would facilitate further application of mass conserving algorithms to ice sheets.
Based on our analysis, we suggest a flightline spacing of one to two times the topographic correlation length in order to best leverage mass conservation based interpolation techniques. For the Jakobshavn Isbræ region, we estimate this spacing to be between 10 and 20 km, but more work is necessary to determine the ideal value for the remainder of the ice sheet.
DB developed the method and wrote the paper. AA and MT provided critical insight into validity of model assumptions and methodology.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
DB was supported by NSF Graduate Research Fellowship grant number DGE1242789. AA was supported by NASA grants NNX13AM16G and NNX13AK27G. MT was supported by NSF grant PLR 1107491. Thanks to Jesse Johnson, Ron Barry, Regine Hock, Christina Carr, and Ed Bueler for discussions and review that led to great improvements to the manuscipt.