A restitution method to reconstruct the 2001–13 surface evolution of Hurd Glacier, Livingston Island, Antarctica, using surface mass balance data

Abstract

The surface restitution method we present reconstructs the evolution of a glacier surface between two time-separated surface topographies using seasonal surface mass balance (SMB) data. A conservative and systematic error analysis is included, based on the availability of surface elevation measurements within the period. The method is applied from 2001 to 2013 at Hurd Glacier (a 4 km² glacier), where we have sufficient SMB and elevation data. We estimate surface elevation changes in two steps: (1) elevation change due to SMB and (2) elevation change due to glacier dynamics. Four different models of the method are compared depending on whether or not accumulation is memorised at each time step and whether they employ balance profiles or SMB maps. Models are validated by comparing a set of surface measurements retrieved in 2007 with the corresponding restituted elevations. Although surface elevation change between 2001 and 2007 was larger than 10 m, more than 80% of the points restituted by the four models showed errors below ±1 m compared to only 33% when predicted by a linear interpolator. As error estimates between models differ by 0.10 m, we recommend the simplest model, which does not memorise accumulation and interpolates SMB by elevation profiles.

List of main symbols

In the observational data

D _aaaa

dataset of surface measurements retrieved during the year aaaa

DSM_aaaa

digital surface model (DSM) of the glacier obtained by geospatial interpolation of D _aaaa. A DSM of a glacier is the digital elevation model (DEM) of its glacier surface elevation

In the restitution method

M_

previous years' accumulations memorised

nM_

previous years' accumulations not memorised

N

number of grid nodes of a DSM

t _i

date of the initial DSM. In the studied case of Hurd Glacier, t _i is 2 January 2001

t _f

date of the final DSM. In the studied case of Hurd Glacier, t _f is 22 January 2013

D _i

dataset of surface elevations at the initial date, t _i. In the studied case of Hurd Glacier, D _i is D ₂₀₀₁

D _f

dataset of surface elevations at the final date, t _f. In the studied case of Hurd Glacier, D _f is D ₂₀₁₃

z _i

elevation of a given point of the glacier surface at t _i

z(t)

elevation of a given point of the glacier surface at a time t

z _f

elevation of a given point of the glacier surface at t _f

Δz(t)

surface elevation change at a given point, from t _i to t

Δz _SMB(t)

component of Δz(t) due to the effect of the cumulative surface mass balance (SMB)

Δz _d(t)

component of Δz(t) due to glacier dynamics

1 Introduction

The use of accurate digital surface models of glaciers, hereafter DSMs, is essential for various glaciological applications, including glacier dynamics modelling, mass-balance studies and climate-forcing analysis. In modelling of glacier dynamics, surface elevation data are used to infer the bed topography from ice thickness data, most often collected using ground-penetrating radar (GPR) (Navarro and others, 2009; Bamber and others, 2013; Fretwell and others, 2013; Lapazaran and others, 2013; Morlighem and others, 2017). When ice thickness data are not available, a DSM of the glacier is essential for the inverse methods that infer the ice thickness distribution from surface topography, velocities, mass balance and glacier thickness changes. For this approach, the error in surface slope is most often the largest source of uncertainty, stressing the importance of an accurate glacier surface topography (Farinotti and others, 2017; Fürst and others, 2017). For ice velocity studies based on synthetic aperture radar (SAR) satellite data, DSMs are used at various data-processing steps. In interferometric methods, the DSM is used for calculating the differential interferogram. In offset-tracking methods, DSMs are used for geocoding and co-registering the imagery. Together with the slant-range and azimuth displacements, DSMs can be used to obtain a 3-D displacement map by combining differential SAR interferometry in the slant-range direction, offset tracking in the azimuth direction and a DSM (Strozzi and others, 2002; Sánchez-Gámez and Navarro, 2017; Gardner and others, 2018).

The comparison of sequential DSMs spaced in time is the most common technique for quantifying glacier volume changes and geodetic mass balance (Lapazaran and others, 2013; Berthier and others, 2014, 2016; Wang and Kääb, 2015; Zemp and others, 2019; Hugonnet and others, 2021). For surface elevation changes and surface mass balance (SMB) studies based on satellite altimetry (e.g. ICESat), it is challenging to compare repeated track profiles due to the relatively large cross-track separation distance between repeating profiles and therefore, DSMs are needed to correct elevation differences caused by the cross-track slope (Moholdt and others, 2010; Kääb and others, 2012; Gardner and others, 2013). The glaciological method for estimating the SMB requires a DSM of the glacier surface for integrating point observations (taken at mass-balance stakes, snow-probing points and snow pits) across the entire glacier surface (Navarro and others, 2013; Zemp and others, 2013).

Contemporaneous DSMs, ice thickness measurements and climate-forcing data can be used together in a single dataset to analyse a glacier's response through time and inform modelling studies of glacier volume evolution under future climate scenarios (Lapazaran and others, 2013; Möller and other, 2016; Morlighem and others, 2017). Distributed SMB modelling (by either temperature index, energy balance or mixed methods) requires a DSM to infer the meteorological conditions over the whole glacier from measurements taken at automatic weather stations (Jonsell and others, 2012). Similarly, DSMs are used for both statistical and dynamical downscaling of regional climate model data to the glacier scale (Mölg and Kaser, 2011; Jarosch and others, 2012; Hannesdóttir and others, 2015).

