RESEARCH NOTE

A note on digital elevation model smoothing and driving stresses

Felicity S. McCormack¹, Jason L. Roberts^2,3, Lenneke M. Jong^2,3, Duncan A. Young⁴ & Lucas H. Beem⁴

¹Institute for Marine and Antarctic Studies, University of Tasmania, Hobart, Tasmania, Australia;

²Australian Antarctic Division, Kingston, Tasmania, Australia;

³Antarctic Climate and Ecosystems Cooperative Research Centre, University of Tasmania, Hobart, Tasmania, Australia;

⁴Institute for Geophysics, Jackson School of Geosciences, University of Texas at Austin, Austin, TX, USA

Abstract

Ice-flow fields, including the driving stress, provide important information on the current state and evolution of Antarctic and Greenland ice-sheet dynamics. However, computation of flow fields from continent-scale DEMs requires the use of smoothing functions and scales, the choice of which can be ad hoc. This study evaluates smoothing functions and scales for robust calculations of driving stress from Antarctic DEMs. Our approach compares a variety of filters and scales for their capacity to minimize the residual between predicted and observed flow direction fields. We find that a spatially varying triangular filter with a width of 8–10 ice thicknesses provides the closest match between the observed and predicted flow direction fields. We use the predicted flow direction fields to highlight artefacts in observed Antarctic velocities, demonstrating that comparison of multiple observational data sets has utility for quality control of continent-scale data sets.

Keywords
Ice-flow direction; smoothing filter; Antarctic ice sheet; ice-sheet dynamics

Abbreviations
DEM: digital elevation models; MEaSUREs v2: Making Earth System Data Records for Use in Research Environments, InSAR-Based Antarctica Ice Velocity Map, Version 2, US National Snow and Ice Data Center (Mouginot et al. 2012; Rignot et al. 2011, 2017); RMS: root mean square

Introduction

Quantifying ice flow is necessary for understanding the dynamics, mass balance and evolution of ice sheets and glaciers. In a steady-state ice sheet or glacier, the net accumulation is equal to the outwards flow distribution, and the ice thickness does not vary in time at any point. The velocity field of the ice when in balance is termed “balance velocity” (Budd & Warner 1996). In the absence of observed velocities, balance velocities have been relied upon for deriving empirical relations for ice flow or sliding properties (e.g., McInnes & Budd 1984), validation of ice-sheet models (e.g., Payne 1999) or fieldwork planning. The calculation of balance velocities requires knowledge of the flow direction, which basic flow theory (Nye 1952, 1957) equates with the direction of steepest surface slope, according to the relation for the gravitational driving stress τ_d:

where r is ice density, g gravitational acceleration, H ice thickness and S surface elevation. This theory applies only in regions where surface and bed elevation gradients are small and for horizontal scales of at least the ice thickness (localized features of the surface slope do not affect the flow).

For horizontal scales between 1 and 20 times the ice thickness, longitudinal coupling of longitudinal stress gradients can substantially modify the flow field (Kamb & Echelmeyer 1986). For Eqn. 1 to apply in this instance, thickness and surface slope fields must be smoothed. Kamb & Echelmeyer (1986) demonstrated that an exponential smoothing function appropriately accounts for longitudinal stresses in Eqn. 1. Theory predicts that the horizontal scale of their smoothing function should be between one and three times the ice thickness for temperate valley glaciers and between four and 10 times the ice thickness for ice sheets (although factors such as bed slipperiness can extend these values, e.g., Gudmundsson [2003]). Le Brocq et al. (2006) compared the effect of different smoothing techniques on flux distribution calculations for the West Antarctic Ice Sheet, finding a variety of smoothing scales that produced feasible outcomes. This result suggests the need for an objective means of selecting a smoothing scale that also minimizes some measure of error when calculating the driving stress.

The aim of this study is to determine optimal and simple smoothing functions and scales for continent-wide DEMs of ice-surface elevation for the Antarctic ice sheet. We consider a variety of filters and examine the residuals between predicted and observed flow direction fields, contrasting the errors associated with using fixed versus spatially varying filter widths. We demonstrate that our technique can facilitate quality control for existing observational data sets of large-scale velocities.

Methods

To determine optimal and simple smoothing functions for continent-wide DEMs, we compared flow directions calculated from smoothed ice-surface DEMs with those of an observed velocity field. Specifically, we considered two ice-surface DEMs: (1) one presented by Bamber et al. (2009), which combines laser altimeter measurements from the ICESat-1 mission with satellite radar altimeter data from the ERS-1 satellite mission; and (2) one developed by Helm et al. (2014), which utilizes data from the CryoSat-2 data acquisition. We used observed surface velocities from the MEaSUREs v2 data set (Rignot et al. 2011, 2017; Mouginot et al. 2012).

