Glacier volume estimation as an ill-posed boundary value problem

Introduction Conclusions References


Introduction
A glacier's volume is its most fundamental geometric property, and knowledge of the total volume of all the world's glaciers and ice caps (GIC) is essential to accurate estimates of near-term sea-level rise (Meier et al., 2007;Bahr et al., 2009;Pfeffer, 2011;Radić and Hock, 2011).Unfortunately, among the hundreds of thousands of glaciers and ice caps around the world, direct knowledge of volume through measurements of ice thickness exists for only a scant few.In almost all cases, a glacier's volume (as well as subsurface geometry, internal velocities, sliding rates and other internal parameters) is inferred from surface parameters such as the surface mass balance, surface velocity, and surface geometry.In particular, estimations of GIC volume have been extrapolated from surface area using volume-area scaling (c.f.Chen and Ohmura, 1990;Bahr et al., 1997;Radić and Hock, 2010) and have been extrapolated more recently by numerical Introduction

Conclusions References
Tables Figures

Back Close
Full inversions of emergence velocities and the surface geometry (e.g.Huss and Farinotti, 2012;Clarke et al., 2013).Both volume-area scaling and numerical inversions take information about the surface boundary of a glacier and extrapolate information to the bottom boundary.For example, a typical scaling analysis uses the surface area to estimate the mean glacier thickness, from which volume (thickness times area) is estimated.In other words, scaling gives an estimate of the location of the bottom boundary.Numerical inversions for volume also estimate the position of the bottom boundary, because no matter how it is calculated, a glacier's estimated volume gives us information about the average thickness of the ice.Although volume-area scaling is not often called an inversion, it is in essence an inversion; surface data is being used to infer unseen parameters within the glacier, such as thickness (e.g., p. 4 Zhdanov, 2002;Lliboutry, 1987, p. 177).
Inversions of this type over-specify data on one boundary (the glacier surface) and under-specify the conditions on another boundary (the glacier bed).This creates an ill-posed boundary value problem (e.g.Lliboutry, 1987, p. 11 and 178).The term "illposed" does not indicate that the boundary conditions are incorrect or improperly specified.Instead, "ill-posed" means that the solution is inherently unstable, is not unique, and is sensitive to small changes in the over-specified boundary.In the case of glaciers, all of the ill-posed information is provided at one boundary (the surface) and all of the mathematically unstable results are extracted at another boundary (the bed).The instability has nothing to do with the type of numerical solution, and will exist whether or not the volume is derived from a finite element analysis, finite difference analysis, volume-area scaling analysis, or by any other approach.The instability of the solution is inevitable and cannot be improved with better models, different models, or new models.It does not matter if the solution is analytical or numerical.Unless there is additional data specified at points below the surface of the glacier, the volume solution will remain ill-posed (e.g.Courant and Hilbert, 1966;Zhdanov, 2002).
Ill-posed boundary value problems are well-studied in engineering, applied mathematics, geophysics and their associated literature.For many ill-posed geophysical Introduction

Conclusions References
Tables Figures