In most of the above applications, the dates of the available DSMs do not always coincide with the dates of collecting relevant variables for the study under consideration. Hence, it is useful to have a method capable of producing DSMs that are fair representations of the glacier surface topography at any given time. For instance, Frappé and Clarke (2007) iteratively generated a time-evolving DSM (and DEM of ice thickness) of Trapridge Glacier (Yukon, Canada) based on Bayesian kriging by merging photogrammetric DSMs and annual elevation data at different sets of surface points. This was achieved using three data types: ground survey profile lines, ground survey flow pole positions and aerial photography. For every time step, they employed the DSM from the previous iteration as a prior model, which was updated as a new model with the current field data (ground survey profile lines and ground survey flow pole positions). This new model became the prior model in the next time step. However, this method can only generate the DSM at the specific date of the field data measurement.

Hence, we present a novel method, hereinafter restitution method, to reconstruct the evolution of a glacier surface between two time-separated DSMs, using seasonal SMB data. The method is referred to as ‘restitution’ because it aims to restore surface elevations (and ice thicknesses) to their value at a previous date during the period. At any point on the glacier surface, the restitution method calculates the local elevation change between two dates in two steps: (1) estimation of the elevation change due to SMB and (2) estimation of the elevation change due to glacier dynamics. The method requires the availability of seasonal SMB data during the period between both DSMs. These SMB data can be obtained by methods such as the glaciological method (Østrem and Brugman, 1991; Hagen and Reeh, 2004; Cogley and others, 2011), comparison of radar and probing data applied to internal radar reflections (e.g. Callens and others, 2016), inversion techniques applied to internal radar reflections (e.g. Köhler and others, 1997) or inversion techniques applied to several DSMs of a glacier and the digital terrain model of its bedrock (e.g. Välisuo and others, 2017). We apply the restitution method to Hurd Glacier for the period 2001–13, for which we have sufficient available data to test and validate various implementations of our method.

2 Geographical setting

Hurd Glacier is one of the main glacier basins of Hurd Peninsula Ice Cap (62°39–42′ S, 60°19–25′ W) located on Livingston Island, South Shetland Islands, Antarctica (Fig. 1). It is a land-terminating glacier with a tapered terminus at Sally Rocks, Argentina, and Las Palmas lobes, and its dominant flow direction is to the southwest. Hurd Glacier covers ~4 km² and spans an elevation range from sea level to ~340 m a.s.l. A local ice divide with elevations of 230–340 m a.s.l. separates Hurd and Johnsons glaciers. The main lobe of Hurd Glacier is ~3 km long and ~2 km wide at the divide with Johnsons Glacier, decreasing to 0.6 km at the terminus. The direction of the prevailing wind in Hurd Glacier is northeast, with an average speed of ~4 m s⁻¹. In the first decade of the 21st century, with average precipitation days per year of 270–290 measured at Faraday/Vernadsky station, the annual temperature in the glacier was −1.1°C, with average summer and winter temperatures of 2.8 and −4°C, respectively, and relative humidity of 80% (Navarro and others, 2013). The maximum ice thicknesses of Hurd Glacier, ~200 m, are found in the accumulation area (e.g. Molina and others, 2007). The Spanish Antarctic station Juan Carlos I (JCI, Fig. 1) is close to Hurd Glacier, which has been the subject of numerous studies during the last two decades. Our research team from Universidad Politécnica de Madrid has undertaken topographic measurements and SMB studies on Hurd Glacier since 2000 (Navarro and others, 2013) and GPR studies since 1999 (Navarro and others, 2009).

Fig. 1. (a) Location of Livingston Island in the South Shetland Islands. (b) Location of Hurd Peninsula on Livingston Island (orthophoto from Sentinel-2 15 November 2016 image). (c) Location of Hurd and Johnsons glaciers in Hurd Peninsula, with the position of Juan Carlos I Antarctic Station (JCI, yellow dot), the mass-balance stakes in December 2015 (red dots) and the surface elevation map (based on a survey during summer 1998/99 and 2000/01). The dashed blue line indicates the ice divide separating Hurd and Johnsons glaciers.

3 Observational data

The primary input data for this study are: surface measurements from field campaigns 2000/01, 2007/08 and 2012/13 (D ₂₀₀₁, D ₂₀₀₇ and D ₂₀₁₃); a sequence of SMB and position measurements retrieved since 2001 at a network of stakes on the glacier (30 stakes in Hurd, 59 in total, as shown by red dots in Fig. 1c, spanning an elevation range of ~30–320 m a.s.l.); and snow and firn density pits and cores retrieved in Hurd Glacier by the authors since 2000 and by Furdada and others (1999) in the divide between Hurd and Johnsons glaciers.

Our surface restitution method requires two DSMs for Hurd Glacier, DSM₂₀₀₀ and DSM₂₀₁₃, defining the limits of the period under study. Accordingly, we generate them with spatial interpolation using ordinary kriging (OK) (see Section 4.5) of our datasets D ₂₀₀₁ and D ₂₀₁₃ (Figs 2a, c). However, two more sets of surface measurements are employed: one (D ₂₀₀₇, Fig. 2b) covers the entire glacier at an intermediate date (2007), and the other a time series of surface elevation measurements at the stakes (red dots in Fig. 1c). The first set is used for validating the model, while the latter is used for error estimation. Neither is applied by the restitution method.

Fig. 2. Black dots represent the datasets of surface elevation GNSS measurements on Hurd Glacier: (a) 308 points in D ₂₀₀₁; (b) 885 points in D ₂₀₀₇ and (c) 7174 points in D ₂₀₁₃. The DSMs at the ends of the period are shown in background: (a) DSM₂₀₀₀ and (c) DSM₂₀₁₃. The axes represent UTM coordinates, zone 20 South.

3.1 Glacier surface measurements