We considered three smoothing filters: a Gaussian filter, a triangular filter proposed by Kamb & Echelmeyer (1986) and a bilinear filter. In their analysis of longitudinal stresses, Kamb & Echelmeyer (1986) demonstrate that the triangular filter is an accurate alternative to the more computationally costly exponential filter; rather than replicate their results, we chose to compare only the triangular filter with the Gaussian and bilinear filters here.

The variable width Gaussian filter is defined by Eqn. 2:

where full width at half maximum σ is given by a multiple of the ice thickness, H, and the distances from the current point in both x and y directions.

The triangular filter is defined by Eqn. 3:

where again the width σ is scaled by the ice thickness.

The local bilinear surface (z as a function of x and y) is defined as follows:

The coefficients a₁ and a₂ are the surface slopes of the DEMs in the x and y directions, respectively, and each of the coefficients a₀–a₃ were chosen to minimize the RMS mismatch between the DEM and the bilinear surface over a moving subwindow centred on the point of interest.

To derive the flow directions from the DEMs, we first applied each of the preceding filters to smooth the DEMs, then calculated the local down-slope directions (q) using the following equation:

Again, a₁ and a₂ in Eqn. 5 are the surface slopes of the DEMs in the x and y directions, respectively. When the local bilinear filter in Eqn. 4 was used, the a₁ and a₂ surface slopes were derived in the same manner as already described; when any of the remaining filters were used, the a₁ and a₂ surface slopes were calculated directly from the filtered DEM. We calculated the uncertainty weighted difference and RMS error between the flow directions predicted from the DEMs, to which each of the filters was applied, and the flow directions of the observed velocity field over the entire ice sheet. We hereafter refer to the difference between the observed and predicted fields as the “bias.” The uncertainties in flow direction for MEaSUREs v2 can be large, for example, the average absolute error in observed flow directions is 33°. To take into account the impact of flow direction uncertainty on the difference statistics, we weighted the bias and RMS error calculations by the uncertainty in the observed flow direction at each point, as follows:

Here, the total uncertainty (a) was calculated using the error estimates in both the x and y directions from MEaSUREs v2 and takes values in the range (0,180)°. The resulting scheme in Eqn. 6 assigns zero weighting to all points with uncertainties larger than 90° and a weighting of one to all points with zero associated uncertainty, linearly weighting between these values. Although the cut-off value of 90° is somewhat arbitrary, uncertainties in flow directions in excess of this threshold are sufficiently large as to warrant exclusion in our comparison between observed and predicted flow fields.

The choice of weighting scheme will influence the magnitude of the biases between the observed and predicted flow directions and where they occur. However, weighting by uncertainties in the observed flow field is particularly useful for validation of driving stresses and balance velocities calculated from DEMs, because the largest weighted biases will be a direct result of the methods used to predict the flow directions and accuracies in the DEMs, rather than a result of the underlying validation data set.

Because of possible errors in the reference surface velocity directions, we were more concerned with minimizing the weighted RMS error than the bias. We considered two cases for the width of the filter: (1) a spatially invariant filter width, calculating biases for the same range of scaling lengths but with no dependence on the ice thickness; and (2) a filter width varying spatially as a function of ice thickness—specifically, directly proportional to ice thickness for a range of smoothing scales.

Results

Weighted RMS errors and biases between the predicted and observed flow directions for each of the Bamber et al. (2009) and Helm et al. (2014) DEMs and filters are compared in Fig. 1. RMS errors ranged between 35° and 60°. In all cases, the weighted biases were approximately 2°–3° and positive, irrespective of type of filter, filter width or DEM used.

Fig. 1  Uncertainty weighted RMS errors and biases for the predicted flow angle (°) between each of the (a, c) Bamber et al. (2009) and (b, d) Helm et al. (2014) DEMs and the observed flow angle from Rignot et al. (2017). In each panel, red curves show the results for a bilinear filter, blue curves for a Gaussian filter and green curves for the triangular filter of Kamb & Echelmeyer (1986). In panels (a) and (b), spatially constant filters are used; in panels (c) and (d) the width of the filters is set according to a scaling function of ice thickness. Solid curves are the biases (°) and dashed curves are the RMS errors.