Back Close
Full problems, approximate rather than exact solutions are possible, though "inverse source problems" of the type discussed here will have an infinite number of possible approximations (Zhdanov, 2002, p. 4 and 18).Within glaciology, the theoretical and practical implications of ill-posed inversions have been evaluated by a number of authors including Hantz and Lliboutry (1981), Lliboutry (1987, p. 177), Balise and Raymond (1985), MacAyeal (1993), Bahr et al. (1994), Joughin et al. (2004), Truffer (2004), Chandler et al. (2006), Gudmundsson and Raymond (2008), Maxwell et al. (2008), andRaymond andGudmundsson (2009) among many others.The theoretical implications are that information is lost exponentially as a calculation progresses deeper into a glacier and becomes increasingly unstable (Bahr et al., 1994;Zhdanov, 2002, p. 26).Therefore, appropriate control of the calculation can damp (though not eliminate) the instability and can lead to reasonable and approximate (though non-unique and less precise) solutions.
For example, in an approach presented by Huss and Farinotti (2012), a numerical inversion uses surface topography to extract the basal shear stress.In turn, this allows an iterative solution for the glacier thickness.When the thickness is calculated at many locations, the glacier volume can be estimated.Therefore, the glacier volume solution depends entirely on ill-posed estimations of the basal shear stress derived from surface parameters and is thus non-unique.Similarly, scaling techniques (e.g.Bahr et al., 1997) use the surface area (top boundary) to infer the average thickness (bottom boundary).The resulting volume estimate is neither unique nor necessarily precise.Analytical methods are no more immune to the ill-posed instability than numerical methods.As pointed out by Lliboutry (1987, p. 180), the best that we can do is say that a particular solution has been derived by a model or a calculation; we cannot claim that the sole solution has been derived.Only with the application of special filtering, a Bayesian analysis, or other techniques (e.g., Raymond and Gudmundsson, 2009;Zhdanov, 2002;Bahr et al., 1994) will the model give a reasonable (though non-unique) approximation to the exact solution.Introduction

Conclusions References
Tables Figures

Back Close
Full While a volume solution derived from surface data will always be ill-posed, there are three important mitigating factors that make both scaling and numerical inversions useful and practical tools.First, ill-posed does not mean unsolvable.A great deal of progress has been made by turning ill-posed problems in seismology into a set of equivalent well-posed problems (Zhdanov, 2002).Similar progress has been made with glacier inversions (e.g., Maxwell et al., 2008;Raymond and Gudmundsson, 2009), and even though these and future improvements will always be an approximation, illposed inversions are best treated with careful calculation and interpretation rather than assumed intractability.It is not our intent to imply that either numerical models or scaling are fundamentally intractable; both can work well within their limitations.
As a second mitigating factor, the random errors in volume will differ for each glacier.Therefore, by the law of large numbers, the sum of many glaciers' volumes will give a reasonable estimation of their total volume.For example, in many sea-level rise applications, only the aggregate volume of the world's glaciers is important (rather than the volume of each individual glacier).Note, however, that the variance of the sum can be very large.As shown in the following section, these variances can grow exponentially large making the total volume estimate much less reliable.Furthermore, each method for estimating glacier volume will have additional sources of error that are separate from the ill-posed instability, such as errors in data, errors from numerical instabilities, and errors in specified parameters (like the scaling constant, or flow law exponent, or sliding law).The impact of these additional errors on the variance of the total volume must be evaluated separately from the ill-posed boundary value problem.However, unless special precautions are taken (as described in the following section), it is the exponentially large variance from the ill-posed instability that has the greatest potential to overwhelm an estimate of total GIC volume.
As a third mitigating factor, the ill-posed instability is dependent on the spatial wavelength or resolution of the solution (e.g.Bahr et al., 1994;Chandler et al., 2006;Zhdanov, 2002, p. 30-31).High resolution model results will be less accurate than low resolution results (where the highest possible resolution is the discretization or sampling interval, and the smallest resolvable spatial wavelength is twice the resolution); for example, a model with 100 m grid cells will have smaller errors than a model interpreted at 10 m intervals.In particular, Bahr et al. (1994) showed that errors from the ill-posed instability grow exponentially as a function of decreasing spatial wavelength.
If H is the glacier thickness and λ is the spatial wavelength of any glacier parameter (velocity, stress, etc.), then where E s is the measurement error of the surface parameter, and E b is the resulting error at the bed of the glacier.For example, E s and E b can represent the surface and basal errors in velocity or stress.The value of r(n) is a constant that depends on Glen's constitutive law parameter n as where Re [. . .] gives the real component of the complex number (Bahr et al., 1994).
For the typical value of n = 3, the function gives r(n) ≈ 0.6, but regardless of the exact value, the growth rate is always exponential.Even at wavelengths of half an ice thickness the errors at the bed will be huge with E b /E s = 1881 when n = 3; in other words basal errors will exceed surface errors by a factor of almost 2000 or by over three orders of magnitude.For reasonable errors that grow by no more than an order of magnitude, one ice thickness is the shortest practical wavelength at the bed (E b /E s = 43).The exponential form of Eq. ( 1) follows from the generic form of the stress equilibrium and strain rate compatibility equations (Bahr et al., 1994), and can be derived exactly without simplifications or assumptions when n = 1.Many other geophysical inverse problems have a similar exponential form for errors (Zhdanov, 2002).In an interesting Monte-carlo inversion, Chandler et al. (2006) have suggested that Eq. ( 1) should