Figure 2 shows the layout of three datasets of surface elevation measurements used in this paper. D ₂₀₀₁, D ₂₀₀₇ and D ₂₀₁₃ were retrieved ~2 January 2001, 18 December 2007 and between 22 January and 16 February 2013, respectively. The latter includes corrections due to the ablation during the period of measurements, thus representing the surface on 22 January 2013. Such corrections were estimated by daily measurements of the surface elevation at a set of control points. The datasets are available in Rodríguez and Navarro (2017) (https://doi.org/10.1594/PANGAEA.873067), although D ₂₀₀₁, D ₂₀₀₇ and D ₂₀₁₃ are the corresponding subsets of points within Hurd Glacier (Fig. 2). The surface elevation data in D ₂₀₀₁, D ₂₀₀₇ and D ₂₀₁₃ are based on differential Global Navigation Satellite System (GNSS) surveys over snow and ice. Different models of Trimble GNSS have been used over the period with all measurements being retrieved in real-time kinematic or processed using the GNSS base station of the Juan Carlos I Antarctic Station. D ₂₀₀₁ was surveyed by foot with horizontal and vertical accuracies of ±7 cm SD_h and ±10 cm SD_v. The 2007 survey was performed by snowmobile with accuracies of ±10 cm SD_h and ±14 cm SD_v. D ₂₀₁₃ was surveyed by foot and snowmobile including photogrammetric restitution at the front of Johnsons Glacier, acquiring final horizontal and vertical accuracies of ±65 cm SD_h,v (Rodríguez, 2014).

Due to the ice roughness and snow compressibility being greater than the measurement accuracy, and based on our experience, we conservatively estimate measurement errors of ±25 cm SD_h,v in both D ₂₀₀₁ and D ₂₀₀₇. However, since D ₂₀₁₃ also includes ablation corrections, we estimate that its measurement error is ±65 cm SD_h,v (thus assuming as standard the maximum error of D ₂₀₁₃).

3.2 Surface mass balance data

We denote mass balance as b to refer to its value at a point on the glacier's surface, while B refers to the value integrated over the entire glacier surface (Cogley and others, 2011). The mass-balance year is assumed to be the hydrological year defined by the World Glacier Monitoring Service (http://www.wgms.ch/downloads/WGMS_GuidelinesforDataSubmission.pdf) for the Southern Hemisphere (1 April to 30 March). Furthermore, we assume that the extended winter season spans from 1 April to 30 November, and the summer season, the remaining months. Thus, the mass-balance year 2013, which starts 1 April 2012 and ends 31 March 2013 is divided into a winter season 2013, from 1 April to 30 November 2012, and a summer season 2013, from 1 December 2012 to 31 March 2013 (e.g. Huss and Hock, 2015). We use the subscripts “w” and “s” to represent winter and summer values. We also add a numeric subscript to represent the corresponding mass-balance year (i.e. B ₂₀₀₇, b ₂₀₀₇, b _w2007 and b _s2007 represent the variables for the mass-balance year 2007). Winter (b _w) and summer (b _s) SMB points were obtained using the direct glaciological method (Cogley and others, 2011).

The SMB dataset is based on a time series of measurements performed concurrently with those in Section 3.3 in a net of stakes (red dots in Fig. 1c) which were installed on Hurd Glacier at the end of January 2001 (middle of summer 2001). The target of these SMB measurements fits with the requirements of the restitution method, i.e. SMB values associated with each summer season's start and end. However, due to logistical constraints, our fieldwork often starts after 1st December and concludes before the actual end of the melting season. Hence, we must estimate a correction of the SMB values to account for mass changes before and after our fieldwork campaign measurements. The post-field season melting is quantified as the discrepancy between the snow depth of the subsequent year and the difference in stake height from the last measurement of the season to the first of the subsequent year (Navarro and others, 2013). The SMB at the beginning of the summer season is corrected to 1st December by linear extrapolation of the summer measurements. Therefore, the corrections affect both b _w and b _s for any given year. A more detailed description of the mass-balance measurements, the net of stakes and its maintenance can be found in Navarro and others (2013). Due to the unavailability of b _s2001 measurements, we have estimated its value as the average of the corresponding values between 2002 and 2009, which was a period with relatively stable SMB data (those following 2009 showed an increasing trend of SMB).

Hurd Glacier outlines have not changed significantly during the period 2000–13. Rodríguez-Cielos and others (2016) analysed the front position changes over the period 2000–12, obtaining maximum retreats in the central flowlines of Argentina, Las Palmas and Sally Rocks tongues of 35, 41 and 83 m, respectively, implying changes in total glacier area of ~1%. Almost all these changes took place during 2000–09, and the fronts remained stable during 2009–13. The near-stationary position of these glacier fronts started in 2009 and lasted until 2016 (unpublished data from Francisco Navarro), coincident with the period of most intense regional cooling and positive SMBs (Oliva and others, 2017). Despite area change being negligible, the glacier-wide SMBs presented here have been computed using updated glacier front positions for each year, thus accounting for any changes in SMB due to glacier area variations. In the period of our analysis, the average glacier-wide SMB has been very close to equilibrium, at −0.13 m w.e. a⁻¹, although with a rather large variability, which are larger for the summer balances (Table 3 of Appendix A). The average equilibrium line altitude (ELA) was 203 m a.s.l. and the average accumulation–area ratio (AAR) was 49%.

Given the somewhat similar characteristics of the various outlets of Hurd Glacier (all outlets are land terminating and flow approximately westwards) and the uniform coverage of the glacier by the net of stakes (Fig. 1), we can assume that the balance profile adequately captures the spatial variability, perhaps except for the steeper surface slopes of Argentina and Las Palmas outlets.

3.3 Surface data points at intermediate dates

We use a time series of surface elevation measurements from 2001 to 2013, at ~30 points (those with stake; red dots in Fig. 1c) on Hurd Glacier. Whenever a mass-balance measurement is performed during the summer season under study, the surface elevation at each stake is also measured using differential GNSS techniques. Thus, a surface elevation measurement is performed at each stake two to four times a year. A detailed description of the mass-balance measurements, the net of stakes and its maintenance can be found in Navarro and others (2013). Since this dataset is not employed to tune the model (Section 4.1.1), it is used as a comparison dataset to measure the restitution error (Sections 4.6, and C1.3 and C2 of Appendix C).

Although these data points are derived from differential GNSS measurements with a precision of ±4 cm SD_h,v, measurements are taken at the top of the stakes, while the surface is estimated by considering the dip angle, the azimuth and the distance from the surface to the top of the stake. Consequently, due to the inclination of some stakes, the roughness of the ice when it is exposed, and the daily variation of the snow cover, we prefer to be conservative, thus considering the measurement error of these surface data to be ±0.5 m SD_v.

4 Methods

We propose an elevation restitution method applied to any (x, y) point on a glacier surface. It estimates the surface elevation of the point (x, y) at any time t within a period between two existing surface topography datasets, D _i and D _f (corresponding to times t _i and t _f), using the seasonal SMB data within this period. It is possible to generate any DSM of the glacier using any grid within the domain at any time within the period, just by applying the restitution method to every point of the grid at the selected time. Applying this methodology to the ice thickness measured at a certain point and time can restitute the ice thickness value on that point of the glacier at a different moment.

Let (x, y) be a point on the glacier surface with elevations z _i and z _f at times t _i and t _f. Then, if (x, y) does not belong to D _i (or to D _f), its elevation z _i (or z _f) is interpolated. Once z _i and z _f are known, the evolution of the elevation change between them is computed to estimate the elevation of (x, y) at time t within the period between t _i and t _f.

4.1 Components of the surface elevation changes

We split the surface elevation change at point (x, y) between t _i and t in two components: the elevation change due to the effect of the cumulative SMB, denoted as Δz _SMB(t); and the elevation change due to glacier dynamics (emergence or submergence), denoted as Δz _d(t)

(1)$$\Delta z( t ) = z( t ) -z_{\rm i} = \Delta z_{{\rm SMB}}( t ) + \Delta z_{\rm d}( t ) .$$

To estimate Δz _SMB(t) at point (x, y), its initial elevation, z _i, is increased and decreased every winter and summer using the accumulation and ablation seasonal data until the time, t ∈ [t _i, t _f] is reached (Fig. 3a). Seasonal SMBs are distributed along their corresponding season by time interpolation (from zero at the beginning of the season to the observed value at the end). Here, we use cubic spline interpolation in time. However, the method has little sensitivity to the choice of the time interpolant for SMB, thus obtaining minimal difference in the results when using piecewise-linear time interpolation. The conversion of SMB data to elevation changes (m w.e. to m) is carried out using the density of the material (Section 4.3).

Fig. 3. Process for estimating the components of the surface elevation change. The dates shown are those of our case study. (a) Elevation change due to SMB. The mass balance (m w.e.) is converted into an elevation change (m) using the density of the material. (b) Surface elevation change between the initial and final stages. (c) The difference between both elevation changes is the dynamic component of the elevation change. To simplify the figure, glacier dynamics is represented with a linear evolution in time. The nature of the glacier dynamics evolution is explained in Section 4.1.1, including the study of the case of the dynamics in Hurd Glacier.

The estimate of Δz _d(t) is performed in two steps. The first step is to estimate Δz _d(t _f) by applying Eqn (1) to t _f as the difference between z _f and the result of summing z _i to the elevation changes due to SMB until t _f is reached (Fig. 3):

(2)$$\Delta z_{\rm d}( {t_{\rm f}} ) = z_{\rm f}-[ {z_{\rm i} + \Delta z_{{\rm SMB}}( {t_{\rm f}} ) } ] .$$

Then, Δz _d(t _f) can be interpolated in time at any date, t, within the period. In the absence of a better approach, we assume a constant dynamics rate between both dates (Fig. 3c), thus obtaining Δz _d(t) by linear interpolation, as

(3)$$\Delta z_{\rm d}( t ) = [ {( {t-t_{\rm i}} ) /( {t_{\rm f}-t_{\rm i}} ) } ] {\rm \;}\Delta z_{\rm d}( {t_{\rm f}} ) .$$

In case of having more information on the dynamics evolution of the glacier (e.g. percentages of movement between summer and winter seasons; intermediate measurements of the surface elevation), a more precise interpolation in time for the dynamic component could be implemented (see Section 4.1.1).

Although Δz _SMB has been calculated up to t _f for estimating Δz _d(t), the method also requires its value at t, so as to obtain Δz _SMB(t). Finally, adding both Δz _SMB(t) and Δz _d(t) to the initial elevation of the point at t _i, we estimate its elevation at date t

(4)$$z( t ) = z_{\rm i} + \Delta z( t ) = z_{\rm i} + \Delta z_{{\rm SMB}}( t ) + \Delta z_{\rm d}( t ) .$$

When using balance profiles by elevation ranges, a problem occurs during the process of estimating Δz _SMB(t) up to t _f. The surface elevation of the point is required by the method at the end of each season to obtain the corresponding seasonal SMB value. However, it will not be possible to know the surface elevation of the point until the dynamic component is determined at the end of the process. To avoid such a recursive problem, since the sensitivity of the SMB to small errors in elevation is negligible, we use a simple interpolation between z _i and z _f in proportion to the date at the end of the corresponding season between t _i and t _f as the input elevations for the profiles of the seasonal SMB. When a large sensitivity of the SMB to small errors in elevation is expected, an iterative process should be implemented in order to obtain a more accurate sequence of seasonal SMB, using as the input elevations for the balance profiles during the next iteration, the elevations already estimated in the previous one. This problem is however absent when using maps of SMB.

4.1.1 Tuning the dynamic component interpolator in time

The simplest time interpolator for the dynamic component of the elevation change is a linear relationship between t _i and t _f, as stated in Eqn (3). To generate a more precise interpolator, we tune its shape to obtain the best-fit estimate for a subset of the data on the glacier surface (Section 3.3). In the shape-tuning process, we only employ surface elevation measurements (Section 3.3) at stakes whose sequence of measurements reconstructs the entire period under study (2001–13) and study the sequence of each stake separately. First, we estimate the dynamic component of the elevation change for each stake measurement, avoiding any interpolation in time. We achieve this using the evolution process explained above, up to the step represented in Eqn (2), but replacing t _f with the date of this measurement. Then, we repeat the above procedure for each datum along the sequence of measurements at this stake to determine the evolution of the dynamic component. To allow the shape comparison between stakes, we normalise the results by the value of the dynamic component in 2013 (usually, the maximum in the sequence). We separately repeat this procedure for all the selected stakes and compare their normalised evolutions. The best shape of the interpolator is given by the best fit (by least-squares fitting) of the normalised evolutions of the dynamic component.

Note that we follow the dynamic movement of the stakes (Lagrangian perspective) in this procedure instead of treating them as fixed points (Eulerian perspective), although we need the latter. In stakes with slight movement, both perspectives give similar results. On the contrary, in stakes with significant movement, their downward displacement could accelerate the dynamic component, thus introducing non-realistic increments of parabolic shapes.

We applied this method using the measurements on Hurd Glacier at intermediate dates (described in Section 3.3). Due to the slight movement of the stakes on Hurd Glacier, there is no significant change in the results owing to the Lagrangian or Eulerian perspectives. We found that the linear interpolator is the best fit interpolator in time for the dynamic component of the elevation change (Fig. 4). Therefore, we use the linear interpolator in time for the dynamic component of the entire glacier between 2001 and 2013. This dataset is employed for the error estimate since none of the intermediate surface data are used to tune the interpolator in time.

Fig. 4. Black dots represent the dynamic component of the surface elevation changes (following the simplest model of the restitution method, nM_BP) at the set of stakes with surface elevation measurement along the period, normalised by the value for the stake in 2013. The red line represents the linear regression estimate. The beginning of winter and summer seasons are shown with blue- and orange-dashed lines, respectively. Dates in the X-axis indicate 1st of December for the corresponding year.

4.2 Density–age law in Hurd Glacier

A density–age law function (DAL) is applied to estimate the density of the material based on the age of its accumulation. The densities determined from the DAL are used to convert SMB to elevation change throughout the time period. The DAL starts with a density of 500 kg m⁻³ for the current year's snow (age is 0) and 520 kg m⁻³ for the first year's firn (mean values obtained from our historical measurements in snow pits at the accumulation zone, excavated to the depth of the last summer layer). Then it reaches 830 kg m⁻³, the density at which firn becomes glacier ice (where voids cease to form a connected network; Cogley and others, 2011) at the same age as in Furdada and others (1999), i.e. in ~24 years, but rises exponentially until year 24 (data from M. Pourchet and J.M. Casas reported in Ximenis, 2001). Thus, the DAL describes the density of the material that is t years old (accumulated between t and t + 1 years ago) with t = 0, initial year's snow: ρ = 500 kg m⁻³ and $\forall t \ge 1\colon \rho = 0.4707\;t^{0.1852}\;{\rm kg}\;{\rm m}^{{-}3}$ (limited to 900 kg m⁻³, which is reached after 30 years).

4.3 Conversion from SMB to elevation change

The restitution method involves a sequence of seasonal SMBs distributed along their corresponding season by time interpolation. Additionally, at the point under study, the cumulative SMB is transformed to elevation change at any time t, after adjusting the density of the accumulated or ablated material.

We compare the results of two different density models to study whether it is worth using a simple model or if it is necessary to use a more precise but more complex one. In both instances, the count of seasonal SMBs starts at the beginning of the winter season of the first year, i.e. the winter's beginning prior to t _i, considering the snow previously accumulated during this first winter.

Both density models differ on whether they use accumulation memory or not. The models with memory implement a memory array, accounting for the ablation of previous accumulations (using their respective densities). However, there should be no differences between models with or without memory, neither in points with permanent accumulation nor in points with permanent ablation. In the former, the positive net accumulation of every year guarantees no ablation of past remains, while in the latter, there is no remaining accumulation from any previous year, thus constantly melting ice. Only at intermediate zones, where net accumulation or ablation depends on the year, memory could provide differences in the results.

4.3.1 Density model nM: no memory – two densities (ice and snow)

The simplest density model only considers two materials: snow (~500 kg m⁻³; see Section 4.2) for the current year accumulation and ice (~900 kg m⁻³) to characterise the whole accumulation from previous years. This model only memorises the accumulation within the current year, which is reset at the end of each summer, thus not memorising any previous snow.

The algorithm assumes snow density for the winter's accumulation, ablating some of this snow during summer. Only in cases of net ablation, this year-round snow column melts completely and continues melting part of the previous material, which is underneath and characterised by ice density.

4.3.2 Density model M: accumulation memory (density–age law)

For any point on the glacier where the restitution method is applied, this model memorises the sequence of their accumulated masses and ages, using the ice density as an initial condition. This memory evolves as accumulation and ablation progress over time, as detailed below.

The algorithm starts similar to the previous nM model: winter's accumulation is estimated with snow density, and only in cases of net ablation, the year-round snow column melts completely and continues melting part of the previous material, which is underneath. However, this underlying material is not characterised here with ice density. Following the DAL, the M model characterises the glacier material (snow, firn or ice) with the density corresponding to its age. This mechanism is implemented using a memory array that stores the amounts of remaining material from previous years at the point under study (see Appendix B for technical details about implementing the accumulation memory and the use of the DAL).

4.4 Feasible SMB data

In order to reconstruct the evolution of the SMB at any point on the glacier surface, there must be enough SMB data to estimate the SMB at any point and time during the period [t _i, t _f]. We use the SMB data in two different formats to compare the differences in model results. On the one hand, we used maps of b _w and b _s for each season obtained by direct interpolation (Section 4.5) of the measured point balance values at the stakes (Fig. 1c). On the other hand, we used balance profiles of b _w and b _s for each season and elevation range (20 m width), generated by averaging the former maps by elevation ranges. Note that the latter is the most common format to represent the SMB of a glacier.

4.4.1 SMB data type MAP: map of SMB

To generate SMB maps, we employ point SMB measurements at a net of stakes in the glacier (red points in Fig. 1c), interpolating (Section 4.5) a map of SMB for each season. Therefore, the SMB map is point dependent rather than elevation dependent like a balance profile.

4.4.2 SMB data type balance profile (BP): interpolation of SMB profiles by elevation ranges

The elevation range of the ice mass is split into a set of ranges 20 m wide, and the values in the balanced profile represent the SMB in the corresponding range. We use a set of these balance profiles for the entire period [t _i, t _f], one for each season. At each balance profile, a simple piecewise-linear interpolation of the SMBs assigned to the centre of the bands is then applied to avoid discontinuities.

4.5 Geospatial interpolation

Most data utilised are poorly distributed, with many in-line points and significant separation between lines (e.g. see Fig. 2). Thus, we apply geospatial interpolation to estimate the glacier surfaces corresponding to the dates t _i and t _f. In selecting the best covariance function, we evaluate the degree of fit of four covariance functions (spherical, stable, Gaussian and exponential) with the experimental semivariogram. We observe that the stable function gives the best fit to our data, especially at low lag distances, which are the most relevant in geospatial interpolation. In selecting the best interpolation method, we considered the OK and universal kriging (UK) methods (e.g. Wackernagel, 2003). Although OK is a robust kriging method and is commonly used, it operates on the assumption of a stationary mean value. When applied to surface elevations on a glacier with surface slope, this assumption of spatial stationarity could be violated. So UK, which accounts for a polynomial trend in the data values, could be a good choice. In our case, our test showed that due to the very gentle surface slope of Hurd Glacier, and due to the different trends depending on the lobe, there is no significant difference between UK and OK. To evaluate it, we generated a synthetic DSM of Hurd Glacier from field measurements, and took samples from the synthetic DSM to generate a synthetic dataset. Both kriging methods were applied with the stable covariance model on the synthetic dataset at the grid nodes of the DSM. Then, comparing the results of both kriging methods to the synthetic DSM, we observed that there was no significant differences, although OK slightly outperformed UK. Hence, we used OK with a stable covariance model on a unique grid of points to interpolate our datasets D ₂₀₀₁ and D ₂₀₁₃, thus respectively obtaining DSM₂₀₀₁ and DSM₂₀₁₃. This type of geospatial interpolation has also been applied to generate the SMB maps (Hock and Jensen, 1999), using the point SMB measurements at the net of stakes (Section 4.4.1).

Following a modified version of the technique described in Lapazaran and others (2016), we systematically estimate both the interpolation error and the bias introduced by each geospatial interpolator. All the interpolated values in this study are bias-corrected using this technique. The technique in Lapazaran and others (2016) is based on performing a successive blanking process (temporary elimination of data) in the dataset, using different blanking radii around each data point. Thus, both the interpolation bias and the interpolation error are estimated as functions of the blanking radius. This paper uses a variant of this technique involving calculating bias and interpolation error as functions of the distance from the nearest data rather than the blanking radius.

4.6 Error analysis

The error analysis of this study provides estimates of all error components at any time during the period [t _i, t _f]. The details of the technical implementation of the error analysis can be found in Section C1 of Appendix C.

In summary, we split the restitution error in two: the error (at any time) due to the errors at the ends of the period; and the error in prediction. The first is linearly propagated from one end to the other, whereas the characterisation of the second is parabolic with a vertical axis, null at both ends of the period and maximum at the centre. Thus, the theoretical shape of the restitution error is known, since the combination of both independent components results from their squared quadratic summation. The square of the restitution error results in a fourth-degree polynomial in time.

We obtained the least-squares fit of the set of squared restitution errors along the period measured, calculated using the glacier surface elevation measurements described in Section 3.3. From this best fit, we obtain the estimate of the parabolic error in prediction. Then we re-estimate the errors at the ends of the period, and split them in two independent errors: the measurement error and the interpolation error at the initial and final surfaces. To estimate the latter, we follow a modified version of the technique described by Lapazaran and others (2016) (see Section 4.5).

5 Results and discussion

We compare the results of four restitution models, with or without memory (M or nM), and with BPs or with maps (MAP): nM_BP, M_BP, nM_MAP and M_MAP. They are provided by the combination of two models of conversion from SMB to elevation change (Section 4.3: nM and M), each using both types of SMB data (Section 4.4: BP and MAP). It is worth noting that the most simplistic model among all four models is nM_BP, with M_MAP being the most sophisticated. We also compare the models with linear interpolation in time (LIN) between the two ends of the period (dates in which the glacier surface elevation was measured).

This section applies the proposed models of the restitution method to Hurd Glacier between 2 January 2001 (t _i) and 22 January 2013 (t _f), thus D ₂₀₀₁ being D _i and D ₂₀₁₃ being D _f. First (Section 5.1), we show the surface evolution during the period in a set of eight points on Hurd Glacier as a result of the restitution method, revealing the different surface oscillations depending on whether the point is in accumulation or ablation zone. Then (Section 5.2), we validate the method by comparing the results of the four restitution models with surface measurements of 2007. We also include an error study (Section 5.3, with details in Section C2 of Appendix C) for the four models of the restitution method, and contrast such error estimates by comparing them with the discrepancies between the surface measurements of 2007 and the restituted values (Section 5.4).

5.1 DSMs of Hurd Glacier at any date

To illustrate the surface evolution of Hurd Glacier throughout the period, we selected eight points spaced ~400 m apart along the main flowline of the glacier, spanning its elevation range. Figure 5 presents the evolution of the elevation change of these eight points, obtained by restitution of their elevation between the initial and final dates using the simplest restitution model, nM_BP. These evolutions are shown by twice-a-month restitutions, connected with small lines. The resulting zigzag lines represent the glacier surface evolution with seasonal oscillation.

Fig. 5. Model nM_BP results depicting the evolution of surface elevation along a centre profile at eight points spaced ~400 m apart. Each of the eight lines is composed of a sequence of restitutions every ~15 d, connected with small lines. The resulting zigzag lines represent the seasonal evolution of the glacier surface elevation at the eight points. Dates in the X-axis indicate 1st of December for the corresponding year.

5.2 Restitution vs DSM₂₀₀₇

We use the dataset D ₂₀₀₇ to evaluate the differences between the restituted values and the surface elevations measured in 2007 by comparing the four restitution models (nM_BP, nM_MAP, M_BP and M_MAP). At the same time, we compare these models with the result of a simple linear interpolation in time between the initial and final elevations (LIN).

Figure 6 illustrates the set of comparisons over the Hurd Glacier map. All four models are capable of correctly reconstructing glacier surface elevation estimates in most areas, while the linear interpolation is worse, with a clear positive bias. However, all models show their most significant errors in the lateral lobes of Hurd Glacier. The differential characteristics of these areas are their large gradients in surface elevation, their thin ice and the proximity to the boundary with surrounding rocks. Thus, both the SMB and the dynamics of the glacier in these areas are expected to vary.

Fig. 6. Errors resulting from the four restitution models, and a linear interpolation in time, when compared to the 2007 surface elevation dataset. M or nM: model with or without previous years' accumulations memory. BP or MAP: model based on the use of balance profiles or maps of SMB. LIN: linear interpolation in time between initial and final elevations of the point (used for comparison, not as a model). The four models give similar errors: nM_BP, −0.28 ± 0.96 m; nM_MAP, −0.42 ± 1.03 m; M_BP, −0.48 ± 0.94 m; M_MAP, −0.45 ± 1.01 m. LIN shows larger errors and a clear positive bias (reddish), 1.27 ± 0.98 m. To make more evident their differences, errors are shown in blue or in red when larger than a threshold of 1 m (negative or positive, respectively). The axes represent UTM coordinates, zone 20 South.

The results in Figure 6 reveal that model nM_BP shows the most significant density of approximately zero-error (yellowish) points among the rest, with 70% of the data falling in the ±0.5 m error range. Consequently, comparing BP models (nM_BP vs M_BP), nM results show errors ~8% lower than M. However, we realise that M_BP errors in the lobes of the glacier are ~2% lower in comparison with nM_BP. Making the same comparison between MAP models (nM_MAP vs M_MAP), we found no significant differences between their results. There is no significant error in the RMS value when nM models are compared to the M models. However, when MAP models are compared to BP models, these latter show RMS error of ~0.1 m less than that of MAP models, although, along the lateral margins, this error increases to ~0.3 m. Hence, MAP models have the largest errors when comparing those using the same density model, nM_BP vs nM_MAP and M_BP vs M_MAP.

To establish an overall comparison of the errors found, we present statistics in Table 1 and two graphs in Figure 7. Looking at Table 1 and Figure 7a, we note a similar range of errors between models when comparing minimum or maximum values. However, the nM_BP model has the lowest bias (absolute value), and BP models generally show a lower dispersion as compared to MAP models (see also Fig. 7a). Observing the RMS value which considers both the bias and the std dev., nM_BP shows the best overall results, while MAP models show the worst ones. The boxplot representation of the errors (Fig. 7a) reveals that nM_BP model shows the smallest bias and dispersion, thus giving the best overall results when validated against 2007 surface data. MAP models capture the zero-error line between quartiles Q1 and Q3, although it is due to their larger dispersion, and the four models capture it between whiskers (whiskers fall at the farthest values from Q2, between Q1–1.5(Q3–Q1), and Q3 + 1.5(Q3–Q1)). On the contrary, linear interpolation does not capture the zero-error line due to its large positive bias.

Table 1. Parameters of the errors (in m) when comparing the restituted values using the four restitution models, with the surface elevations measured in D ₂₀₀₇

For each model, minimum error, maximum error, std dev. of the error, average error (restitution bias) and RMS error are given.

Fig. 7. (a) Boxplot representation of the errors for the four restitution models and the linear interpolation in time. (b) Percentage of predicted points with lower absolute error than a given threshold for a continuous sequence of thresholds and for the four restitution models plus the linear interpolation in time. Taking a vertical line at any threshold, it gives the percentage of the points predicted by each model in which the absolute error is lower than this threshold.

We have counted the percentage of points with absolute error falling under a sequence of thresholds, showing the results in Figure 7b as a continuous graph for the represented range of error threshold. For instance, taking a vertical line at 1 m in Figure 7b, it gives that only 33% of the points predicted by LIN have error lower than ±1 m, while more than 80% of the points predicted by each of the four models have error lower than ±1 m. Figure 7b reveals that the best restitution model when compared to 2007 measurements is nM_BP, although M_MAP has a larger number of points with almost no error. Overall, BP models show less error than the MAP models. LIN produces the greatest errors (Fig. 7).

5.3 Related errors, obtained from sets of discrepancies

The detailed results of the error estimate for the four models when applied to the case study are shown in Section C2 of Appendix C. The four models give similar restitution error estimates with maximum differences between models ≲0.1 m, somewhat before the middle of the period, in 2006 (being an error ~1 m larger in the case of LIN). The estimates of the restitution error of the four models start in all cases at 1.28 m in 2001, increase to ~1.4 m in 2006 and diminish to 0.75 m in 2013.

It is valuable to compare the error estimates for the four models with the discrepancies found between the 2007 elevation measurements and the corresponding restituted elevations.

5.4 Estimated errors vs discrepancies with 2007 data

Table 2 shows the errors measured by comparing restitution vs 2007 surface data (Section 5.2; Fig. 7), and those estimated by their characterisation in Section C2 of Appendix C. On the one hand, we have considered that the std dev. of the uncertainty in any random variable is used to characterise its error. On the other hand, we know that ~68.3% of a standard normal distribution falls closer than 1 std dev.. Thus, a good way to characterise the measured error in Section 5.2 is to find the smallest error threshold which guarantees that at least 68.3% of the points have absolute restitution error lower than this threshold.

Table 2. For each model: the smallest threshold (with a precision of 1 cm) which guarantees that at least a 68.3% of the points are restituted with absolute error lower than such a threshold, and the restitution errors (ɛ_R) estimated in Section C2 of Appendix C

We found out that the errors measured by comparing restitution vs 2007 surface data are much lower (although of the same order) than those we are estimating, showing an overestimation tendency on the estimated errors for the four models. Thus, we can be confident that the method gives conservative error estimates.

5.5 Application of the restitution method to other glaciers

The restitution method presented here can be applied to any glacier with stable dynamics (i.e. not necessarily close to equilibrium, but avoiding surges during the period under study) where the following measurements are available: a surface elevation model at each of the ends of the period, SMB data (balance profiles or maps) within the dates of the two surface elevation models and a series of surface measurements within the period (the more frequent and uniform the data coverage is over the entire glacier, the better) to be able to characterise the shape of the dynamics. In order to apply any model with accumulation memory (M models), it is necessary to estimate the DAL for the glacier. Additional surface elevation measurements throughout the period are necessary to carry out the error study. SMB measurements should at least cover the entire elevation range of the glacier under study. To identify the seasonal wave, SMB will need to be measured twice per year (one at the end of each balance season). However, if surface elevation changes are only intended to be monitored annually, measuring SMB once per year could be enough.

6 Conclusions

We demonstrate the ability to reconstruct the continuous evolution of a glacier surface topography between two timestamps defined by an initial and a final DSM, using a novel restitution method on Hurd Glacier between 2000 and 2013, which requires seasonal SMB data during the period. Our restitution method enables the acquisition of a DSM (as well as a DEM of ice thickness) using any grid within the glacier at any date between the two timestamps. We compare four different models of the restitution method, which differ on how glacier surface elevations are interpolated from SMB data and whether or not accumulation is memorised at each time step. As the error between all models differs by 0.10 m, we recommend the simplest model, nM_BP, that relies on the balance profiles with no memory of the previous years' accumulations.

The restitution method implements the estimate of the changes on the glacier surface elevation between the initial and final DSMs in two steps: (1) estimate the elevation change due to SMB and (2) estimate the elevation change due to glacier dynamics. The method employs a sequence of surface elevation measurements to infer the shape of the evolution of the glacier dynamics. In the case of Hurd Glacier, the results are consistent with linear modelling of glacier dynamics evolution. Furthermore, it allows extrapolating the restitution method from the period between DSMs if information on the corresponding seasonal SMB is available. Although extrapolations could increase the error, it would not be the case for small extensions since the prediction error is null at the ends of the period (dates with known surface topography). Consequently, this method also generates valuable restitutions at dates close to that period. However, our error model does not estimate the error in such a case.

We have validated the four models of the restitution method, measuring their errors by comparing their surface elevation restitutions to 2007 measurements (Section 5.2). The four models are effective at predicting the surface at any date between surface constraints with satisfactory errors. Error estimates for our preferred model are <0.8 m ~80% of the time. Errors for the four restitution models are statistically indistinguishable, but the data suggest that the most sophisticated model produces the smallest errors. All restitution models significantly outperform linear interpolation (Table 1).

Our error analysis provides estimates of each component (Section 4.6, detailed in Section C1 of Appendix C). When applied to Hurd Glacier, the error analysis estimates the restitution errors of the four models based on the comparison with surface elevation measurements throughout the period. All the four models present similar error estimates, with a difference <0.1 m (Fig. 9 in Section C2 of Appendix C). Therefore, all the four models are capable of reconstructing reasonable glacier surface elevation estimates in most areas.

Comparing our errors to the 2007 measurements reveals that the restitution method is biased towards overestimating the error, supporting the application of conservative error estimates.