The patterns in RMS errors and biases calculated using the triangular filter were similar for both DEMs. For a spatially invariant smoothing scale (Fig. 1a, b), the RMS errors monotonically decreased with increasing filter width; for a filter width that varied as function of ice thickness (Fig. 1c, d), the RMS errors decreased rapidly for smaller spatial scales to an optimal filter width of 8–10 ice thicknesses. The weighted biases showed a slight increase with increasing filter width for both DEMs. The bilinear filter produced similar results to the triangular case for each DEM and smoothing scale but with RMS errors on average 0.5° greater. However, the bilinear filter came at greater computational expense than the triangular and Gaussian filters, requiring more than an order of magnitude greater compute time.

The patterns in RMS errors and biases calculated using the Gaussian filter differed markedly from those of the triangular and bilinear filters. For spatially invariant filter widths (Fig. 1a, b), the Gaussian filter RMS errors were minimized for filter widths of 8000–10 000 m with the Bamber et al. (2009) DEM and 4000–6000 m with the Helm et al. (2014) DEM. For Gaussian filter widths that varied as a function of the ice thickness (Fig. 1c, d), the RMS errors were minimized for two to three ice thicknesses with both DEMs, and the RMS errors increased beyond this width.

Uncertainty weighted biases between the flow directions predicted from the Helm et al. (2014) DEM and the Rignot et al. (2017) MEaSUREs v2 velocities using a triangular filter of eight ice thicknesses are illustrated in Fig. 2a. The biases, ranging between –90° and 90°, were present almost everywhere across the Antarctic continent, although positive biases were more spatially widespread. Lower absolute weighted biases arose at latitudes greater than 82.5°S and in some areas of slow flow, for example, along drainage divides in Queen Maud Land, East Antarctica, where the predicted flow directions in MEaSUREs v2 were more uncertain. By contrast, larger absolute weighted biases were apparent towards regions of streaming flow and especially marked in glaciers with channelized underlying bed topography features (e.g., Denman Glacier, Lambert/Fisher/Mellor glaciers, Thwaites Glacier, Pine Island Glacier, Foundation Ice Stream and Recovery Glacier; Figs. 2a, 3). The absolute weighted biases were also large—upwards of |60|°—in areas with strong contrasts in bed topography, steep bed slopes or deep topographies (<–1000 m), such as south of Law Dome in the Aurora Subglacial Basin (consistent with Roberts et al. 2011), although these bed features are not necessarily clear from the Bedmap2 bed elevation (Fig. 3; Fretwell et al. 2013). Large absolute weighted biases were associated with both along (azimuth) and across (range) track directions of the satellite swaths incorporated into MEaSUREs v2 (e.g., inland of the Denman, Frost and Cook glaciers in Wilkes Subglacial Basin; Aurora Subglacial Basin; and Princess Elizabeth Land).

Fig. 2  (a) Uncertainty weighted and (b) absolute differences between the flow direction (°) calculated from a triangular filter with a spatially-varying width (eight ice thicknesses) CryoSat-2 DEM by Helm et al. (2014) and the surface flow direction from MEaSUREs v2 (Rignot et al. 2017). The drainage basin divides defined by Zwally et al. (2012) are overlaid in grey.

Fig. 3  Bedmap2 bed elevation (m) (Fretwell et al. 2013).

The weighted biases were influenced by two factors: (1) the level of agreement between the predicted and observed flow directions (Fig. 2b) and (2) the level of uncertainty in the MEaSUREs v2 velocities (Fig. 4). For example, very low absolute biases at latitudes greater than 82.5°S and along drainage divides were a result of large uncertainties in the observed flow direction (Fig. 4), which outweighed large biases between the observed and predicted flow directions in these regions. Large absolute weighted biases occurred in regions of poorer satellite coverage, such as in the approach southwards in the range 80°–82.5°S, where there were larger absolute biases between the observed and predicted flow directions and relatively low uncertainties in the observed flow directions.

Fig. 4  Uncertainties in flow direction (°) calculated from the MEaSUREs v2 dataset (Rignot et al. 2011, 2017).

Discussion and concluding remarks

Our study examined the impact of filter type and smoothing scale on the biases between predicted and observed flow direction fields for the Antarctic continent. We found that Kamb & Echelmeyer’s (1986) triangular filter was most appropriate for minimizing the RMS errors and biases between observed and predicted flow direction fields (Fig. 1), regardless of whether the DEM by Bamber et al. (2009) or by Helm et al. (2014) was used to predict flow direction.