TCD Introduction Conclusions References
Tables Figures

Back Close
Full have a power-law rather than an exponential form.If this is the case, then the derivations that follow in this paper will be conservative, and practical solutions at the bed of a glacier could include some shorter spatial wavelengths than indicated here.However, although postulated, no analytical power-law formulation for error growth has been derived, so we will use the exponential form which is guaranteed to give a reasonable and conservative lower bound for the limiting short wavelengths.
In practice, short spatial wavelengths can be removed with a low-pass frequencydomain filter (apply a Fourier transform to the thickness solution and then truncate high frequencies in the frequency domain).If we assume smooth basal topography, for example, then the resulting low-resolution, long wavelength solution will be very accurate.
In reality, such an assumption of smooth basal topography may not be realistic for an arbitrary glacier.Many glaciers have ice falls, over-deepened valleys, barely submerged nunataks, and other high spatial frequency features which can be buried under ice and cannot be known a priori.The sum of many short (spatial-domain) wavelengths can have a surprisingly large macroscopic effect on the spatial distribution of ice thickness.
For example, consider a glacier with an ice fall created from basal topography that is approximated with a Heaviside step function.The Fourier transform of a step function has many high frequency components that cannot be removed (filtered) without turning the ice fall into a shallow ramp.Conversely, a glacier with a shallow ramp can be inadvertently turned into an ice fall by the instabilities at short spatial wavelengths.
The difference between a ramp and a step function in the bed topography can have a significant influence on the glacier's volume.
To avoid filtering explicitly, other numerical techniques have differing strategies for removing the unphysical short wavelength oscillations.For example, an inversion with Bayesian inference assigns probabilities to each solution, so that the short wavelength errors become a possible but improbable solution (Raymond and Gudmundsson, 2009).Regardless of the particular technique (filtering, Bayesian, or other), the improbable high amplitude short wavelength errors must be identified, minimized, and/or removed.In this paper, we do not develop a new technique for inversions, nor do we derive a new technique for error supression or removing short wavelength errors.Instead, we assume and use existing theory in the form of Eq. (1) to discuss glacier volume estimations as an ill-posed inversion.
As the following analysis shows, volume-area scaling inherently averages over long wavelengths and does not over-interpret the ill-posed results at shorter wavelengths.
In contrast, numerical inversion with fixed grid-spacing can over-interpret the volume estimation by including excessively short wavelengths.As a result, the published total GIC volume estimations from volume-area scaling (e.g.Radić and Hock, 2010) do a good job of controlling ill-posed instabilities and will be inherently more precise (less variance) than the total volume estimates from published numerical models that do not filter or otherwise minimize short wavelength errors in the manner(s) described above.However, both methods of volume estimation are invaluable and essential to building an ensemble solution; therefore the following section also details the appropriate wavelengths for filtering this new class of numerical inversions (c.f.Farinotti et al., 2009;Huss and Farinotti, 2012;Clarke et al., 2013).Many published inversions already acknowledge the need for averaging, filters, or error suppression of some kind (e.g., Huss and Farinotti, 2012;Clarke et al., 2013), but we argue that with ill-posed errors, additional precision and careful reasoning are essential for accurate GIC volume estimates.

Error analysis for ill-posed instability
To avoid ill-posed errors, the wavelength λ must be selected so that the right hand side of Eq. ( 1) grows as either a constant or as an even slower function of the glacier size.In other words, ideally the exponent in Eq. ( 1) should be 0, but practically this means keeping its value much less than 1:

