ORIGINAL RESEARCH article
Reconstruction of Past Glacier Changes with an Ice-Flow Glacier Model: Proof of Concept and Validation
- 1Institute of Geography, University of Bremen, Bremen, Germany
- 2Institute of Hydrology and Water Resources Management, Leibniz Universität Hannover, Hannover, Germany
- 3Department of Atmospheric and Cryospheric Sciences, University of Innsbruck, Innsbruck, Austria
- 4MARUM - Center for Marine Environmental Sciences, University of Bremen, Bremen, Germany
Estimations of global glacier mass changes over the course of the 20th century require automated initialization methods, allowing the reconstruction of past glacier states from limited information. In a previous study, we developed a method to initialize the Open Global Glacier Model (OGGM) from past climate information and present-day geometry alone. Tested in an idealized framework, this method aimed to quantify how much information present-day glacier geometry carries about past glacier states. The method was not applied to real-world cases, and therefore, the results were not comparable with observations. This study closes the gap to real-world cases by introducing a glacier-specific calibration of the mass balance model. This procedure ensures that the modeled present-day geometry matches the observed area and that the past glacier evolution is consistent with bias-corrected past climate time series. We apply the method to 517 glaciers, spread globally, for which either mass balance observations or length records are available, and compare the observations to the modeled reconstructed glacier changes. For the validation of the initialization method, we use multiple measures of reconstruction skill (e.g., MBE, RMSE, and correlation). We find that the modeled mass balances and glacier lengths are in good agreement with the observations, especially for glaciers with many observation years. These results open the door to a future global application.
Compared to the Greenland and Antarctic ice sheets, the ice mass stored in glaciers seems negligibly small (<1%) (e.g., Church et al., 2013), but glaciers contributed between 25 and 30% of the observed global mean sea-level rise during 1993–2015 (Zemp et al., 2019). Along with ice sheet mass loss and thermal expansion of sea water, glaciers will continue to be a major sea-level rise contributor in the 21st century (Slangen et al., 2017; Oppenheimer et al., 2019; Marzeion et al., 2020). Independent from the Representative Concentration Pathway (RCP), glacier mass loss may exceed the contributions from Greenland or Antarctica ice sheets until 2100 (Oppenheimer et al., 2019). However, glacier changes also have an impact on local and regional scales: mountain glaciers affect the water availability in many regions of the world (e.g., Kaser et al., 2010; Biemans et al., 2019; Immerzeel et al., 2020; Viviroli et al., 2020) and both retreating and advancing glaciers can lead to increased geohazards (see Kääb et al., 2005, for an overview). It is consequently relevant to improve the knowledge of past and future glacier mass change (Marzeion et al., 2012).
To estimate glacier mass change on a global scale, in situ measurements of mass and length changes (e.g., Leclercq et al., 2014; Zemp et al., 2015), remote-sensing techniques (Jacob et al., 2012; Gardner et al., 2013; Bamber et al., 2018; Wouters et al., 2019), or mass balance modeling driven by climate observations (e.g., Radić and Hock, 2011; Giesen and Oerlemans, 2012; Marzeion et al., 2012, Giesen and Oerlemans, 2013; Marzeion et al., 2014; Radić and Hock, 2014; Huss and Hock, 2015; Maussion et al., 2019) can be used. However, most of these model-based studies limit their application to the recent past (ignoring past glacier geometry change) or to the future. Nevertheless, glacier mass change in the past should not be neglected: a large fraction of future glacier mass loss (36
So far, only one model was able to provide global estimates of glacier mass changes over the course of the entire 20th century, including a simple representation of glacier geometry change (Marzeion et al., 2012). For more complex models (e.g., Huss and Hock, 2015; Maussion et al., 2019), this process is impeded by initialization issues, as more and more variables need to be initialized. That is, flow-line models require input data along the coordinates of the flow lines (e.g., bed topography, surface elevations, and/or widths) and thus more complex initialization methods are required. Zekollari et al. (2019) added an ice-flow model to Huss and Hock (2015) and developed an automated initialization for glaciers in 1990 (i.e., prior to the inventory date) in order to avoid spin-up issues. The ensemble approach from Eis et al. (2019) can be used for an initialization further back in time (e.g., 1850), but the larger time difference to the inventory date leads to nonunique solutions of the initial glacier geometry. Hereafter, we refer to glacier geometry (i.e., glacier outlines and surface elevation of a glacier at a given point in time) as the “state” of the glacier.
The approach from Eis et al. (2019) works as follows: first, a large number of physically plausible glacier states for a given year in the past is generated, using a spin-up run of the model forced by random climate representatives of the conditions around the year of initialization. In order to create a large set of different states, we additionally vary a temperature bias β. From this set of possible past glacier states, a subset of equidistant and approximately uniformly distributed states is selected and called “glacier candidate states.” All these candidate states are then modeled forward to the present-day and evaluated by comparing the result of the forward runs to the present-day glacier states. In order to separate uncertainties of the initialization method alone from all other sources of uncertainty, Eis et al. (2019) relied on synthetic experiments. These synthetic experiments were generated by running the same glacier model forward such that a perfectly known glacier evolution until the present day was produced (including boundary conditions). The synthetic experiment state at the present day was then used to reconstruct the synthetic experiment state in 1850.
While this procedure allows an exact evaluation of the uncertainties of the method alone, it does not allow an external validation, e.g., against in situ measurements of glacier mass balance or length changes. The results from Eis et al. (2019) thus do not provide information about actual past glacier changes, because the synthetic experiment states in 2000 can be different to real states. Consequently, the reconstructed evolution from 1850 onward does not necessarily correspond to reality either, and uncertainties in a real-world application will necessarily be larger than in the synthetic experiments.
In order to allow the application to real-world cases, it is necessary to modify the method of Eis et al., 2019 by introducing a glacier-specific calibration, which is presented in this study. This modification allows a validation of Eis et al. (2019) based on real glacier observations. We first introduce all relevant features of the Open Global Glacier Model (OGGM) in The Open Global Glacier Model. Next, we present the glacier-specific calibration of the mass balance model of the OGGM by introducing the setup of calibration runs (see Glacier-Specific Calibration of the Mass Balance Model). The performance of these calibration runs is assessed in Optimized Mass Balance Offsets. We then apply the initialization method from Eis et al. (2019) by using the results of the calibration runs to 517 glaciers with mass balance observations or length records. The results of the initialization are shown in Reconstructions in the Year 1917. Afterwards, we present the validation of the initialization method from Eis et al. (2019) by using glacier observations (see Validation With Glacier Observations). To this end, a statistical analysis based on mass balance observations provided by the World Glacier Monitoring Service (WGMS, 2018) can be found in Mass Balance Validation. Additionally, our method is validated based on glacier length change observations from Leclercq et al. (2014) (see Glacier Length Validation). Finally, we summarize the results and discuss the limitations of the method, including future applications in Discussion and Conclusion.
The Open Global Glacier Model
The open-source numerical framework Open Global Glacier Model (OGGM; Maussion et al., 2019) is able to model the evolution of all glaciers worldwide. OGGM utilizes the outline of any glacier, which can, e.g., be obtained from the Randolph Glacier Inventory (RGIv6.0; Pfeffer et al., 2014). A digital elevation model (DEM), covering the glacier outlines, is automatically downloaded and subsequently interpolated to a local grid. The resolution
In order to allow a glacier to also grow downstream of the current glacier geometry, the OGGM calculates a downstream line. A routing algorithm, minimizing the distance between the glacier terminus and the border of the map, including the total elevation gain, ensures that the glacier will follow the valley floor. Cross-sections at each of the downstream line grid points are obtained by fitting a parabola to the valley topography. As this part of the glacier is not covered by ice at the inventory date of the glacier outline, the bed topography along the downstream line can directly be derived from the topography file. Consequently, all information along the downstream line is constrained more narrowly than along the flow line, as we do not rely on the ice thickness inversion at this part of the glacier.
Utilizing monthly temperature and precipitation data, mass balances at the flow-line grid points are calculated. A number of different sources of climate data can be used to force the mass balance model, e.g., past climate observations (gridded), reanalyses, projections, or randomized time series. The evolution of the glacier geometry under the specific climate forcing is calculated by a dynamical flow-line model, which solves the shallow-ice approximation. More information about the OGGM can be found under http://docs.oggm.org and in Maussion et al. (2019).
Input Data and Default Model Parameters
The glacier outlines used in this study are taken from the Randolph Glacier Inventory (RGI v6.0, region 11; Pfeffer et al., 2014). For each of the glacier outlines, a transverse Mercator projection, centered on the glacier, is defined in order to preserve map projection properties (e.g., area, distances, and angles). Next, the topographical data are automatically downloaded by the OGGM. The acquisition dates of the DEMs vary from 2000 (for SRTM) to 2009 and roughly match that of the RGI. The climate dataset we use for this approach is the Climatic Research Unit (CRU) TS v4.01 dataset (Harris et al., 2014, released on September 20, 2017). It is a gridded dataset at 0.5° resolution that covers the period from 1901 to 2016. Further downscaling to a resolution of 10 arc minutes is achieved by applying the 1961–1990 anomalies to the CRU CL v2.0 gridded climatology (New et al., 2002), see also Maussion et al. (2019) for more details. Monthly temperature and precipitation time series are extracted for each glacier from the nearest CRU CL v2.0 grid point and then converted to the local temperature according to a temperature gradient (default: −6.5 K km−1). As the initialization procedure requires climate data of the 15 preceding years of the initialization year
Identical to the study of Maussion et al. (2019) and Eis et al. (2019), only the mass balance model is calibrated while the creep parameter A and the sliding parameter
Glacier-Specific Calibration of the Mass Balance Model
Eis et al. (2019) showed that the use of synthetic experiments provides many advantages (e.g., the understanding of the limitations and errors of the developed method itself through the isolation of the analysis from real-world uncertainties), but it also causes one main disadvantage. The geometry of synthetic experiment states at the inventory date usually does not coincide with the geometry of the observed states. Thus, the reconstructed time series does not correspond to real past glacier states, which implies that a validation of the method with glacier observations is not possible.
In this study, we aim to reduce the geometry mismatch between the model results and the real present-day states by introducing the so-called calibration runs. This improvement allows the application of the initialization method from Eis et al. (2019) to real-world cases followed by the validation of the reconstructed time series with glacier observations. To this end, we introduce a mass balance offset
Before motivating the choice to introduce and optimize the mass balance offset
With the help of the calibration runs, we iteratively obtain an optimal mass balance offset
To this end, we generate a glacier state for the year of the initialization
FIGURE 1. Workflow of the calibration runs and the ensemble approach from Eis et al. (2019). (A,B) Calibration runs; the gray lines indicate the iteration procedure and the brown line the result obtained with the optimal mass balance offset
As the mass balance offset will be subtracted from the modeled mass balance, the mass balance offset is increased if the area of
This procedure is repeated until the maximum number of iterations
The optimal mass balance offset
We apply our method to all glacier outlines from the RGI (version 6) that can either be linked to mass balance observations provided by the World Glacier Monitoring Service (WGMS, 2017) or to the glacier length records from Leclercq et al. (2014). We exclude all Antarctic and Subantarctic glaciers (RGI region 19), because no CRU data are available for this region. Additionally, we exclude glaciers that are strongly connected to the Greenland ice sheet (connectivity Level 2 in Rastner et al., 2012). Because of the complex dynamical response of water-terminating glaciers, they are excluded as well. We detect them by using the corresponding RGI attribute. This leads to a total sample size of 519 glaciers. For 253 glaciers, there are mass balance time series, and for 317 glaciers, there are length records (51 glaciers have both). The identification of glaciers in the Leclercq dataset in the RGI is done as in Parkes and Goosse (2020). In our study, we use the two suggested types of matches: the positive matches (the (lon,lat) pair from the Leclercq data lies within the RGI outline or within rounding distance to the outline, and the areas in both datasets are similar to each other) and the best effort matches (the glacier cannot be identified uniquely, as, e.g., the (lon,lat) pair is on the border between two glaciers and none or both glaciers satisfy the area criteria). During the calibration runs, two glaciers failed due to processing errors, and thus, we were able to initialize the glacier states in
Optimized Mass Balance Offsets
We perform the calibration runs for all tested glaciers as described in Glacier-Specific Calibration of the Mass Balance Model and obtain for each of the glaciers an optimal mass balance offset
FIGURE 2. (A) Distribution of the derived optimal mass balance offsets
Figure 2B demonstrates that
Reconstructions in the Year 1917
We use the calibration runs to identify optimal mass balance offsets for each of the glaciers and apply the initialization method from Eis et al. (2019) to reconstruct their states in year 1917. The derived mass balance offsets
with the uncertainty
FIGURE 3. Results for the Hintereisferner (Ötztal, Austria). Evaluation of the glacier candidates by fitness function based on geometry differences to (A,B) the OGGM state
TABLE 1. Glacier length, area, and volume of the median state in 1917 and 2003 of Hintereisferner (Ötztal, Austria).
TABLE 2. Relative differences (in
Validation with Glacier Observations
The uncertainty analysis in Eis et al. (2019) showed that the uncertainty of our method reduces quickly over time. As a result, a meaningful validation of the method, which also covers the time period with the largest uncertainties, should ideally be performed with observations starting before or at the beginning of the initialization year (at the beginning of the 20th century). Remote-sensing techniques are available for a relatively short time period and are consequently not of sufficient length for validation purposes. More suitable are direct observations of glacier mass change, which start to increase in number from the year 1950 onward. Here, we use glacier-wide mass balance observations provided by the World Glacier Monitoring Service (WGMS, 2018). The Fluctuations of Glaciers (FoG) database contains annual mass balance values for several hundreds of glaciers worldwide, but only 40 reference glaciers have record longer than 30 years. Thus, we additionally include glacier length change observations to our validation, as they date back furthest in time. They are usually obtained through annual measurements of the terminus position, but they can also be reconstructed from geomorphologic landforms (e.g., moraines) and biological evidence (e.g., overridden trees and lichens) (e.g., Bushueva et al., 2016) or from historical paintings, photos, and early maps (e.g., Nussbaumer and Zumbühl, 2012). The dataset of Leclercq et al. (2014) includes long-term glacier length records that start no later than 1950 and span at least 40 years. The temporal frequency and spatial accuracy of the dataset increase rapidly from 1900 onward (Leclercq et al., 2014).
As the initialization method results in nonunique solutions, we use the median state
FIGURE 4. Distribution of the reconstructability values of all glaciers used in this study. (A) The values from the glaciers used for the validation with the WGMS dataset and (B) from the glaciers used for the validation with the Leclercq et al. (2014) dataset.
The median state is then used as initial condition for a forward run, using the past climate from 1917 until the RGI date
The validation is performed based on five different error criteria: the model error, the mean bias error (MBE), the mean absolute error (MAE), the root mean square error (RMSE), and the correlation r. A summary of the performance can be found in Table 3. When referring in the following to one of the mass balance errors, a subscripted
TABLE 3. Summary of the comparisons of modeled and measured mass balances and glacier length records.
Mass Balance Validation
Before presenting the mass balance validation in detail, we show the results of three different example glaciers and the corresponding observations from WGMS (2018) in Figure 5. Figure 5A shows an example with an almost perfect match to the observation, whereas the other two figures show two examples with a bad performance. In the example in Figure 5B, mass balances are well matched until 1975, but afterwards the modeled mass balances are high compared to the observation. The glacier shown in Figure 5C fits the observation on average, but the modeled amplitude is far too small. Here, the smaller variability in the modeled mass balances originates from the calibrated precipitation correction factor
FIGURE 5. Reconstructed annual mass balances from 1917 onward of three example glaciers. Figures from all glaciers used in this study can be found in the supplement. (A) RGI60–11.00804, region 11 (Central Europe), (B) RGI60–13.06359, region 13 (Central Asia), and (C) RGI60–10.05085, region 10 (North Asia). The lines are colored by their fitness values. The red line with dots indicates the observed annual mass balances from WGMS (2018). The dots mark years with observations. The vertical dashed line marks
Mass Balance Model Error
The mass balance model error is calculated for each of the 5,325 pairs of annual modeled and measured mass balances. For each pair, the measured mass balance is subtracted from the modeled mass balance at the same year.
Table 3 shows that the mean over all mass balance model errors is with a value of −0.13 m w.e. close to zero, and also, the standard deviation (0.79 m w.e.) indicates a good fit between modeled and observed mass balances. Figure 6 shows a binned scatter plot of the mass balance model errors over time and thus provides information about their temporal distribution. In addition, the temporal distribution of the mass balance measurements (Figure 6A) and the distribution of the mass balance model errors (Figure 6C) are shown. With an increasing number of measurements since 1960, the counts raise, especially in the middle range between −1 and 1 m w. e. Most values are in close proximity to zero, but single errors reach values up to
FIGURE 6. Mass balance model error. (A) Histogram of mass balance measurements per year. (B) Binned scatter plot of model errors over time. The colors indicate the counts of errors in each bin. The numbers in the upper left corner indicate the total number of errors and the correlation with the time. (C) Distribution of the modeled mass balance errors.
Mean Bias Error (MBEmb)
The mean bias error (also referred as bias, e.g., in Marzeion et al. (2012)) is calculated for each of the 253 glaciers with mass balance measurements:
where n is the total number of years with mass balance measurements for the specific glacier. Thus, the MBEmb provides information about the performance of individual glaciers. Figure 7A shows that the MBEmb is small for the majority of the glaciers. This is also reflected by the small mean (0.18 m w.e.) and median (0.14 m w.e.) values. Nevertheless, a few outliers range up to 2 m w.e. The distribution is slightly positive skewed, indicating that the mass balance is marginally overestimated by the model on average.
FIGURE 7. Histogram of (A) MBEmb (mm w. e.), (B) MAEmb (mm w. e.), and (C) RMSEmb (mm w. e.). Each histogram has a sample size of Nmb = 253.
Mean Absolute Error (MAEmb) and RMSEmb
The root mean square error (RMSEmb; Eq. 3) and the mean absolute error (MAEmb; see Eq. 4), calculated for each glacier as well, do not consider the direction of the error. Thus, they provide more information about the magnitude of the errors than the MBEmb:
The RMSEmb gives a relatively strong weight to large errors since the errors are squared before they are averaged and thus the values are larger than the ones of the MAEmb (see Figures 7B,C; Table 3). Their distributions are quite similar to each other, and also, their mean values are in close proximity to each other (0.60 m w.e. for MAEmb and 0.71 m w.e. for RMSEmb). Also, their standard deviation does not differ significantly. Both measures show that the magnitude of the error is small for most glaciers. Even the few outliers fit the observation with an RMSE of about 2 m w.e. well.
Beside the different types of errors, we calculate the correlation between the modeled and measured mass balances over all glaciers. The negatively skewed distribution of the correlation coefficients in Figure 8A shows that for most glaciers, the correlation is high. Nevertheless, for some single glaciers, the modeled mass balances do not correlate with the measurements. Some glaciers even have a negative correlation, which distinctly influences the mean value across all glaciers. Figure 8B shows that there is a relation between the correlation coefficients
FIGURE 8. Correlation coefficients in relation to number of mass balance measurements per glacier. Red dots or bars mark glaciers with more than 30 years of mass balance measurements and blue marks glaciers with less than 30 years. (A) Histogram of the correlation coefficients of all glaciers. (B) Scatter plot of the correlation against the number of mass balance measurements. A large dot indicates a correlation value that is significant at the 5% level and a small dot a not significant value. The horizontal, black line indicates the mean value over all glaciers, the orange line the average correlation weighted number of mass balance measurements per glacier, and the dotted, black line the median value. (C) Histogram of number of mass balance measurement per glacier. Note that the sample size is N = 250, because the WGMS dataset contains three glaciers with constant mass balance series, and due to the zero variance, it is consequently not possible to calculate the correlation.
Glacier Length Validation
In addition to the validation with the mass balance time series from the WGMS dataset, in this section, we compare modeled and observed glacier length changes. As before, we will use the modeled glacier length from the median state of the reconstruction. A clear advantage of glacier length fluctuations is the better temporal coverage compared to the mass balance observations. Some of the length observations even start prior to 1850. Eis et al. (2019) showed that largest uncertainties of the reconstruction method occur at the beginning of the initialization (in this study 1917) and consequently a validation with many early observations is particularly helpful. We transform the glacier length changes from Leclercq et al. (2014) into respective glacier length values. To this end, we first determine all glacier length changes relative to the RGI inventory date
In some cases,
FIGURE 9. Reconstructed glacier lengths from 1917 onward of three example glaciers. Figures from all glaciers used in this study can be found in the supplement. (A) RGI60–11.02584, region 11 (Central Europe), (B) RGI60–05.11966, region 05 (Greenland), and (C) RGI60–11.00090, region 12 (Caucasus and Middle East). The lines are colored by their fitness values. The red line with dots indicates the observed glacier length from Leclercq et al. (2014). The dots mark years with observations. The vertical dashed line marks
Table 3 shows a summary of the performance of all tested glaciers described by the length error, the mean bias error (MBEL), the mean absolute error (MAEL), the root mean square error (RMSEL), and the correlation rL.
Similar to the mass balance model error in Mass Balance Validation, we calculate the length error for each of the 10,420 pairs of modeled and observed glacier length. In opposition to the other errors (MBEL, MAEL, RMSEL, and rL), the length error is not calculated per glacier, but every data point contributes in equal proportions.
Figure 10B shows the distribution of the length errors over time. Most of the values range between ±1 km (see also Figure 10C), and also the small mean value over all length errors of −0.007 km indicates a good agreement with the observations. Nevertheless, a lot of outliers range up to ±3 km and very few values have a length error about ∼6 km (see Figure 10B). This also influences the standard deviation which has a value of 0.62 km. A correlation of the length errors with time is not identifiable (r = −0.118), and the data coverage is good over the whole time period (see Figure 10A).
FIGURE 10. Glacier length error. (A) Histogram of glacier length observations per year. Note that the x-axis is displayed in time difference to
Mean Bias Error (
The mean bias error (
FIGURE 11. Histograms of (A) the mean bias errors (km), (B) the mean absolute errors (km), (C) the root mean square error (km), (D) the mean bias errors (
Mean Absolute Error (
) and RMSEL
The MAEL and RMSEL provide more insights about the magnitude of the errors than the MBEL, as negative and positive values cannot be compensated. Thus, the mean overall MAEL and RMSEL values are with values of 0–35 km (for MAEL) and 0.45 km (for RMSEL) noticeable higher than the ones from the length errors or the MBEL (see Table 3). The standard deviations are slightly smaller (0.31 km for MAEL and 0.4 km for RMSEL), but this can be explained due to the reduced range of possible values (only positive). The histograms in Figure 11 show that the majority of the errors is small. Considering, e.g., the MAEL, then half of the errors are smaller than 0.26 km, which indicates a very good fit between the modeled and observed glacier lengths. Outliers range again up to ∼2 km for MAEL and up to ∼2.4 km for RMSEL, but they occur very rarely. In general, the RMSEL errors are a little bit higher than those of the MAEL which is due to the relatively strong weight to large errors during the RMSE calculation. When considering the errors relative to the glacier lengths, the majority of the MAEL (in %) and RMSEL (in %) values are smaller than 25%. The 95th percentiles are located in both cases under 50% (24% for the MAEL and 34% for the RMSEL), which confirm the good fit.
The results from the correlation analysis of the modeled and observed glacier lengths, calculated for each glacier individually, do not indicate a good agreement for all tested glaciers. The distribution of the individual correlation coefficients for each glacier is slightly bimodal (see Figure 12A). Most correlation values range from 0.7 to 1, but 87 out of 302 glaciers have a negative correlation (some of them close to −1). In contrast to the correlation analysis from the mass balances (see Figure 8B), not all negative correlation values can be explained by a low number of observations or because the correlation is not significant at the 5% level. Figure 12B shows that 12 glaciers with more than 30 length observations have a negative correlation, and four of them are even in close proximity to −1.
FIGURE 12. Correlation coefficients in relation to total number of glacier length observations per glacier. Red dots or bars mark glaciers with more than 30 years of measurements and blue marks glaciers with less than 30 years. (A) Histogram of the correlation coefficients of all glaciers. (B) Scatter plot of the correlation against the number of glacier length observation. Horizontal, black line indicates the mean value over all glaciers, the orange line the average correlation weighted by the number of observations per glacier, and the dotted, black line the median value. (C) Histogram of the number of glacier length observations per glacier. Note that the sample size is N = 302, because for 13 glaciers, the modeled glacier length is constant, and due to the zero variance, it is consequently not possible to calculate the correlation.
An example of one of these outliers can be found in Figure 9C. The reconstructed state in 1917 is far too small, requiring that the glacier needs to advance continuously to reach the present-day state. In contrast to that, the observed glacier retreated over the entire period, which leads to a strong negative correlation. In the figure, one can see that glacier candidates exist which would have reflected the observation very well. However, these candidates have high fitness values, which could, e.g., stem from a bad representation of the OGGM state or from a poorly chosen parameter setup of the mass balance model or dynamical model for individual glaciers.
Nevertheless, the histogram and the scatter plot in Figures 2B, 12A show that a majority of the glaciers with negative correlation values can be attributed to a low number of length observations. The negative correlation values largely influence the mean value (0.36), for what reason the mean weighted over the number of observations per glacier is with a value of 0.65 more representative. This is also underlined by the even higher median value (0.73).
Discussion and Conclusion
This study extends Eis et al. (2019) with the addition of mass balance calibration runs, to adapt the method to real-world use cases. After the standard calibration of the mass balance model, we additionally perform a glacier-specific calibration of the mass balance offset
Water-terminating glaciers are excluded from this study, because a dynamical forward model for tide-water glaciers would be required, which is still under development and validation. For their identification, the attribute from the RGI database is used, but this characteristic only refers to the inventory date
We validate our initialization method by comparing modeled glacier changes to mass balance observations (WGMS, 2018) and to length records (Leclercq et al., 2014). For this purpose, different types of errors (model error, mean bias error, mean absolute error, root mean square error, and correlation) are taken into account. As the initialization can result in nonunique solutions, we use the median state of the 5th percentile for the validation. This is the best representative for the possibly large set of acceptable glacier states (see Eis et al., 2019, for more details), which is additionally justified by the high reconstructability values for glaciers in this study. Unfortunately, this also means that glaciers with a low reconstructability are underrepresented in our study. Eis et al. (2019) showed a negative correlation between the reconstructability of a glacier and its slope. We speculate that this could be one reason for the high reconstructability values in our study, as observations are biased towards the accessibility of a glacier. Another reason could be the positive correlation between the response time of a glacier and its reconstructablity. The considered time interval in this study from 1917 to the present day does not seem to exceed the response time by a large amount for the majority of the investigated glaciers.
Altogether, the validation shows a good agreement between the modeled glacier states and the observations. Considering the mass balance observations, a mean model error (calculated over each pair of modeled and measured mass balances) of 0.13
Regardless of all these different sources of uncertainty, we can summarize that both datasets, the mass balance observations from WGMS (2018) and the glacier length records from Leclercq et al. (2014), are well matched by the modeled glacier states and thus the initialization method from Eis et al. (2019) is verified for the considered time period.
Applications to longer time scales, e.g., over the last centuries/millennia, would in theory be possible under the condition that the climate is well known. However, the reconstructablity value decreases in time and will thus be a limitation factor for this type of application. Taking into account more large-scale observational data (e.g., multitemporal outlines) could improve the reconstructability further and thus would allow a reconstruction further back in time.
Future work will include a global application, which would be the first global glacier reconstruction using a flow-line model. However, a global run would have a very high computational cost, meaning it will likely take multiple months on our small cluster. Computational costs could potentially be reduced by decreasing the model resolution or the number of the glacier candidates. Additional challenges of global reconstructions originate from the poor representation of glacier complexes and ice caps in flow-line models, which stem from their nontrivial geometries (Maussion et al., 2019). Furthermore, the mentioned problems regarding water-terminating glaciers need to be solved beforehand. Such a global reconstruction would also neglect disappeared glaciers, which are uncharted today (Parkes and Marzeion, 2018). The missing but mandatory present-day outlines make it impossible to reconstruct them, and only multitemporal outlines on a large scale could provide a remedy.
Furthermore, a similar approach to the calibration runs could be used in order to optimize other fundamental dynamical parameters in the OGGM, to improve the behavior of the dynamical model in general. For this purpose, it would be sufficient to consider a much shorter time period than in this study (as it was done, e.g., in Zekollari et al., 2019). The sliding parameter and creep parameter, which in the OGGM by default are not glacier-specific, are, e.g., potential parameters for such an approach. Obviously, also other criteria than a match in present-day area could be taken into account (e.g., glacier length, mean ice thickness, and surface elevation at the head of the glacier). Such an approach could not only further improve the reconstruction of glaciers, but also future projections.
Hardware Requirements and Performance
For this study, we used a small cluster comprising two nodes with two 14-core CPUs each, resulting in 112 parallel threads (two threads per core). As the mass balance offset calibration converges within 12 iterations, running the calibration runs is cheap. The application of the initialization method from (Eis et al., 2019) is more expensive, as hundreds of dynamical model runs need to be performed and the dynamical runs are the most expensive computations in the OGGM. The size of the glacier and the required time stepping to ensure a numerical stability (adaptive time step following the CFL criterion) strongly influence the required computation time. The computation time needed to apply the initialization procedure to one glacier varies from glacier to glacier. In total, initializing the 517 glaciers in 1917 takes approximately 24 h on our small cluster. This implies that substantial computing resources will be needed in order to apply our method to, e.g., initialize a global model run.
The OGGM software together with the initialization method is coded in the Python language and licensed under BSD-3-Clause License. The latest version of the OGGM code is available on GitHub (https://github.com/OGGM/oggm), the documentation is hosted on ReadTheDocs (http://oggm.readthedocs.io), and the project webpage for communication and dissemination can be found at http://oggm.org. OGGM version used in this study is v1.2. The code for the initialization module is available on GitHub (https://github.com/OGGM/initialization).
Data Availability Statement
All data and code used in this study can be found on a public repository at https://github.com/OGGM/initialization, and the version used for this paper is archived with a DOI at http://doi.org/10.5281/zenodo.4625362. The OGGM version used for this study is 1.2, archived at http://doi.org/10.5281/zenodo.3597756.
JE is the main developer of the initialization module and wrote most of the paper. LL supported the development of this study and worked on methodological concepts and testing. BM and FM are the initiators of the OGGM project and helped to conceive this study. FM is the main OGGM developer.
BM and JE were supported by the German Research Foundation, grants MA 6966/1-1 and MA 6966/1-2. LL was supported by the German Research Foundation, grant FO 1269/1-1. FM acknowledges support from the Austrian Science Fund (FWF) Grant P30256-NBL.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.595755/full#supplementary-material.
Biemans, H., Siderius, C., Lutz, A., Nepal, S., Ahmad, B., Hassan, T., et al. (2019). Importance of snow and glacier meltwater for agriculture on the indo-gangetic plain. Nat. Sustainability 2, 594–601. doi:10.1038/s41893-019-0305-3
Church, J., Clark, P., Cazenave, A., Gregory, J., Jevrejeva, S., Levermann, A., et al. (2013). Sea level change, Cambridge, United Kingdom and New York, NY, USA: Cambridge University Press, 13. 1137–1216. doi:10.1017/CBO9781107415324.026
Eis, J., Maussion, F., and Marzeion, B. (2019). Initialization of a global glacier model based on present-day glacier geometry and past climate information: an ensemble approach. The Cryosphere 13, 3317–3335. doi:10.5194/tc-13-3317-2019
Farinotti, D., Huss, M., Bauder, A., Funk, M., and Truffer, M. (2009). A method to estimate the ice volume and ice-thickness distribution of alpine glaciers. J. Glaciology 55, 422–430. doi:10.3189/002214309788816759
Farinotti, D., Huss, M., Fürst, J. J., Landmann, J., Machguth, H., Maussion, F., et al. (2019). A consensus estimate for the ice thickness distribution of all glaciers on earth. Nat. Geosci. 12, 168–173. doi:10.1038/s41561-019-0300-3
Gardner, A. S., Moholdt, G., Cogley, J. G., Wouters, B., Arendt, A. A., Wahr, J., et al. (2013). A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009. Science 340, 852–857. doi:10.1126/science.1234532
Giesen, R. H., and Oerlemans, J. (2013). Climate-model induced differences in the 21st century global and regional glacier contributions to sea-level rise. Clim. Dyn. 41, 3283–3300. doi:10.1007/s00382-013-1743-7
Gregory, J. M., White, N. J., Church, J. A., Bierkens, M. F. P., Box, J. E., van den Broeke, M. R., et al. (2013). Twentieth-century global-mean sea level rise: is the whole greater than the sum of the parts?. J. Clim. 26, 4476–4499. doi:10.1175/JCLI-D-12-00319.1
Harris, I. C., Jones, P. D., Osborn, T., and Lister, D. (2014). Updated hight-resolution grids of montly climatic observations- the cru ts3.10 dataset. Int. J. Climatology 34, 623–642. doi:10.1002/joc.3771
Immerzeel, W. W., Lutz, A. F., Andrade, M., Bahl, A., Biemans, H., Bolch, T., et al. (2020). Importance and vulnerability of the world’s water towers. Nature 577, 364–369. doi:10.1038/s41586-019-1822-y
Kaser, G., Grosshauser, M., and Marzeion, B. (2010). Contribution potential of glaciers to water availability in different climate regimes. Proc. Natl. Acad. Sci. USA 107, 20223–20227. doi:10.1073/pnas.1008162107
Kienholz, C., Rich, J. L., Arendt, A. A., and Hock, R. (2014). A new method for deriving glacier centerlines applied to glaciers in Alaska and northwest Canada. The Cryosphere 8, 503–519. doi:10.5194/tc-8-503-2014
Leclercq, P. W., Oerlemans, J., Basagic, H. J., Bushueva, I., Cook, A. J., and Le Bris, R. (2014). A data set of worldwide glacier length fluctuations. The Cryosphere 8, 659–672. doi:10.5194/tc-8-659-2014
Marzeion, B., Hock, R., Anderson, B., Bliss, A., Champollion, N., Fujita, K., et al. (2020). Partitioning the uncertainty of ensemble projections of global glacier mass change. Earth’s Future 8, e2019EF001470. doi:10.1029/2019EF001470
Marzeion, B., Kaser, G., Maussion, F., and Champollion, N. (2018). Limited influence of climate change mitigation on short-term glacier mass loss. Nat. Clim. Change 8, 305–308. doi:10.1038/s41558-018-0093-1
Marzeion, B., Leclercq, P. W., Cogley, J. G., and Jarosch, A. H. (2015). Brief communication: global reconstructions of glacier mass change during the 20th century are consistent. The Cryosphere 9, 2399–2404. doi:10.5194/tc-9-2399-2015
Maussion, F., Butenko, A., Champollion, N., Dusch, M., Eis, J., Fourteau, K., et al. (2019). The open global glacier model (oggm) v1.1. Geoscientific Model. Develop. 12, 909–931. doi:10.5194/gmd-12-909-2019
Nussbaumer, S. U., and Zumbühl, H. J. (2012). The little ice age history of the glacier des bossons (mont blanc massif, France): a new high-resolution glacier length curve based on historical documents. Climatic Change 111, 301–334. doi:10.1007/s10584-011-0130-9
Pfeffer, W. T., Arendt, A. A., Bliss, A., Bolch, T., Cogley, J. G., Gardner, A. S., et al. (2014). The randolph glacier inventory: a globally complete inventory of glaciers. J. Glaciology 60, 537–551. doi:10.3189/2014JoG13J176
Radić, V., and Hock, R. (2014). Glaciers in the earth’s hydrological cycle: assessments of glacier mass and runoff changes on global and regional scales. Surv. Geophys. 35, 813–837. doi:10.1007/s10712-013-9262-y
Rastner, P., Bolch, T., Mölg, N., Machguth, H., Le Bris, R., and Paul, F. (2012). The first complete inventory of the local glaciers and ice caps on Greenland. The Cryosphere 6, 1483–1495. doi:10.5194/tc-6-1483-2012
Recinos, B., Maussion, F., Rothenpieler, T., and Marzeion, B. (2019). Impact of frontal ablation on the ice thickness estimation of marine-terminating glaciers in Alaska. The Cryosphere 13, 2657–2672. doi:10.5194/tc-13-2657-2019
Slangen, A. B. A., Adloff, F., Jevrejeva, S., Leclercq, P. W., Marzeion, B., Wada, Y., et al. (2017). A review of recent updates of sea-level projections at global and regional scales. Surv. Geophys. 38, 385–406. doi:10.1007/s10712-016-9374-2
Viviroli, D., Kummu, M., Meybeck, M., Kallio, M., and Wada, Y. (2020). Increasing dependence of lowland populations on mountain water resources. Nat. Sustainability 3, 917–928. doi:10.1038/s41893-020-0559-9
Zekollari, H., Huss, M., and Farinotti, D. (2019). Modelling the future evolution of glaciers in the european alps under the euro-cordex rcm ensemble. The Cryosphere 13, 1125–1146. doi:10.5194/tc-13-1125-2019
Zemp, M., Frey, H., Gärtner-Roer, I., Nussbaumer, S. U., Hoelzle, M., Paul, F., et al. (2015). Historically unprecedented global glacier decline in the early 21st century. J. Glaciology 61, 745–762. doi:10.3189/2015JoG15J017
Zemp, M., Huss, M., Thibert, E., Eckert, N., McNabb, R., Huber, J., et al. (2019). Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016. Nature 568, 382–386. doi:10.1038/s41586-019-1071-0
Keywords: glaciers, modeling, reconstruction, initialization, validation
Citation: Eis J, van der Laan L, Maussion F and Marzeion B (2021) Reconstruction of Past Glacier Changes with an Ice-Flow Glacier Model: Proof of Concept and Validation. Front. Earth Sci. 9:595755. doi: 10.3389/feart.2021.595755
Received: 17 August 2020; Accepted: 28 January 2021;
Published: 30 March 2021.
Edited by:Daniel Farinotti, ETH Zurich, Switzerland
Reviewed by:Harry Zekollari, Delft University of Technology, Netherlands
Jordi Bolibar, Université Grenoble Alpes, France
Copyright © 2021 Eis, van der Laan, Maussion and Marzeion. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Julia Eis, firstname.lastname@example.org