Fixed-width filters produced similar results to filters that varied as a function of ice thickness for each of the filters. For the triangular filter, the RMS errors decreased more rapidly with a varying-size filter. Given minimal increases in both code complexity and computational expense with a varying-size filter, and reduced sensitivity to smoothing scale, we see no benefit in using a fixed-size filter for this application.

Our results suggested that the Kamb & Echelmeyer (1986) triangular filter of 8–10 ice thicknesses and the Gaussian filter of 2–4 ice thicknesses performed similarly in minimizing the RMS errors. However, in their development of a method to compute balance velocities in ice sheets—which relies on accurate driving stress calculations— Brinkerhoff & Johnson (2015) demonstrated that filter widths greater than four ice thicknesses were required to minimize over-convergent flow. Concordant with Kamb & Echelmeyer (1986), the Brinkerhoff & Johnson (2015) findings suggest that a filter less than four ice thicknesses wide may not be optimal in regions where longitudinal stresses control ice-sheet flow. Applying this to our results suggests that the triangular filter of 8–10 ice thicknesses is optimal for achieving driving stress calculations from DEMs that agree best with velocity observations.

This study has highlighted implications for the calculation of balance velocities from DEMs. Using a triangular filter of width 8–10 thicknesses, we compared biases in observed and predicted flow directions. Weighting the bias by the uncertainties in the observed flow direction allowed us to attribute large biases to the methods used to calculate flow directions and/or the accuracy of the underlying DEMs. Our method assumes that ice flows along the direction of steepest surface slope, which may not necessarily be the case. However, we found large weighted biases towards fast-flowing regions (e.g., Pine Island, Thwaites, Lambert/Fisher/Mellor and Totten glaciers) where bed topographies were highly channelized, as well as in areas with strong contrasts in bed topography, steep bed slopes or deep topography (<–1000 m). Balance velocities calculated from DEMs using this method may also differ markedly from observed velocities in such regions. This finding has potential consequences for applications using balance velocities, for example, mass balance estimates.

Our results also highlighted a potential bias in the observed velocity data. Figure 1 showed a positive bias in predicted flow direction for both the Bamber et al. (2009) and Helm et al. (2014) DEMs when each of the triangular, Gaussian and bilinear filters were used. The bias was slightly lower for the Bamber et al. (2009) data set, which may be a consequence of the Rignot et al. (2017) MEaSUREs v2 velocities being calibrated against Bamber et al. (2009) rather than Helm et al. (2014). However, Young et al. (2017) compared the Bamber et al. (2009) and Helm et al. (2014) DEMs with the laser altimetry system from the International Collaboration for Exploration of the Cryosphere through Aerogeophysical Profiling II project, finding that the Helm et al. (2014) data set was slightly more accurate than the Bamber et al. (2009) data set at Dome C. Even after weighting the biases by uncertainties in the velocity observations (Fig. 4), large absolute biases between the Helm et al. (2014) DEM and Rignot et al. (2017) were correlated with satellite swaths (Fig. 2), suggesting potential artefacts in MEaSUREs v2.

We have determined that a triangular filter with an optimal width of 8–10 ice thicknesses minimizes the residual between observed flow direction fields and those predicted from continent-wide DEMs of ice-surface elevation for the Antarctic ice sheet. Such a filter can be used to improve the accuracy of driving stress—and hence balance velocity (Brinkerhoff & Johnson 2015)—calculations by ensuring self-consistent smoothing of continent-wide DEMs. Furthermore, by facilitating the comparison of ice dynamic fields between independent data sets, our technique has applications in identifying potential artefacts in observed fields and may therefore serve as a means for quality control.

References

Bamber J., Gomez-Dans J. & Griggs J. 2009. A new 1 km digital elevation model of the Antarctic derived from combined satellite radar and laser data—part 1: data and methods. The Cryosphere 3, 101–111, https://doi.org/10.5194/tc-3-101-2009.

Brinkerhoff D. & Johnson J. 2015. A stabilized finite element method for calculating balance velocities in ice sheets. Geoscientific Model Development 8, 1275–1283, https://doi.org/10.5194/gmd-8-1275-2015.

Budd W. & Warner R. 1996. A computer scheme for rapid calculations of balance–flux distributions. Annals of Glaciology 23, 21–27, https://doi.org/10.1017/S0260305500013215