TCD Introduction Conclusions References
Tables Figures

Back Close
Full For n = 3, the constant terms become 2πr(n) = 3.77, and the shortest acceptable wavelength will be approximately four times the glacier thickness.For example, if a model estimates the velocity at the bed as a function of position, then the computed basal velocity should be filtered to remove all wavelengths shorter than four times the thickness, and preferably filtered at much longer wavelengths as well.As a result, a low-resolution estimate of the entire glacier volume is generally manageable, but an estimate with a spatial resolution of better than four ice-thicknesses is not.Consider the special case where a glacier's volume is estimated from volume-area scaling.The thickness can be replaced by a function of the glacier surface area S. Using scaling relationships, where V is volume, γ = 1.375 for glaciers, γ = 1.25 for ice caps, and c = 0.034 is a scaling parameter (Bahr et al., 1997;Bahr, 1997a).Similarly, the area can be expressed as a function of the glacier length L (Bahr, 1997b).

S = bL q+1
(5) where b ≈ 1 is a scaling constant and q = 0.6 for glaciers and q = 1 for ice caps.Solving for L gives the characteristic glacier length associated with volume-area scaling.
In other words, for volume estimates that are derived from volume-area scaling, the wavelength is λ = 2L.From Eq. ( 5), Substituting Eqs. ( 4) and ( 6) into (1), the errors from the ill-posed instability will grow as TCD Introduction

Conclusions References
Tables Figures

Back Close
Full Note that γ − 1 − (1/(q + 1)) = −0.25 for both glaciers and ice caps.Also, e πr(n)cb 1/(q+1) = e 0.015 ≈ 1.0 (8) The parameters b and c have a distribution of values (Bahr, 1997a) and can vary significantly from glacier to glacier.However, even with order of magnitude estimates of b and c, Eq. ( 8) remains near 1.0, and Therefore, for volume-area scaling, the ill-posed instabilities decrease with increasing glacier and ice cap size.In fact, for glaciers larger than 1 km 2 , the exponential rapidly approaches 1, and E b ≈ E s .Appropriately short wavelengths are automatically filtered.
No additional averaging or filtering is necessary to control the ill-posed instability when using volume-area scaling.
The same filtering of short wavelengths should be done for any numerical inversion.Consider for example, a hypothetical numerical model with a grid spacing of dx along the surface of a glacier.The wavelength is λ = 2dx, and from Eq. ( 1) If dx is fixed by the model, then the errors in Eq. ( 10) grow exponentially with glacier thickness and size.If dx is selected to appropriately filter the world's largest glaciers, then that same grid size will be too large to model the world's smallest glaciers.Instead, the model should use a dx (grid size) tuned to match Eq. ( 3) according to each glacier's size.
As an example, Huss and Farinotti (2012) use a fixed grid spacing of dz = 10 m elevation bands.With a surface slope α this corresponds to a grid spacing of dx = dz/ sin α along the surface of each glacier.In that case, errors grow as Or equivalently, H dz 1.9 sin α (13) With dz fixed at 10 m, Eq. ( 13) is possible but unlikely.For typical surface slopes of 1 to 4 degrees, the glacier thickness must be substantially less than 300 m to 75 m respectively.
For example, consider a large glacier with an area of 1000 km 2 .From volume-area scaling (Eq.4) we know that the corresponding predicted volume will be roughly 453 km 3 and the average thickness will be roughly 453 m.For a slope of 1 degree, Eq. ( 13) shows that the grid spacing dz should be substantially greater than 15 m rather than the 10 m used in Huss and Farinotti's model.In addition, for any given surface slope, the errors will become exponentially worse for progressively larger glaciers.For a very large glacier of 10 000 km 2 at the same slope of 1 degree, the volume will be 10 751 km 3 , the average thickness will be 1075 km, and the grid spacing dz should be much greater than 36 m.With a substantial fraction of the global GIC volume contained in the largest glaciers, the greater ill-posed instabilities associated with these large ice masses could be significant.
Using fixed elevation bands of dz is substantially better than using a fixed horizontal grid spacing dx, but using a low-pass filter tuned to Eq. ( 3) is necessary to reduce errors.As an approximation of the potential volume errors, consider the previous scenario with a 10 000 km 2 glacier and a 1 degree surface slope.In other words, at the highest resolution, errors from the ill-posed instability could be up to roughly 32 % of the actual thickness.From Eq. ( 4), that implies potential volume errors of up to 32 % for large glaciers.
In defense of Huss and Farinotti (2012), we do not expect their model to generate ill-posed errors nearly so large.Their model averages stresses over 10 to 20 ice thicknesses during each iteration (c.f.Kamb and Echelmeyer, 1986), and they indicate that the final ice thickness distribution is smoothed by an unspecified amount.Averaging and smoothing will act as low-pass filters, though with a non-linear constitutive relationship, the averaging of stress (in each model iteration) may not propagate equivalently to all other relevant parameters like thickness (as suggested for example in the disconnect between velocity and topography in numerical inversions by Gudmundsson and Raymond, 2008).This could account for some of the 30 % root-mean square mismatch between their modeled thicknesses and the observed glacier thicknesses.Proper low-pass filtering or Bayesian statistics (Raymond and Gudmundsson, 2009) or other well-reasoned error supression (e.g., Maxwell et al., 2008;Zhdanov, 2002) should give more defensible results.
For any analytical or numerical technique that does not control the ill-posed instability, the errors will grow exponentially with the size of the glacier.This is mitigated somewhat when calculating the total world-wide GIC volume -with enough glaciers, the sum of many random errors will trend toward zero.However, glaciers have a power-law distribution of sizes (Bahr and Radić, 2012).When summing the many glaciers in a small size bin, their total volume will control the random error effectively.Unfortunately, there are very few glaciers in the largest bins, and their errors will be exponentially larger than the errors of the small glaciers.The sum of these very few large glaciers may contain a sizable fraction of the total GIC volume (Bahr and Radić, 2012), but the law of large numbers will not apply, and their contribution to the total GIC volume error will be large.Therefore, if the ill-posed instability is not carefully controlled, the total GIC volume will have extremely large variances dominated by the exponential errors of the few largest glaciers.

