List of main symbols
In the observational data
- D aaaa
dataset of surface measurements retrieved during the year aaaa
- DSMaaaa
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 2001
- D f
dataset of surface elevations at the final date, t f. In the studied case of Hurd Glacier, D f is D 2013
- 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, Reference Navarro2009; Bamber and others, Reference Bamber2013; Fretwell and others, Reference Fretwell2013; Lapazaran and others, Reference Lapazaran2013; Morlighem and others, Reference Morlighem2017). 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, Reference Farinotti2017; Fürst and others, Reference Fürst2017). 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, Reference Strozzi, Luckman, Murray, Wegmüller and Werner2002; Sánchez-Gámez and Navarro, Reference Sánchez-Gámez and Navarro2017; Gardner and others, Reference Gardner2018).
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, Reference Lapazaran2013; Berthier and others, Reference Berthier2014, Reference Berthier, Cabot, Vincent and Six2016; Wang and Kääb, Reference Wang and Kääb2015; Zemp and others, Reference Zemp2019; Hugonnet and others, Reference Hugonnet2021). 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, Reference Moholdt, Nuth, Hagen and Kohler2010; Kääb and others, Reference Kääb, Berthier, Nuth, Gardelle and Arnaud2012; Gardner and others, Reference Gardner2013). 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, Reference Navarro, Jonsell, Corcuera and Martín-Español2013; Zemp and others, Reference Zemp2013).
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, Reference Lapazaran2013; Möller and other, Reference Möller, Navarro and Martín-Español2016; Morlighem and others, Reference Morlighem2017). 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, Reference Jonsell, Navarro, Bañón, Lapazaran and Otero2012). Similarly, DSMs are used for both statistical and dynamical downscaling of regional climate model data to the glacier scale (Mölg and Kaser, Reference Mölg and Kaser2011; Jarosch and others, Reference Jarosch, Anslow and Clarke2012; Hannesdóttir and others, Reference Hannesdóttir2015).
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 (Reference Frappé and Clarke2007) 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, Reference Østrem and Brugman1991; Hagen and Reeh, Reference Hagen, Reeh, Bamber and Payne2004; Cogley and others, Reference Cogley2011), comparison of radar and probing data applied to internal radar reflections (e.g. Callens and others, Reference Callens, Drews, Witrant, Philippe and Pattyn2016), inversion techniques applied to internal radar reflections (e.g. Köhler and others, Reference Köhler, Moore, Kennett, Engeset and Elvehøy1997) or inversion techniques applied to several DSMs of a glacier and the digital terrain model of its bedrock (e.g. Välisuo and others, Reference Välisuo, Zwinger and Kohler2017). 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 km2 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−1. 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, Reference Navarro, Jonsell, Corcuera and Martín-Español2013). The maximum ice thicknesses of Hurd Glacier, ~200 m, are found in the accumulation area (e.g. Molina and others, Reference Molina, Navarro, Calvet, García-Sellés and Lapazaran2007). 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, Reference Navarro, Jonsell, Corcuera and Martín-Español2013) and GPR studies since 1999 (Navarro and others, Reference Navarro2009).
(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 2001, D 2007 and D 2013); 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 (Reference Furdada, Pourchet and Vilaplana1999) in the divide between Hurd and Johnsons glaciers.
Our surface restitution method requires two DSMs for Hurd Glacier, DSM2000 and DSM2013, 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 2001 and D 2013 (Figs 2a, c). However, two more sets of surface measurements are employed: one (D 2007, 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.
Black dots represent the datasets of surface elevation GNSS measurements on Hurd Glacier: (a) 308 points in D 2001; (b) 885 points in D 2007 and (c) 7174 points in D 2013. The DSMs at the ends of the period are shown in background: (a) DSM2000 and (c) DSM2013. 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 2001, D 2007 and D 2013 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 (Reference Rodríguez and Navarro2017) (https://doi.org/10.1594/PANGAEA.873067), although D 2001, D 2007 and D 2013 are the corresponding subsets of points within Hurd Glacier (Fig. 2). The surface elevation data in D 2001, D 2007 and D 2013 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 2001 was surveyed by foot with horizontal and vertical accuracies of ±7 cm SDh and ±10 cm SDv. The 2007 survey was performed by snowmobile with accuracies of ±10 cm SDh and ±14 cm SDv. D 2013 was surveyed by foot and snowmobile including photogrammetric restitution at the front of Johnsons Glacier, acquiring final horizontal and vertical accuracies of ±65 cm SDh,v (Rodríguez, Reference Rodríguez2014).
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 SDh,v in both D 2001 and D 2007. However, since D 2013 also includes ablation corrections, we estimate that its measurement error is ±65 cm SDh,v (thus assuming as standard the maximum error of D 2013).
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 (Reference CogleyCogley 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, Reference Huss and Hock2015). 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 2007, b 2007, 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, Reference Cogley2011).
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, Reference Navarro, Jonsell, Corcuera and Martín-Español2013). 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 (Reference Navarro, Jonsell, Corcuera and Martín-Español2013). 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 (Reference Rodríguez-Cielos2016) 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, Reference Oliva2017). 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−1, 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 (Reference Navarro, Jonsell, Corcuera and Martín-Español2013). 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 SDh,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 SDv.
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)
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).
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):
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
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
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.
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−3 for the current year's snow (age is 0) and 520 kg m−3 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−3, the density at which firn becomes glacier ice (where voids cease to form a connected network; Cogley and others, Reference Cogley2011) at the same age as in Furdada and others (Reference Furdada, Pourchet and Vilaplana1999), i.e. in ~24 years, but rises exponentially until year 24 (data from M. Pourchet and J.M. Casas reported in Ximenis, Reference Ximenis2001). 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−3 and $\forall t \ge 1\colon \rho = 0.4707\;t^{0.1852}\;{\rm kg}\;{\rm m}^{{-}3}$
(limited to 900 kg m−3, 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−3; see Section 4.2) for the current year accumulation and ice (~900 kg m−3) 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, Reference Wackernagel2003). 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 2001 and D 2013, thus respectively obtaining DSM2001 and DSM2013. This type of geospatial interpolation has also been applied to generate the SMB maps (Hock and Jensen, Reference Hock and Jensen1999), 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 (Reference Lapazaran, Otero, Martín-Español and Navarro2016), 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 (Reference Lapazaran, Otero, Martín-Español and Navarro2016) 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 (Reference Lapazaran, Otero, Martín-Español and Navarro2016) (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 2001 being D i and D 2013 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.
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 DSM2007
We use the dataset D 2007 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.
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.
Parameters of the errors (in m) when comparing the restituted values using the four restitution models, with the surface elevations measured in D 2007
For each model, minimum error, maximum error, std dev. of the error, average error (restitution bias) and RMS error are given.
(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.
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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2021.104
Acknowledgements
We thank the Scientific Editor, Shad O'Neel, the Associate Chief Editor, Frank Pattyn and the reviewers for their valuable comments and suggestions that helped improving the quality of this manuscript. We also thank Francisco Navarro for his support in providing and refining the SMB source data used in this study. This research has been supported by the Ministerio de Economía, Industria y Competitividad, Gobierno de España (grant no. CTM2017-84441-R).
Author contributions
JL led the study. JL and JO developed the core of the method at its first steps. DM, JL and JO developed the final version of the method. DM developed all the calculations and programmes. CR-B provided SMB data. All the authors contributed to the discussion of the results. JL and DM wrote the initial draft of the paper. JL made the figures. All co-authors contributed to and approved the final manuscript.
Appendix A: Surface mass balance (SMB) on Hurd Glacier 2002–13
See Table 3.
Glacier-wide winter, summer and annual SMB, ELA and accumulation area ratio for Hurd Glacier over the study period, together with their averages and std. dev.s
Appendix B: Implementation of the memory model, M
Section 4.3.2 describes the density model, M, which considers an accumulation memory, with variable density depending on its age (using a DAL). Here, we detail how it is technically implemented.
The application of the density model, M, to a selected point on the glacier surface implies the storage of an accumulation memory specific for this point, using an LIFO (last in, first out) cell-array structure. Each cell contains the accumulation of a certain age. Both the number of cells and their corresponding ages increase by one each year. At the bottom of the array, the ablation of initial mass at this point is memorised, if necessary.
The algebraic structure of the LIFO memory consists of an array with as many cells as memorised ages plus one (cell #0 to cell #n, for a memory of n years). The position and content of each cell in the array, respectively, represent the age and amount of accumulated material (in m w.e.) of this age (age 0 means accumulation at the present year), thus containing only positive or null values. The only exception is the last cell (#n), which contains a null or negative value, representing the total ablation (negative accumulation, in m w.e.) of the initial mass at this point. The content of the cells evolve continuously in time, since accumulation grows during winter seasons and then ablation progresses during summer seasons.
The size of this array grows each year by the addition of a new cell #0 (with zero value), thus incrementing by one the corresponding position (and age) of the pre-existent cells. The present year's accumulation, if positive, is stored in cell #0. However, in cases of net ablation, cell #0 preserves its zero value (no accumulation of present year's snow), and ablation is propagated to the next cell with a non-zero value, diminishing it as follows. If the propagated ablation is larger than the value in the cell, the entire material of this age is ablated, thus storing a zero in the corresponding cell and propagating the remaining ablation to older layers (next cells). If the entire column of accumulation is melted (diminishing to zero the values in cells #0 to #(n − 1)), the remaining ablation diminishes (increments in negative) the stored value in cell #n, thus ablating the initial material, since it is not covered by any memorised accumulation.
Note that although the cell values are not static and could be modified each year, only cell #0 can increment its value (when there is net accumulation during the current year), while the values of the other cells can only stay or decrease (in years with net ablation, deeper and older layers could be ablated if the previous layers are entirely ablated).
Once the array of accumulated materials (amounts and ages) is known for a point and date, we use the DAL (Section 4.2) to convert these accumulated masses into elevation changes, Δz SMB(t), as the summation of the accumulated mass values of each year divided by the corresponding density related to its age. In the case of ablation of the initial mass, it is assumed with ice density.
Appendix C: Error analysis
List of additional symbols
- ɛR(t)
restitution error to a date t
- ɛΩ(t)
component of ɛR(t) due to the wrong estimates of the elevations at the ends of the period, t i and t f, characterised as a straight line between $\varepsilon _\Omega ( {t_{\rm i}} )$
and $\varepsilon _\Omega ( {t_{\rm f}} )$- ɛIi
component of $\varepsilon _\Omega ( {t_{\rm i}} )$
: error of the spatial interpolation of the DSM at t i- ɛMi
component of $\varepsilon _\Omega ( {t_{\rm i}} )$
: propagation of measurement errors to the DSM at t i- ɛIf
component of $\varepsilon _\Omega ( {t_{\rm f}} )$
: error of the spatial interpolation of the DSM at t f- ɛMf
component of $\varepsilon _\Omega ( {t_{\rm f}} )$
: propagation of measurement errors to the DSM at t f- ɛP(t)
component of ɛR(t) due to the wrong prediction of the surface evolution in time. Characterised in a parabolic shape with vertical axis, null at t i and at t f, and maximum, ɛPmax, at the centre of the period
- z M(t C)j
j-th of the elevation measurements at points of the glacier surface at a date t C
- z R(t C)j
restituted elevation at a date t C corresponding to the (x, y) point of z M(t C)j
- D(t C)j
restitution discrepancy, D(t C)j, by subtracting z M(t C)j from the corresponding restituted elevation, z R(t C)j
- ɛD(t C)
standard deviation (estimated by the RMS value) of the overall random variable of restitution discrepancies, D(t C), for the glacier at this date
C1 Method of the error analysis
Section 4.6 describes the main lines of the procedure we use to analyse the restitution error. Here, we detail how it is technically implemented.
The restitution error to a date t, ɛR(t), can be split into two: one, ɛΩ(t), due to the wrong estimates of the elevations at the ends of the period, t i and t f; and second, ɛP(t), due to the wrong prediction of the surface evolution in time.
Both errors can be treated as independent, and can consequently be combined as
C1.1 Propagation of the elevation errors at the ends
The elevation errors at the ends of the period ($\varepsilon _\Omega ( {t_{\rm i}} )$
and $\varepsilon _\Omega ( {t_{\rm f}} )$
), affect the restitution to any date t. In fact, the error ɛΩ(t) is the propagation in time of the errors on point elevations at t i and at t f. These errors mostly depend on the quality and distribution of the elevation data and the method used for spatial interpolation.
The propagation in time of $\varepsilon _\Omega ( {t_{\rm i}} )$
starts from its complete value at t i and diminishes to zero at t f (the elevation errors at t f correspond to another dataset), and vice versa in the case of $\varepsilon _\Omega ( {t_{\rm f}} )$
. Moreover, their combination, ɛΩ(t), would not be larger than the largest of both errors at the ends, ɛΩ(t i) and $\varepsilon _\Omega ( {t_{\rm f}} )$
. Thus, we set a linear interpolation between both values to propagate the error from the ends to any intermediate date t:
where $m = ( {\varepsilon_\Omega ( {t_{\rm f}} ) -\varepsilon_\Omega ( {t_{\rm i}} ) } ) /( {t_{\rm f}-t_{\rm i}} )$
and $n = ( {\varepsilon_\Omega ( {t_{\rm i}} ) t_{\rm f}-\varepsilon_\Omega ( {t_{\rm f}} ) t_{\rm i}} ) /( {t_{\rm f}-t_{\rm i}} )$
.
C1.2 Prediction error
The prediction error, ɛP(t), is the error on the prediction of the point elevation change between t i and t (any date). This would be the only error if the restitution method were applied to points measured, both at t i and at t f, with no errors in elevation (no interpolation error or measurement error). It comes from the errors on SMB, on densities and on the assumptions about the evolution of the dynamic component of the elevation change.
We characterise ɛP(t) in a parabolic shape with vertical axis, being null at t i and at t f, and assuming that it is maximum at the centre of the period, ɛPmax = ɛP((t i + t f)/2). Thus, its equation must be:
where γ = −4 ɛPmax/(t f − t i)2.
C1.3 Best fit
To estimate ɛR(t), we propose a method based on the discrepancies between a set of elevation measurements of the glacier surface (see Section 3.3) at a sequence of intermediate dates, t C, and their corresponding restituted elevations. We consider this set of measurements to be the comparison dataset.
First, we get the restitution discrepancy at each point of the set of surface measurements retrieved at the same comparison date t C. We denote z M(t C)j as the j-th measured elevation of the set of elevation measurements at points of the glacier surface at a date t C. We calculate the restitution discrepancy, D(t C)j, by subtracting z M(t C)j from the corresponding restituted elevation, z R(t C)j:
Then, we take the RMS value of the discrepancies at a date t C as an estimate of ɛD(t C), the std dev. of the overall random variable of restitution discrepancies, D(t C), for the glacier at this date. We prefer using their RMS value instead of their std dev., because if discrepancies are biased (have non-zero mean), the use of their std dev. would become a global underestimation of the discrepancies (small value of ɛD(t C)).
The more random and widespread the distribution on the glacier of the points with elevation measurement at t C, the better the estimation of ɛD(t C). If this distribution is not random and widespread enough, we recommend overestimating ɛD(t C).
From Eqn (C4), we see that the discrepancies are affected by both measurement errors and restitution errors. Since measurement and restitution are independent random variables, ɛD(t C) can be split as
where ɛR(t C) represents the restitution error of the surface points at t C, and ɛM(t C) is the measurement error of the surface data points measured at t C (assumed that they are evenly distributed). Note that when using ɛD(t C) as an estimate of ɛR(t C), the contribution of the measurement errors is to overestimate the restitution error.
The dataset of surface measurements used to establish the discrepancies contain data at a sequence of different intermediate dates t C. However, we aim for a continuous function in time to estimate ɛR(t) at any time within the period. This is carried out following Eqn (C1) and characterising $\varepsilon _{\rm R}^2 ( t )$
by its combination with Eqns (C2) and (C3), thus obtaining a fourth-degree polynomial:
where
The polynomial of Eqn (C6) would be obtained by least-squares fitting of the $\varepsilon _{\rm R}^2 ( {t_{\rm C}} )$
values estimated using ɛD(t C) as an estimate of ɛR(t C) at each of the different dates t C. The system has a total of three unknown parameters. Although we show the system expressions using γ, m and n, we solve it using ɛPmax, ɛΩ(t i) and ɛΩ(t f), due to their smaller sensitivity to computational errors and their physical interpretation.
C1.4 Independent re-estimates of the errors at the ends
The proposed best-fit procedure includes the estimate of both the parabolic prediction error and the linear propagation of the errors at the ends (initial and final DSM). However, the resulting errors at the ends turn out to be dependent on the applied restitution method. We proceed to re-estimate the errors at the ends using a method that depends only on the measurement and interpolation errors of the corresponding DSMs. The resulting $\varepsilon _\Omega ( {t_{\rm i}} )$
and $\varepsilon _\Omega ( {t_{\rm f}} )$
are then propagated in time following Eqn (C2) to obtain ɛΩ(t).
Note that although we re-estimate ɛΩ(t), its previous estimate by the least-squares fitting cannot be suppressed. The use of a complete least-squares fitting for ɛΩ(t) and ɛP(t) together is still required. It should be considered estimating ɛP(t) by subtracting (squared quadratic subtraction, from Eqn (C1)) the values resulting from the ɛΩ(t) deduced by an independent method, from the total restitution errors (using ɛD(t C) as an estimate of ɛR(t C)). However, it is not possible, since the method to re-estimate the errors at the ends aims at giving overestimated errors, so such a subtraction gives impossible results each time the overestimated value of ɛΩ(t) is larger than the restitution error.
Below we describe the procedure to create a new ɛΩ(t) by re-estimating $\varepsilon _\Omega ( {t_{\rm i}} )$
and $\varepsilon _\Omega ( {t_{\rm f}} )$
.
To re-estimate the errors at the ends, $\varepsilon _\Omega ( {t_{\rm i}} )$
and $\varepsilon _\Omega ( {t_{\rm f}} )$
, we follow Lapazaran and others (Reference Lapazaran, Otero, Martín-Español and Navarro2016). It implies splitting the error at each DSM, ɛΩ, into two independent errors: ɛI, the error of the spatial interpolation of the DSM; and ɛM, the measurement error at the DSM, which is the result of the spatial propagation of the errors of the elevation measurements involved in the interpolations of the DSM. For the DSMs at the ends, t i and t f, it results
To estimate the interpolation error at the initial and final surfaces, ɛIi and ɛIf, we follow the modified technique introduced by Lapazaran and others (Reference Lapazaran, Otero, Martín-Español and Navarro2016). This technique assesses the bias and std dev. of the interpolation error at any point, as a function of its distance to the nearest data point.
To propagate the measurement errors to the interpolation point to obtain ɛMi and ɛMf, we also follow Lapazaran and others (Reference Lapazaran, Otero, Martín-Español and Navarro2016). Thus, since the interpolated elevation is obtained as a weighted combination of the elevations at some neighbouring data points, the error propagation from the data points to the interpolation point is calculated using the same weighted combination, but now applied to the errors instead of the elevations. In almost all cases, ɛM is a known and unique value for the entire dataset (i.e. one at t i and another at t f). In such cases, the outcome of the propagation is also this value.
The previous descriptions facilitate the estimate of point dependent errors for ɛI, ɛM and ɛΩ. To characterise the error of the restitution process for the entire glacier, we evaluate both ɛI and ɛM at the nodes of a regular grid covering the entire glacier under study, at the dates with surface measurements (t i and t f). Thus, following Eqn (C7), ɛΩ can be estimated at every grid node, both at t i and at t f. The RMS value of the errors at the grid nodes characterises the error for the whole glacier:
where N is the number of grid nodes and the subscript k represents the k-th grid node.
Then, ɛΩ averaged for the entire glacier can be estimated at any time t, following Eqn (C2), using the results of Eqn (C8), to obtain $\varepsilon _\Omega ( t )$
.
Although the method used to re-estimate the error ɛΩ(t) aims to be conservative (i.e. showing an error overestimating trend), it could be that the values of ɛΩ at the ends, $\varepsilon _\Omega ( {t_{\rm i}} )$
or $\varepsilon _\Omega ( {t_{\rm f}} )$
, resulting from Section C1.3 were larger than their re-estimations. In such a case, with the criterion that it is better overestimating an error than underestimating it, we propose the use of the largest value at each extreme.
C1.5 Restitution error
The final estimate of ɛR(t) is given by the combination, following Eqn (C1), of ɛP(t) resulting from Section C1.3 and ɛΩ(t) resulting from Section C1.4. It is a conservative estimate of ɛR(t), in which ɛP(t) depends on the applied restitution method, but ɛΩ(t) does not.
C1.6 Errors of the SMB and dynamic components
The proposed method starts evaluating Δz SMB(t f), which, added to z i, is taken as part of the elevation at t f which is only explained by the SMB. Then, the dynamic component, Δz d(t f), is obtained as its difference to z f (Eqn (2)). However, such a procedure implies that the error produced in the estimate of the SMB component at t f, Δz SMB(t f), is included in the dynamic component at t f, Δz d(t f), and later propagated to any other date t. Therefore, the errors of both components are dependent, both affected by the errors in SMB, although in opposite signs. However, the aim of this study does not include the separate evaluation of the errors of both components of the elevation change.
C2 Error estimate for the case study on Hurd Glacier
Section 5.3 summarises the main results of the error analysis when applied to the case study on Hurd Glacier. Here, we detail this error analysis and give the complete set of results.
As stated in Section C1 of Appendix C, the estimates of the errors associated with the restitution at any date t for the four models are implemented following the described error model of a straight line for ɛΩ(t) (Eqn (C2)) and a parabola for ɛP(t) (Eqn (C3)), both quadratically joined as in Eqn (C6) to characterise the squared restitution error. Although ɛΩ(t) and ɛP(t) are both calculated by least-squared fitting of squared restitution errors, we re-estimate ɛΩ(t) using a conservative method only dependent on the measurement and interpolation errors of the initial and final DSMs.
To obtain a cloud of squared restitution errors with which we could generate the best fit by least-squares, we use a set of surface data points at intermediate dates, well distributed, both spatially on Hurd Glacier and in time along the period (Section 3.3). We use Eqn (C4) to compare the set of surface measurements with the results of the restitution process at each point and date, thus obtaining a set of discrepancies D(t C)j, each one at its date, t C. Then, we join discrepancies in batches by their date of measurement, t C. Since we detect that the discrepancies are not bias-free, we characterise the restitution error of each batch of discrepancies as the RMS value of this batch, instead of its std dev..
In Figure 8, we show the set of estimates of squared restitution errors at the different dates, t C, with batch of discrepancies, after the calculations explained above. To extend the restitution error to any date t, we follow Eqn (C6), once the best values of the parameters ɛPmax, $\varepsilon _\Omega ( {t_{\rm i}} )$
and $\varepsilon _\Omega ( {t_{\rm f}} )$
are obtained by mean-squared fitting.
Estimates of the squared restitution errors for the four restitution models. Dots represent squared restitution errors at the batches of surface elevation measurements. The lines are the best fit of Eqn (C6), as described in Section C1.3, used to extend the restitution error at any date within the period. Dates in the X-axis indicate 1st of December for the corresponding year.
Splitting the restitution error into ɛΩ(t) and ɛP(t), due to the values of the selected best-fit parameters, ɛPmax, $\varepsilon _\Omega ( {t_{\rm i}} )$
and $\varepsilon _\Omega ( {t_{\rm f}} )$
, as shown in Table 4, gives us important information. Note that errors at t i are estimated to be (Table 4) between 0.70 and 0.82 m, and errors at t f fall between 0.31 and 0.61 m. In practice, it shows that the glacier surface at t i is generated with larger error than that at t f (due to the number and distribution of measured points, their measurement errors and the interpolator used).
Parameters of the error estimates (in m) for the four models as resulting from the best fit of the squared restitution errors obtained when comparing the results of each model with the set of surface measurements described in Section 3.4
However, theoretically, $\varepsilon _\Omega ( {t_{\rm i}} )$
and $\varepsilon _\Omega ( {t_{\rm f}} )$
must not depend on the applied restitution method. To conservatively unify their values, we re-estimate them using an independent method described in Section C1.4. The resulting interpolation errors are ε Ii = 1.25 m and ε If = 0.38 m. Combining them using Eqn (C7) with the measurement errors described in Section 3.1 for the datasets D 2001 and D 2013, respectively ε Mi = 0.25 m and ε Mf = 0.65 m, we obtain ε Ω(t i) = 1.28 m and ε Ω(t f) = 0.75 m. As predicted, the use of such a conservative method gives larger errors than those estimated using any model.
Figure 9 shows the comparison of the whole error estimates for the four models. Looking at the prediction errors, both in Table 4 and in Figure 9, we find that the largest error estimates are for nM_BP (ε Pmax = 0.97 m), and the smallest, for M_MAP (ε Pmax = 0.84 m). In decreasing order of maximum values of the restitution error (somewhat before the middle of the period, in 2006), we find 1.41, 1.38, 1.38 and 1.34 m for the respective models nM_BP, nM_MAP, M_BP and M_MAP (being of 2.27 m for LIN, not shown).
Solid lines show, for each model, the estimated glacier-wide restitution error at any date within the period. Dotted parabolas show the estimated prediction errors of the models, while the straight-dashed line is the propagation in time of the DSM errors at the ends of the period.
It may seem contradictory that nM_BP presents the largest error estimates (Fig. 9) while it has been the best model when comparing restitution vs 2007 measurements (Table 1 and Fig. 7). However, it just means that the procedure for estimating errors generates the largest ɛP(t) at nM_BP. Note that the error estimates of the different models only differ in their prediction error, since we use a common ɛΩ(t) for all the models. This prediction error estimate has been inferred by splitting the restitution errors (at the batches of surface elevation measurements) in ɛΩ(t) (not used) and ɛP(t), during the least-squared fitting. Thus, the resulting ɛP(t) depends on how the best fit has distributed such a split.
Open access