Fretwell P., Pritchard H.D, Vaughan D.G., Bamber J.L., Barrand N.E., Bell R., Bianchi C., Bingham R.G., Blankenship D.D., Casassa G., Catania G., Callens D., Conway H., Cook A.J., Corr H.F.J., Damaske D., Damm V., Ferraccioli F., Forsberg R., Fujita S., Gim Y., Gogineni P., Griggs J.A., Hindmarsh R.C.A., Holmlund P., Holt J.W., Jacobel R.W., Jenkins A., Jokat W., Jordan T., King E.C., Kohler J., Krabill W., Riger-Kusk, M., Langley K.A., Leitchenkov G., Leuschen C., Luyendyk B.P., Matsuoka, K., Mouginot J., Nitsche F.O., Nogi Y., Nost O.A., Popov S.V., Rignot E., Rippin D.M., Rivera A., Roberts J., Ross N., Siegert M.J., Smith A.M., Steinhage D., Studinger M., Sun B., Tinto B.K., Welch B.C., Wilson D., Young D.A., Xiangbin C. & Zirizzotti A. 2013. Bedmap2: improved ice bed, surface and thickness datasets for Antarctica. The Cryosphere 7, 375–393, https://doi.org/10.5194/tc-7-375-2013.

Gudmundsson G. 2003. Transmission of basal variability to a glacier surface. Journal of Geophysical Research—Solid Earth 108, article no. 2253, https://doi.org/10.1029/2002JB002107.

Helm V., Humbert A. & Miller H. 2014. Elevation and elevation change of Greenland and Antarctica derived from CryoSat-2. The Cryosphere 8, 1539–1559, https://doi.org/10.5194/tc-8-1539-2014.

Kamb B. & Echelmeyer K. 1986. Stress-gradient coupling in glacier flow: I. Longitudinal averaging of the influence of ice thickness and surface slope. Journal of Glaciology 32, 267–284, https://doi.org/10.3189/S0022143000015604.

Le Brocq A., Payne A. & Siegert M. 2006. West Antarctic balance calculations: impact of flux-routing algorithm, smoothing algorithm and topography. Computers and Geoscience 32, 1780–1795, https://doi.org/10.1016/j.cageo.2006.05.003.

McInnes B. & Budd W. 1984. A cross-sectional model for West Antarctica. Annals of Glaciology 5, 95–99, https://doi.org/10.1017/S0260305500003566

Mouginot J., Scheuchl B. & Rignot E. 2012. Mapping of ice motion in Antarctica using synthetic-aperture radar data. Remote Sensing 4, 2753-2767, https://doi.org/10.3390/rs4092753.

Nye J. 1952. The mechanics of glacier flow. Journal of Glaciology 2, 82–93, https://doi.org/10.1017/S0022143000033967

Nye J. 1957. The distribution of stress and velocity in glaciers and ice-sheets. Proceedings of the Royal Society of London 239, 113–133, https://doi.org/10.1098/rspa.1957.0026

Payne A. 1999. A thermomechanical model of ice flow in West Antarctica. Climate Dynamics 15, 115–125, https://doi.org/10.1126/science.1208336

Rignot E., Mouginot J. & Scheuchl B. 2011. Ice flow of the Antarctic ice sheet. Science 333, 1427–1430, https://doi.org/10.1126/science.1208336.

Rignot E., Mouginot J. & Scheuchl B. 2017 - present. MEaSUREs InSAR-based Antarctica ice velocity map, version 2.0. NASA Distributed Active Archive Center at the National Snow and Ice Data Center: Boulder, CO. https://doi.org/10.5067/D7GK8F5J8M8R Accessed on the internet at https://nsidc.org/data/nsidc-0484/versions/2 on 20 June 2017.

Roberts J.L., Warner R.C., Young D., Wright A., van Ommen T.D., Blankenship D.D., Siegert M., Young N.W., Tabacco I.E., Forieri A., Passerini A., Zirizzotti A. & Frezzotti M. 2011. Refined broad-scale sub-glacial morphology of Aurora Subglacial Basin, East Antarctica derived by an ice-dynamics-based interpolation scheme. Cryosphere 5, 551–560, https://doi.org/10.5194/tc-5-551-2011.

Young D.A, Roberts J.L., Ritz C., Frezzotti M., Quartini E., Cavitte M.G.P., Tozer C.R., Steinhage D., Urbini S., Corr H.F.J., van Ommen T. & Blankenship D.D. 2017. High resolution boundary conditions of an old ice target near Dome C, Antarctica. The Cryosphere 11, 1897–1911, https://doi.org/10.5194/tc-11-1897-2017.

Zwally H.J., Giovinetto M.B., Beckley M.A. & Saba J.L. 2012. Antarctic and Greenland drainage systems. Accessed on the internet at https://icesat4.gsfc.nasa.gov/cryo_data/ant_grn_drainage_systems.php on 23 February 2019.