Conclusions
Errors from the ill-posed boundary value instability are by far the most significant source of error in glacier volume estimation derived from surface observations.Other significant sources of error exist, but unless proper spatial filtering (or another equivalent error supression technique) is applied, the ill-posed boundary value instability can grow exponentially and quickly swamp all other errors.Before considering any other source of errors (from data, from model parameters, or from numerical instabilities), the ill-posed instability must be controlled by filtering all wavelengths that are less than at least four times the glacier thickness.Filtering at longer wavelengths will damp the ill-posed instability further.
Models that do a better job of filtering small wavelengths will have smaller variances and greater accuracy.Volume-area scaling intrinsically eliminates ill-posed errors by filtering at the correct wavelength.Recent GIC numerical inversion techniques are a welcome and very important alternative to scaling methods and have distinct advantages over power-law scaling; for example, power-law scaling is best used to derive aggregate volume for statistically large samples of glaciers, while numerical inversions like that of Huss and Farinotti (2012) can be applied with greater accuracy to individual glaciers.Nevertheless, numerical inversions must explicitly filter at appropriately short wavelengths to avoid exponential growth of errors; doing so would be relatively easy and will add a powerful new technique to glacier volume estimations.
When calculating total world-wide GIC volume, the summation over all glaciers reduces any source of random error, including those from the ill-posed instability.However, the variance of such estimates is still dependent on the magnitude of the original errors, particularly for large glaciers.Techniques that control the large ill-posed instability will have an exponentially smaller variance associated with the total GIC volume.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Huss and Farinotti estimate surface DEM errors on the order of E s = 10 m.At a wavelength of 20 m (twice the model grid spacing of dz = 10 m), the thickness error in meters will be Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |