the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Initialization of a global glacier model based on presentday glacier geometry and past climate information: an ensemble approach
Fabien Maussion
Ben Marzeion
To provide estimates of past glacier mass changes over the course of the 20th century, an adequate initial state is required. However, empirical evidence about past glacier states at regional or global scales is largely incomplete, both spatially and temporally, calling for the use of automated numerical methods. This study presents a new way to initialize the Open Global Glacier Model from past climate information and presentday glacier states. We use synthetic experiments to show that even with these perfectly known but incomplete boundary conditions, the problem of model initialization is an illposed inverse problem leading to nonunique solutions, and we propose an ensemble approach as a way forward. The method works as follows: we generate a large set of physically plausible glacier candidates for a given year in the past (e.g., 1850 in the Alps), all of which are then modeled forward to the date of the observed glacier outline and evaluated by comparing the results of the forward runs to the presentday states. We test the approach on 2660 Alpine glaciers and determine error estimates of the method from the synthetic experiments. The results show that the solution is often nonunique, as many of the reconstructed initial states converge towards the observed state in the year of observation. We find that the median state of the best 5 % of all acceptable states is a reasonable best estimate. The accuracy of the method depends on the type of the considered observation for the evaluation (glacier length, area, or geometry). Trying to find past states from only presentday length instead of the full geometry leads to a sharp increase in uncertainty. Our study thus also provides quantitative information on how well the reconstructed initial glacier states are constrained through the limited information available to us. We analyze which glacier characteristics influence the reconstructability of a glacier, and we discuss ways to develop the method further for realworld applications.
 Article
(4802 KB) 
Supplement
(6790 KB)  BibTeX
 EndNote
Glaciers contributed significantly to past sealevel rise (SLR; e.g., Gregory et al., 2013; Slangen et al., 2017a; WCRP Global Sea Level Budget Group, 2018; Wouters et al., 2019; Zemp et al., 2019), and they will continue to be a major contributor in the 21st century (e.g., Church et al., 2013; Slangen et al., 2017b; Hock et al., 2019). A large fraction of this contribution will be caused by the ongoing adjustment of glaciers to previous climate change (Marzeion et al., 2014, 2018). Reconstructions of past glacier mass change are therefore not only necessary to determine the budget of past sealevel change (Gregory et al., 2013) and to increase the confidence in projections (by allowing the quantification of the agreement with observations; Marzeion et al., 2015); they also enable us to quantify the pattern of the ongoing adjustment of glaciers to presentday climate. Estimates of global glacier mass change are based on in situ measurements in mass and length changes (e.g., Zemp et al., 2015; Leclercq et al., 2011), on remote sensing techniques (e.g., Gardner et al., 2013; Jacob et al., 2012; Bamber et al., 2018; Wouters et al., 2019), or on mass balance modeling driven by climate observations (Marzeion et al., 2012, 2015). Since observations of temperature, and to a smaller degree precipitation, are more ubiquitous (e.g., Harris et al., 2014) than glacier observations (WGMS, 2018), reconstructions of glacier change produced by forcing a glacier model with climate observations have the potential to increase the understanding of past glacier behavior. Finally, reconstructing glacier change based on climate model output allows us to test the skill of climate models (Goosse et al., 2018).
A number of global glacier models were developed in the past (e.g., Radić and Hock, 2011, 2014; Giesen and Oerlemans, 2012, 2013; Marzeion et al., 2012, 2014; Huss and Hock, 2015; Maussion et al., 2019a). The more recent and complex of these models (e.g., Huss and Hock, 2015; Maussion et al., 2019a) require digital elevation models (DEMs) and outlines from the Randolph Glacier Inventory (RGI; Pfeffer et al., 2014) to derive the initial surface hypsometry. Hence, their starting date of a glacier evolution simulation depends on the recording date of the DEM and the outline, which typically do not coincide with one another, nor with the required starting date of a projection. The model of Huss and Hock (2015) indicates a high sensitivity to the initial ice volume. Similarly, Maussion et al. (2019a) remark that great uncertainties, especially on local and regional scales, derive from unknown initial conditions.
Despite the importance of glacier contribution to past sealevel rise, so far only one model was able to provide estimates of glacier mass changes over the course of the entire 20th century on the global scale (Marzeion et al., 2012). All other global modeling studies limit their application to the recent past and future projections. The reconstruction by Marzeion et al. (2012) was possible because of the highly parameterized representation of ice dynamics and glacier geometry change, applying a volume–area–time scaling to translate mass change into surface area and elevation range changes. Based on this approach, it was sufficient to iteratively optimize one variable (glacier size in the year of interest, e.g., 1850) such that when run forward to the year of the observed glacier outline, the modeled glacier area agreed with the observed glacier area.
An increase in model complexity impedes the process as more and more variables are required for initialization. Flow line models require input data along the coordinates of the flow line (e.g., bed topography, surface elevations, and/or widths), and thus more complex initialization methods are needed. For example, van Pelt et al. (2013) developed an iterative inverse method to reconstruct distributed bedrock topography and simultaneously initialize an ice flow model. Zekollari et al. (2019) added an ice flow model to Huss and Hock (2015), which required an automated initialization for glaciers in 1990 (prior to the glacier inventory date) to avoid spinup issues and so that the reconstructed initial states fit the glacier geometry at the inventory date after being modeled forward. By choosing a decadelong initialization, they avoid problems of nonuniqueness (as we discuss below), but they raise the question of how arbitrarily this date can be chosen. Similar approaches exist for the initialization of ice sheet models, where most work focuses on estimating the presentday state of ice sheets in order to make accurate projections of future ice sheet change (e.g., Heimbach and Bugnion, 2009; Lee et al., 2015; Mosbeux et al., 2016). Goelzer et al. (2018) divide the existing initialization approaches into three methods: spinup, assimilation of velocity, and assimilation of surface elevation. Spinup procedures are typically used for longterm and paleoclimate simulations; the required spinup time is unknown and can be relative long. Additionally, the reconstruction cannot be expected to represent effects from internal climate variability correctly. The data assimilation approaches typically determine model parameters (e.g., basal parameters like basal friction or bedrock topography) that reduce the mismatch between observed and modeled velocities or surface elevations.
In this study, we aim to identify fundamental limitations that narrow the reconstruction of past glacier states from presentday geometries, under the assumption of perfectly known boundary conditions and a perfect glacier model. Specific research questions are as follows.

To which degree does the past evolution of a glacier constrain its presentday geometry?

How much information does the presentday glacier geometry contain about its past states?

Is it possible to reconstruct past glacier states from the partial information available to us?

How far can we go back in time to have an initial geometry that still determines the presentday glacier geometry?

Which glacier attributes influence the answers to the questions above?
To this aim, we present a new method estimating past glacier states and apply it to synthetic numerical experiments, and we show the obstacles that need to be overcome before applying our method to realworld problems. After introducing the relevant features of the Open Global Glacier Model (OGGM; Maussion et al., 2019a) in Sect. 2.1, we describe the design of the synthetic experiments in Sect. 2.3. The synthetic framework serves to test the skill of our approach in a surrogate model world where everything is known, and it allows us to apply data denial experiments to address the questions listed above. The initialization method is presented in Sect. 2.4. The developed method consists of three steps: generation of plausible glacier states, identification of glacier state candidates, and their evaluation based on the misfit between the modeled and the observed geometry at the year of the observation. We applied our approach to 2660 Alpine glaciers and present the results for the reconstructed initial states in the year 1850 in Sect. 4.1. The influence of the considered type of observation (e.g., glacier length, area, or geometry) is shown in Sect. 4.2 and a statistical analysis of glacier attributes that influence the reconstructability of a glacier is presented in Sect. 4.3. Finally, we summarize the results and discuss the limitations of the method and its applicability to realcase studies, as well as needed and possible future developments in Sect. 6.
2.1 The Open Global Glacier Model
The Open Global Glacier Model (OGGM; Maussion et al., 2019a) is an opensource numerical framework that allows the modeling of past and future changes of any glacier in the world. Starting with a glacier outline, provided by the Randolph Glacier Inventory (RGIv6.0; Pfeffer et al., 2014), a suitable surface DEM is automatically downloaded and interpolated to a local grid. The size of the local grid is given by a border parameter, which is the number of grid points outside the glacier boundaries. We choose a border value of 200 grid points to ensure that large glacier states can also be generated. The resolution of the map topography dx depends on the size of the glacier ($\mathrm{d}x=a\sqrt{S}$, with a=14 m km^{−1} and S the area of the glacier in square kilometers) and is constrained to $\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\le \mathrm{d}x\le \mathrm{200}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$. After the preprocessing, glacier centerlines are computed using a geometrical routing algorithm (adapted from Kienholz et al., 2014). They are then considered to be glacier flow lines, and grid points are generated using a fixed, equidistant grid spacing, which is twice that of the underlying 2D map topography. Surface elevations along the flow line coordinates are then obtained from the underlying topography file, and glacier section widths are computed by intersecting the normal of the flow line to the boundaries of the glacier. By making assumptions about the shape of the bed (parabolic, rectangular, or a mix of both), OGGM estimates the ice thickness with a massconservation approach (Farinotti et al., 2009, 2017, 2019). Information on bed topography at each grid point results from the calculated ice thickness and the surface elevation. From this, the glacier length, area, and volume can be determined. These values depend strongly on the surface topography and are based on the (often wrong) assumption that the recording date of DEM and that of the outline coincide. The dynamical flow line model of OGGM can then be used to determine the evolution of the glacier under any given climate forcing by solving the shallow ice approximation along the flow lines.
The mass balance is computed at each grid point using climate data (monthly temperature and precipitation). Climate data can be used from different sources, including gridded observations or reanalyses for past climate, projections for future climate, or randomized climate time series. The purpose of forcing the mass balance model with randomized climate is to easily produce a great number of realistic climate forcings representative of a given time period, characterized by a center year y_{0} and a window size h (typically 31 years). All climate years $\in [{y}_{\mathrm{0}}\frac{h\mathrm{1}}{\mathrm{2}},{y}_{\mathrm{0}}+\frac{h\mathrm{1}}{\mathrm{2}}]$ are then shuffled infinitely in the next step. Additionally, it is possible to set a temperature bias β, which shifts all values of the temperature series towards warmer or colder climates.
Identically to the study of Maussion et al. (2019a), we only calibrate the mass balance model, while the creep parameter A and the sliding parameter f_{s} are the same for each glacier and set to their default values ($A=\mathrm{2.4}\times {\mathrm{10}}^{\mathrm{24}}$ s^{−1} Pa^{−3}, f_{s}=0, no lateral drag). The following massbalancerelated parameter values were used in this study: p_{f}=1.75, ${T}_{\mathrm{Melt}}=\mathrm{1.5}$ ^{∘}C, T_{Liquid}=2.0 ^{∘}C, and $\mathrm{\Gamma}=\mathrm{6.5}$ K km^{−1}. This parameter set was determined with a crossvalidation performed with the HISTALP data set and tested for the 41 Alpine glaciers with more than 5 years of mass balance observation. For more details concerning the glacier model (e.g., the mass balance calibration or sensitivities to the dynamical parameters of the model), please refer to Maussion et al. (2019a) and http://docs.oggm.org (last access: 10 December 2019).
2.2 Problem description
Here, we define a glacier state (hereinafter referred to as state) as follows.
 Definition 1.

Let m∈ℕ be the total number of grid points of all flow lines of a glacier. Then ${s}_{t}=({z}_{t},{w}_{t},b)$ is a glacier state at time t, with surface elevation ${z}_{t}\in {\mathbb{R}}_{+}^{m}$, widths ${w}_{t}\in {\mathbb{R}}_{+}^{m}$, and bed topography $b\in {\mathbb{R}}_{+}^{m}$. The set ${\mathcal{S}}_{{t}_{i}}=\mathit{\left\{}{s}_{t}\rightt={t}_{i}\mathit{\}}$ contains all physically plausible glacier states at time t_{i}.
The construction of an initial state is an inverse problem and can be defined in opposition to the direct problem. The direct problem corresponds to a forward model run: given an initial state ${s}_{{t}_{\mathrm{0}}}\in {\mathcal{S}}_{{t}_{\mathrm{0}}}$ at time t_{0}, the state s_{t}∈𝒮_{t} at time t>t_{0} can be computed by
where ${G}_{\mathrm{past}}:{\mathcal{S}}_{{t}_{\mathrm{0}}}\to {\mathcal{S}}_{t}$ is an operator representing the equations of OGGM, using known climate time series as the boundary condition.
For inverse problems, the solution is known by direct observations: ${s}_{{t}_{e}}={s}_{{t}_{e}}^{\mathrm{obs}}$, whereas the desired initial state ${s}_{{t}_{\mathrm{0}}}$ is unknown. The inverse problem consists of finding the initial state ${s}_{{t}_{\mathrm{0}}}\in {\mathcal{S}}_{{t}_{\mathrm{0}}}$, such that the forward modeled solution at time t_{e} fits the observations from the same year t_{e}:
Unfortunately, we do not have an explicit formulation for ${G}_{\mathrm{past}}^{\mathrm{1}}$ in our case. A backwards reconstruction is impeded by the nonlinear interaction between glacier geometry, ice flow, and mass balance. Optimization methods can be used to solve inverse problems. To this end, we introduce a minimization problem such that the forward modeled state is as close as possible to the observation:
with
This function calculates the averaged difference in surface elevation and width between the observed and forward modeled glacier states. Differences in bed topography can be neglected, as we assume the bed topography to remain the same over the inspected time period.
In many cases, however, OGGM forward integrations of different initial states result in very similar states at time t_{e}. This implies that there exist many local minima of the function $j\left({s}_{{t}_{\mathrm{0}}}\right)$. As uncertainties of the model can safely be assumed to be larger than the differences between those states at time t_{e}, it is impossible to identify the global minimum of $j\left({s}_{{t}_{\mathrm{0}}}\right)$. That is, the solution of our inverse problem is nonunique.
The objective of our approach is therefore to identify the set ${\mathcal{S}}_{{t}_{\mathrm{0}}}^{\mathit{\u03f5}}$ of all states, which correspond to the observed state ${s}_{{t}_{e}}^{\mathrm{obs}}$ within a given uncertainty ϵ after being modeled forward. We call this condition the acceptance criterion:
The function $J\left({s}_{{t}_{\mathrm{0}}}\right)$ is called the fitness function in the following. Assuming a vertical error of 5 m in x and a horizontal error of 10 m in w, we propose setting $\mathit{\u03f5}=(\mathrm{5}\phantom{\rule{0.125em}{0ex}}\mathrm{m}{)}^{\mathrm{2}}+(\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m}{)}^{\mathrm{2}}=\mathrm{125}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}$. These numbers can be changed easily, and in a realworld application they should be based on the vertical uncertainty of the reconstructed ice thickness and the horizontal uncertainty of the used outline. All states ${s}_{{t}_{\mathrm{0}}}\in {\mathcal{S}}_{{t}_{\mathrm{0}}}^{\mathrm{125}}$ that have a fitness value smaller than 1 are called acceptable states. The first expectation would be that the glacier candidate with the smallest fitness value is also the best solution. However, due to uncertainties that derive from the model integration itself, this is not always the case. As an alternative, we determine the 5th quantile of all states in ${\mathcal{S}}_{{t}_{\mathrm{0}}}^{\mathit{\u03f5}}$. This set contains the best solutions of all acceptable states referring to their fitness values. We choose the median state as a representative of this set and compare the state with the minimal fitness value and the median state in Sect. 4.1.
2.3 Synthetic experiments
2.3.1 Design
We create a time series of glacier states, which range from the target year of initialization t_{0} (e.g., 1850) to the present date t_{e} (see Fig. 1 for an example). These are the glacier states that we aim to reconstruct with our initialization method, using only partial information (here the observed state at present day). This type of experiment is sometimes called “inverse crime” in the inverse problem literature (e.g., Colton and Kress, 1992; Henderson and Subbarao, 2017), and we explain this rationale below. To generate the synthetic experiments, we apply a random climate scenario (window size h=31 years and center year ${y}_{\mathrm{0}}={t}_{\mathrm{0}}=\mathrm{1850}$) and run the model 600 years forward (see Fig. 1c). The temperature bias is set to $\mathit{\beta}=\mathrm{1}$ K to ensure that a relatively large 1850 glacier state is created (as expected for most real glaciers at the end of the Little Ice Age). The resulting state is defined to be the synthetic experiment state in year t_{0} (see Fig. 1a). We model this state forward, applying the past climate time series from t_{0} until t_{e} (here 2000) (see Fig. 1d) and obtain the observed state of the synthetic experiment (see Fig. 1b). Thanks to the initial temperature bias of $\mathit{\beta}=\mathrm{1}$ K, these synthetic states in t_{e} are very close to the real observed states in 2000 on average (total area difference for the Alps of about 1 %, but individual glaciers can vary), and the total synthetic glacierized area in 1850 fits well to an estimate of Zemp et al. (2006) (see Appendix A for more details). We call the states derived from the synthetic experiment ${s}_{t}^{\mathrm{exp}}$.
2.3.2 Rationale
These synthetic states therefore provide a realistic setting with a strong advantage over actual observations: they are perfectly known. As stated in the introduction, reconstructing past glacier states is a complex inverse problem, the accuracy of which will depend on (i) the uncertainties in the boundary conditions (climate, glacier bed, etc.), (ii) the uncertainties in the glacier model itself, and (iii) a theoretical lower bound (termed “reconstructability” in this study) tied to the characteristics of the glacier itself (slope, size, the past climate, etc.). The main objective of the synthetic experiments is to separate these issues from one another and to focus on point (iii) only. This allows us to isolate and understand the limitations and errors of the developed method itself, as opposed to uncertainties that derive from unknown boundary conditions and model parameters. They also allow us to realize data denial experiments and detect which kinds of observations are necessary to reduce the uncertainties of our reconstruction (Sect. 4.2) and to determine which glacier characteristics affect the reconstructability of a glacier (Sect. 4.3).
2.4 Reconstruction of initial glacier states
Our initialization method consists of three main steps: generation of a set of physically plausible glacier states ${\mathcal{S}}_{{t}_{\mathrm{0}}}$, identification of glacier candidates ${s}_{{t}_{\mathrm{0}}}\in {\mathcal{S}}_{{t}_{\mathrm{0}}}$, and their evaluation based on the fitness function $J\left({s}_{{t}_{\mathrm{0}}}\right)$ (see Sect. 2.2).
2.4.1 Generation of potential glacier candidates
In a first step, we generate a set of different, physically plausible states from which we will pool our candidates (Fig. 2a). For this purpose, we utilize a random mass balance model with a window size of h=31 years and the center year y_{0}=t_{0} to create different climate conditions. Obviously, we do not use the same permutation as for the creation of the synthetic experiments (Sect. 2.3). This procedure generates a climate representative of a given time period with an interannual variability uncorrelated to that of the original period. For each random climate a different way of permutation is used. This ensures that all generated climate time series differ from each other, but at the same time all represent the climate conditions around t_{0} (and an associated temperature bias β). The infinite permutation is necessary to obtain a time series that is long enough for the glaciers to reach an equilibrium (while maintaining the impact of interannual climate variability) with the forcing climate (here 600 years). To create a large set of states, we additionally vary the temperature bias β. Glaciers respond differently to changes in climate, and thus the required temperature biases vary from glacier to glacier and have to be inferred. We start with temperature biases $\mathit{\beta}\in [\mathrm{2},\mathrm{2}]$ K. If β=2 K is not large enough to result in a presentday glacier with zero ice thickness, higher values will be used. If $\mathit{\beta}=\mathrm{2}$ K is not small enough to result in a glacier that reaches the boundary of the local grid (200 grid points outside of the glacier outline), smaller values will be used.
2.4.2 Selection of candidates
Figure 2a shows the evolution of the volume of the generated glacier states over time. In the first years, the time series clearly diverge (mostly caused by the temperature bias β), but after a certain time all time series begin to fluctuate around an equilibrium value. We refer to the period of fluctuations around the assumed equilibrium as the stagnation period. During the stagnation period the glacier volume does not increase or decrease strongly in comparison to the total volume change since the beginning of the simulation. We define t_{stag} as the point in time when all trajectories have reached this stagnation period and choose the upper 10 volume trajectories, corresponding to the lowest temperature biases, to determine t_{stag}. To this end, we smooth each of the 10 curves with a 10year rolling window and calculate the time point of their first maximum. t_{stag} is defined as the latest of all previously determined time points (see Fig. 2b).
Defining t_{stag} is necessary because we determine initial glacier states at t_{0}=1850 and the searched glaciers are assumed to be in equilibrium with the climate around 1850. Hence, each state that fluctuates around an equilibrium value is a potential glacier state candidate (in the following referred to as “candidate”). This holds true for all states s_{t} with t>t_{stag}. Depending on t_{stag} and the number of successfully completed random climate runs n_{r} (number of grey lines in Fig. 2), the sample size is n_{r}(600−t_{stag}) (glacier states are stored yearly). The sample size is sufficiently large for all cases, e.g., in the case of the Guslarferner (Fig. 2) the sample contains approximately 44 500 members. In order to avoid testing very similar states, we classify all states by their volume and select one candidate per class. We choose n equidistantly and approximately uniform distributed classes, where n (default n=200) is the number of candidates to evaluate in step three.
2.4.3 Evaluation
The last step evaluates all previously selected candidates. Each candidate is used as the initial condition for a forward run, using observed past climate time series, e.g., from t_{0}=1850 until t_{e}=2000. All runs use the same model parameter set, except for the initial condition and exactly the same climate time series (e.g., no temperature bias β is applied for the past climate runs). Afterwards, we compare the resulting modeled states ${s}_{{t}_{e}}$ with an observed state ${s}_{{t}_{e}}^{\mathrm{obs}}$ (here taken from the synthetic experiment) by applying the fitness function $J\left({s}_{{t}_{\mathrm{0}}}\right)$ (Eq. 5). This function calculates the averaged difference between the glacier geometries at the grid points, more specifically between the surface elevations ${z}_{{t}_{e}}$ and the widths ${w}_{{t}_{e}}$, of the observed and the modeled glacier. In Fig. 2c the candidates are colored by their fitness value.
We tested our approach on Alpine glaciers. The glacier outlines are taken from the Randolph Glacier Inventory (RGI v6.0, region 11; Pfeffer et al., 2014). We use topographical data from the Shuttle Radar Topography Mission (SRTM) 90 m Digital Elevation Database v4.1 (Jarvis et al., 2008). The SRTM acquisition date (2000) matches that of the RGI well (2003 for most glaciers). The climate data set we use for this approach is the HISTALP database (Auer et al., 2007, http://www.zamg.ac.at/histalp, last access: 10 December 2019). The temperature time series covers the period 1780 to 2014 and the precipitation time series 1801 to 2014. Both data sets are available on a regular grid of 5 min resolution (approximately 9.3 km in the Alps).
We generate synthetic experiments (see Sect. 2.3) for all glaciers in the Alps, and we determine their glacier states in 1850 if the area of the observed synthetic state ${s}_{\mathrm{2000}}^{\mathrm{exp}}$ is larger than 0.01 km^{2}. This value is consistent with the minimumarea threshold of the RGI. The condition is satisfied for 2660 synthetic experiments of the 3927 glaciers included in the Randolph Glacier Inventory in central Europe (region 11).
Here we show the results for two example glaciers in 1850 as well as an error analysis for all tested glaciers (Sect. 4.1), the influence of the choice of the fitness function on the quality of our results (Sect. 4.2), and a statistical analysis of glacier attributes (including glacier response time) that influence the reconstructability of a glacier (Sect. 4.3).
4.1 Initial glacier states in 1850
Following the method described in Sect. 2.4, we determine reconstructed initial glacier states in t_{0}=1850. Figures 3 and 4 show the results of Guslarferner, as an example glacier with a large set of accepted candidates. A second case with a more narrow set of acceptable states, Hintereisferner, is shown in Figs. 5 and 6. More examples can be found in the Supplement.
In particular the result of Guslarferner shows clearly that the determination of past states is not unique (see Fig. 3). Multiple initial states (violet and blue) merge to the observed state in the year of observation. The fitness values, which means the averaged difference between the forward modeled states and the observation at t_{e}=2000, are extremely small for most candidates. The fitness values of all candidates range from 1.08 $\times {\mathrm{10}}^{\mathrm{6}}$ to 7.98. Only 16 of the 200 candidates have a fitness value higher than 1 and thus do not fulfill the acceptance criterion (Eq. 5); for these states, the glacier in 1850 is too small to reach the volume of the observed glacier within 150 years. Also, Fig. 4 illustrates the diversity of the different acceptable solutions (grey, shadowed area). The length of all states in ${S}_{\mathrm{1850}}^{\mathrm{125}}$ varies between 0.98 and 8.1 km. The acceptance criterion in this example is not strong enough to provide any information about the searched state in t_{0}=1850, as any of the candidates would lead to an acceptable result. Figure 4 also shows the 5th percentile of all acceptable states (blue shadowed area). This set contains the 5 % best solutions, based on the fitness value. All candidates of the 5th percentile are in close proximity to the synthetic experiment. The range of fitness values of all candidates of the 5th percentile is $[\mathrm{1.08}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{7.95}\times {\mathrm{10}}^{\mathrm{5}}]$, and the length of the states in 1850 only varies from about 3.6 to 5.3 km. All these candidates match the synthetic experiment in t_{e}=2000 very well and converge to the synthetic experiment by 1900 at the latest, which can be seen in Fig. 4c. As a representative of this set, we choose the median state of the 5th percentile of ${S}_{t}^{\mathrm{125}}$ (in the following referred to as ${s}_{t}^{\mathrm{med}}$). Figure 4a shows that the surface elevation of ${s}_{\mathrm{1850}}^{\mathrm{med}}$ in 1850 corresponds very well to the synthetic experiment, whereas the state with the minimum fitness value (in the following referred to as ${s}_{t}^{\mathrm{min}}$) mismatches the synthetic experiment at the tongue of the glacier. Regarding the volumes, ${s}_{t}^{\mathrm{med}}$ exactly matches the volume of the synthetic experiment in 1850, whereas the volume of ${s}_{t}^{\mathrm{min}}$ differs by 0.4 km^{3}.
In the case of Hintereisferner (see Fig. 5) the fitness values of most candidates are large compared to those of Guslarferner. Only a few candidates have extremely small fitness values and the past state is thus much more narrowly confined. The different states need more time to adapt to the climate conditions and therefore they do not converge as quickly to one state. As a result, the differences between the forward modeled states and the observed state in 2000 are larger. The fitness values of all candidates range between $\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{5}}$ and 43.
A total of 36 candidates fulfill the acceptance criterion (Eq. 5). Figure 6 shows that the acceptance criterion in this case confines the result better than in the case of Guslarferner. The length of all glaciers in ${S}_{\mathrm{1850}}^{\mathrm{125}}$ ranges from 8.4 to 12.3 km. In this case the 5th quantile of ${S}_{t}^{\mathrm{125}}$ is again in close proximity to the synthetic experiment, and all candidates of the 5th quantile have extremely small fitness values (between $\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{5}}$ and $\mathrm{5.4}\times {\mathrm{10}}^{\mathrm{4}}$). The length of the candidates of the 5th quantile in 1850 only varies from 9.1 to 10.3 km and is thus more precise than in the Guslarferner example. In this example, all candidates of the 5th percentile converge no later than 1920 to the state of the synthetic experiment. Here ${s}_{t}^{\mathrm{min}}$ matches the surface elevation of the synthetic experiment in 1850, as well as the volume trajectory over time, slightly better than ${s}_{t}^{\mathrm{med}}$, but the volume differences to the synthetic experiments in 1850 are also very small in this example (0.004 km^{3} for ${s}_{t}^{\mathrm{min}}$ and −0.08 km^{3} for ${s}_{t}^{\mathrm{med}}$).
For both examples we were able to show that our method is able to recover the state in t_{0}=1850 of the synthetic experiment by only using information about the observed state of the synthetic experiment in t_{e}=2000 and combining it with information about the past climate evolution. ${s}_{t}^{\mathrm{med}}$ as well as ${s}_{t}^{\mathrm{min}}$ match the synthetic experiment in t_{0}=1850 extremely well. In the following, we provide an error analysis including all glaciers in the Alps to which we applied our method. For each of the 2660 glaciers we have calculated the absolute volume error to the synthetic experiment:
where v^{exp}(t) is the volume of the synthetic experiment in year t, and v^{med∕min}(t) is the volume of ${s}_{t}^{\mathrm{med}}$ or ${s}_{t}^{\mathrm{min}}$ in the same year t. Figure 7a shows the absolute volume errors in cubic kilometers for ${s}_{t}^{\mathrm{med}}$, as well as for ${s}_{t}^{\mathrm{min}}$.
Whereas the absolute volume errors in 1850 vary widely from approximately −1.1 to 2.9 km^{3}, they reduce rapidly within 60 years. In 1910, the errors range from approximately −0.25 to 0.17 km^{3}. The range of errors in the first 60 years is largely influenced by a few single outliers. Differences between ${s}_{t}^{\mathrm{min}}$ and ${s}_{t}^{\mathrm{med}}$ are small. Figure 7b shows the median and the range of the 5th–95th percentiles of ${e}_{\mathrm{abs}}^{\mathrm{med}}$ and ${e}_{\mathrm{abs}}^{\mathrm{min}}$ over time, indicating the robustness of our method. The median of e_{abs} of both analyzed states is very small; 0.00028 km^{3} for ${s}_{t}^{\mathrm{min}}$ and 0.00076 km^{3} for ${s}_{t}^{\mathrm{med}}$ in 1850. The improvement with time can also be seen here: the median of ${e}_{\mathrm{abs}}^{\mathrm{med}}\left(\mathrm{1910}\right)$ is of the order of 10^{−5} km^{3}, and that of ${e}_{\mathrm{abs}}^{\mathrm{min}}\left(\mathrm{1910}\right)$ is of the order of 10^{−6} km^{3} .
As our test site contains large and small sized glaciers, we also evaluate relative errors (%):
Figure 8a shows the histogram of the relative errors in 1850, whereas the evolution from 1850 to 2000 of the median and the 5th–95th percentiles of the relative errors is shown in Fig. 8b.
The median of the relative volume errors in 1850 is −0.97 % for ${s}_{t}^{\mathrm{med}}$ and −2.69 % for ${s}_{t}^{\mathrm{min}}$. The 95th percentile value of ${e}_{\mathrm{rel}}^{\mathrm{med}}$ is 70 %. With 48 % the value of ${e}_{\mathrm{rel}}^{\mathrm{min}}$ is smaller for ${s}_{t}^{\mathrm{min}}$. Whereas ${s}_{t}^{\mathrm{min}}$ has a slightly smaller 5th–95th percentile range than ${s}_{t}^{\mathrm{med}}$ in 1850, the median error of ${s}_{t}^{\mathrm{med}}$ is slightly smaller than that of ${s}_{t}^{\mathrm{min}}$. Both states fit the synthetic experiment well. In many cases, ${s}_{t}^{\mathrm{med}}$ is equal to ${s}_{t}^{\mathrm{min}}$, but for some glaciers either ${s}_{t}^{\mathrm{min}}$ or ${s}_{t}^{\mathrm{med}}$ has a clearly better performance. In all cases, the uncertainties quickly decrease after around 1900 to 1930.
Figure 8a also shows that the error distribution is skewed, and our method has a slight tendency to underestimate the glacier volume. Although 64 % of the relative errors have a negative sign, a few large positive outliers influence the mean error and shift it to a positive value of 16 % (in 1850) for the minimum states or 23 % (in 1850) in the case of the median states.
4.2 Impact of the fitness function
For the evaluation of the glacier candidates we used a fitness function based on differences in the geometry of the glacier (see Eq. 5). In this section we want to test the influence of limited information on glacier geometry on the reconstructability of past glacier states. Thus, we additionally evaluate the candidates by only using information about the glacier area or length.
For the glacierareabased evaluation, we used the following fitness function:
where ${A}_{{t}_{e}}$ is the glacier area at time t_{e}. The fitness function that takes only information about the glacier length l(t_{e}) at time t_{e} into account is similar:
For each glacier at our test site, we also evaluate the 200 candidates with the fitness functions J_{A} and J_{l}. For each evaluation method (geometry, area and length based), we determine the state with the minimal fitness function^{1}and calculate the relative volume error between it and the synthetic experiment.
Figure 9 shows the relative errors of all three evaluation methods. Figure 8a shows the distribution of the relative errors of the 5th–95th percentiles in 1850, whereas the evolution from 1850 to 2000 of the median and the 5th–95th percentiles of the relative errors is shown in Fig. 8b. The more information taken into account for the evaluation, the smaller the errors. The greatest uncertainties are associated with using the glacierlengthbased fitness function (Eq. 9), whereas the differences between the areabased evaluation (Eq. 8) and the geometrybased evaluation (Eq. 5) are small. While the median errors in 1850 of the geometry and the areabased evaluations are close (−2.69 % for the geometry and −2.83 % for the area approach), the median error in 1850 of the glacier length evaluation has, with 107 %, the worst performance. This also applies for the values of the 95th percentile; 95 % of the tested cases have a relative volume error smaller than 1043 % in 1850 if the lengthbased fitness function is used for the evaluation. In contrast to the other two evaluations, this approach overestimates the volume. For the areabased evaluation, 95 % of the tested glaciers have an error smaller than 90 %, and for the geometrybased fitness function the error is smaller than 49 %. This shows that the advantage of using the geometry instead of the glacier area to evaluate the candidates is not very high; both evaluations show a very good performance. Especially if the states are modeled forward (e.g., to 1900), both approaches perform well. However, it is not advisable to use the glacierlengthbased evaluation.
4.3 Reconstructability
The examples from Sect. 4.1 as well as in the Supplement indicate a high variation in the number of viable reconstructed candidates between glaciers. This number can range from a few viable solutions in a welldefined range to many solutions without any constraints (all tested candidates have the same fitness value). In other words, some glaciers can be reconstructed easily, and some cannot.
We define a new measure of reconstructability r, where we set the volume range of the 5th percentile in relation to the volume range of all acceptable states of the glacier:
For a glacier with a unique solution, this measure is equal to 1. If all accepted candidates have exactly the same fitness value, the measure will be zero (this occurs if all candidates converge to exactly the same state before the year 2000). Thus, a small measure represents a glacier with low reconstructability, and a measure close to 1 implies a higher reconstructability of the glacier. For example, r is equal to 0.857 for Hintereisferner and 0.879 for Guslarferner. The similarity of the two values can be explained by the similar proportion of the range of the 5th percentile to the range of all acceptable states in both cases (see Figs. 3 and 6). A histogram of the reconstructability values of all 2660 tested glaciers in the Alps is shown in Fig. 10a. The distribution is bimodal and slightly skewed towards a high reconstructability. Values in the middle range are rare.
What glacier characteristics will influence this reconstructability? The working hypothesis is that it is likely to be associated with the concept of glacier “response time” (here formulated qualitatively). Glaciers with a short response time tend to be less sensitive to initial conditions and will “forget” their initial state after a short period of time. This will probably lead to low reconstructability values. Inversely, glaciers with a long response time will be easier to reconstruct.
To test this hypothesis, we used the efolding approach (as defined in Oerlemans, 1997, 2001) and calculated the time response to a step function. To this end, we first run the 1850 state of the synthetic experiment glacier into an equilibrium state by using a constant climate (mean climate of the years 1835–1865, temperature bias $=\mathrm{1}$ K). We choose the same settings that were used for the generation of the synthetic experiments in order to obtain an equilibrium state s_{eq1} close to our synthetic experiment in 1850. Next, we apply to s_{eq1} a constant climate obtained by the mean climate of the years 1850–1880 using no temperature bias and receive the corresponding equilibrium state s_{eq2} (i.e., a step change of 1 K). We calculate the efolding time of these two states for each glacier, but we exclude the glaciers where the volume of s_{eq2} reaches zero (which was the case for approximately 500 glaciers out of 2660^{2}).
The scatter plot in Fig. 10b indicates a relation between the reconstructability measure and response time. The variance of the response time increases for reconstructability values close to 1. Dependencies with the reconstructability could also be detected for the position of the equilibrium line altitude (ELA) (Fig. 10c), the mean surface slope in 2000 (Fig. 10d), and the mean surface slope in 2000 of the last third of the glacier (Fig. 10e).
Furthermore, we calculated correlations of both reconstructability and response time with the following variables: glacier length (in 2000), area (in 2000), volume (in 2000), equilibrium line altitude (ELA) in 2000, equilibrium line altitude change from 1850 to 2000, mean surface slope (in 2000), and mean surface slope over the lowest third of the glacier (in 2000) (Fig. 11).
The variable explaining reconstructability best is the glacier response time (correlation of 0.54). Both values correlate with the same glacier characteristics. Contrary to a common misunderstanding, glacier length, area, and volume do not correlate well with the reconstructability measure or with the response time. The variable having the main influence is slope: generally, the larger the slope, the lower the reconstructability measure or the response time of a glacier. These findings coincide with results from Lüthi (2009), Zekollari and Huybrechts (2015), and Bach et al. (2018), who concluded that response times depend more on the steepness of the surface than on the glacier size. The correlation of the mean surface slope can be further increased by taking the lowest third of the glacier. In addition, the position of the ELA in 2000 also influences the reconstructability, whereas the ELA change from 1850 to 2000 only plays a minor role.
Taken alone, these correlation values remain quite low and do not provide enough predictive power to create a statistical model of reconstructability. However, they provide a good indication about which factors should be taken into account for future applications.
For this study we used a small cluster comprising two nodes with two 14core CPUs each, resulting in 112 parallel threads (two threads per core). Our method requires running hundreds of dynamical model runs for each single glacier, and, as described in Maussion et al. (2019a), the dynamical runs are the most expensive computations. The size of the glacier and the required time stepping to ensure a numerical stability strongly influence the required computation time. The computation time needed to apply our initialization procedure to one glacier varies from 30 s to 26 min. In total, initializing the 2660 glaciers in 1850 takes about 3.75 d on our small cluster.
In this study, a new method to initialize past glacier states is presented and applied to synthetic experiments. Assuming a perfectly known world allows us to identify the errors of our method alone and to separate them from uncertainties in observations and errors introduced by model approximations, a task impossible to realize in realworld applications. However, the synthetic experiments do not allow external validation, e.g against past outlines derived from moraines, historical maps, or remote sensing (such as provided by GLIMS; Raup et al., 2007). Model uncertainties will have to be accounted for and will have to be compared to and added to the theoretical lower bound discussed in this study. Similarly, our results do not provide information about actual past glacier mass change. Since in our synthetic experiments glaciers states in 2000 may be different from the real ones, the modeled initial glacier states in 1850 do not correspond to reality either. The past states determined in this study can only serve to verify the functionality of the developed method.
Our results have shown that the solutions are not unique. Multiple candidates match the observation in t_{e}=2000, sometimes with a large spread. This raises interesting questions about the use of past glacier change information to reconstruct climate variations, which we do not address here. In the context of model initialization, this nonuniqueness impedes the reconstruction. We evaluated the candidates with a fitness function based on averaged geometry differences between the forward modeled and the observed state in t_{e}=2000. The threshold value ϵ=125 m^{2} was derived by assuming a typical error of 5 m in surface elevations and 10 m in glacier width, but how these values should be chosen depends on the specific glacier setting. Especially in cases where many of the candidate states have extremely small fitness values, a more strict acceptance criterion can help to narrow the results. On the other hand, an ϵ that is too strict could lead to none of the candidates fulfilling the criterion.
Due to uncertainties that derive from the model integration, the glacier state with the minimal fitness value is not always close to the synthetic experiment. As a more robust alternative, we propose using ${s}_{t}^{\mathrm{med}}$, the median state of the 5th percentile of all acceptable states as the best estimate. In Sect. 4.1 we compared the errors of both approaches. The median error of ${s}_{t}^{\mathrm{med}}$ is slightly smaller than that of ${s}_{t}^{\mathrm{min}}$, and the total range of absolute errors is smaller for ${s}_{t}^{\mathrm{med}}$ in 1850, too. Modeling the reconstructed initial states forward in time approximately 60 years leads to a rapid reduction of the error, and ${s}_{t}^{\mathrm{min}}$ performs a bit better than ${s}_{t}^{\mathrm{med}}$. By making use of the knowledge about the past climate, the number of candidates at later stages is, through this forward run, more constrained than by initializing them directly at a later time (see Appendix B for a more detailed description of the inverted approach at different times).
By comparing different fitness functions for the candidate evaluation, we showed that using limited information only (glacier area and glacier length) leads to an increase in the errors in 1850. This indicates what kind of observation is needed to be able to reconstruct past glacier states from today's state. The differences between the geometrybased evaluation and the areabased evaluation are small, but the differences to the lengthbased evaluation are significant. But this effect is also influenced by the spatial resolution of the model grid: a higher resolution of the grid would lead to more variability in fitness values and hence to a more precise initialization. At the same time, a higher resolution would increase the computational demands of the initialization method. We strongly recommend using either the geometry or the areabased fitness function for the evaluation. In this study, we only take the observation of the year t_{e} into account. Multitemporal outlines are likely to greatly reduce uncertainties at prior times.
Our results are relevant for future glacier evolution modeling studies, as they indicate that at least for some glaciers the time needed to converge to a similar evolution regardless of the 1850 state is comparatively short. Our study might also be useful to determine a good starting point of a past simulation, e.g., to improve the initialization date in Zekollari et al. (2019). A correlation analysis of the reconstructability and glacier characteristics showed the position of the ELA as well as the slope (especially in the lower part of the glacier) influence the reconstructability, whereas attributes like the glacier size do not have a strong impact. We could also show that the reconstructability measure correlates well with a separately obtained response time of the glacier.
Future work will include the application of the method to realworld cases, which will come with additional challenges. For example, we will have to consider the merging of neighboring glaciers when growing. Importantly, the effect of uncertainties in the boundary conditions (in particular the glacier bed and its outlines and uncertainties in the climate forcing) will have to be quantified. This also includes testing the influence of the choice of climate conditions on the accuracy of our method. Here again, the synthetic framework will be useful by allowing data denial and data alteration experiments. To ensure the robust reconstruction of realworld glacier states, additional changes and model developments are necessary. This includes the development of a glacierindividual calibration method for dynamical parameters (e.g., sliding parameter, creep parameter) as well as of the mass balance model.
The OGGM software together with initialization method are coded in the Python language and licensed under the GPLV3 free software license. The latest version of the OGGM code is available on GitHub (https://github.com/OGGM/oggm, Maussion et al., 2019b), the documentation is hosted on Read the Docs (http://oggm.readthedocs.io, last access: 10 December 2019), and the project web page for communication and dissemination can be found at http://oggm.org (last access: 10 December 2019). The OGGM version used in this study is v1.1. The code for the initialization module is available on GitHub (https://github.com/OGGM/initialization, Eis et al., 2019). The OGGM version used for this study is available in a permanent DOI repostiory (https://doi.org/10.5281/zenodo.2580277, Maussion et al., 2019c).
For the generation of the synthetic experiment state in 1850, we use a temperature bias of −1 K in order to create a relatively big glacier state. To justify the choice of this value, we have tested different temperature biases: the results are summarized in Fig. A1. This figure shows that applying positive or small negative temperature biases to the synthetic experiments results in large area differences to the RGI in 2000, and the total glacierized area in 1850 is also too small. The sample size is reduced because fewer glaciers fulfill the area threshold criteria of 0.01 km^{2}. Negative temperature biases that are too large also reduce the sample size because some runs fail (the glacier becomes larger than the underlying grid). The experiments with a temperature bias of −1, −1.25, or −1.5 K perform best regarding the area difference to the RGI in 2000. But only the experiment with the temperature bias of −1 K performs well regarding the estimation in 1850 of Zemp et al. (2006), where however it needs to be taken into account that the dots only represent a subset (the small glaciers in 2000 are missing) of the glaciers considered in Zemp et al. (2006).
We applied our method to different starting times (1850, 1855, …, 1965) to test how far one can go back in time to obtain a good initial state for this glacier. While this inverted setup is computationally very expensive, unfortunately it does not lead to improved results. See Fig. B1 for two different examples.
For each tested starting year, we determined the median state and conducted an uncertainty analysis (similar to the one in Sect. 4.1). We find that the uncertainties of the median states at the different starting points are higher than when performing the initialization for the year 1850 (only) and running this state forward in time. While this is counterintuitive, the main reason is that by starting in 1850 even with a very large number and range of candidates, the very unrealistic ones are quickly forced to converge by the boundary conditions (i.e., by climate), effectively reducing the number of potential candidates for a later date. In other words, we make use of our knowledge about past climate to reduce the number of candidates at each later stage. In realworld applications, results might be different since uncertainties in past climate are large. While this should be explored further, because of the computational cost it is hard to imagine an eventual applicability on the global or even large regional scale.
The supplement related to this article is available online at: https://doi.org/10.5194/tc1333172019supplement.
JE is the main developer of the initialization module and wrote most of the paper. BM and FM are the initiators of the OGGM project and helped to conceive this study. FM is the main OGGM developer and participated in the development of the initialization module.
The authors declare that they have no conflict of interest.
We thank the editor Andreas Vieli and the two anonymous referees for their comments that helped to improve our paper.
Ben Marzeion and Julia Eis were supported by the German Research Foundation (grant nos. MA 6966/11 and MA 6966/12).
The article processing charges for this openaccess publication were covered by the University of Bremen.
This paper was edited by Andreas Vieli and reviewed by two anonymous referees.
Auer, I., Böhm, R., Jurkovic, A., Lipa, W., Orlik, A., Potzmann, R., Schöner, W., Ungersböck, M., Matulla, C., Briffa, K., Jones, P., Efthymiadis, D., Brunetti, M., Nanni, T., Maugeri, M., Mercalli, L., Mestre, O., Moisselin, J.M., Begert, M., MüllerWestermeier, G., Kveton, V., Bochnicek, O., Stastny, P., Lapin, M., Szalai, S., Szentimrey, T., Cegnar, T., Dolinar, M., GajicCapka, M., Zaninovic, K., Majstorovic, Z., and Nieplova, E.: HISTALP—historical instrumental climatological surface time series of the Greater Alpine Region, Int. J. Climatol., 27, 17–46, https://doi.org/10.1002/joc.1377, 2007. a
Bach, E., Radić, V., and Schoof, C.: How sensitive are mountain glaciers to climate change? Insights from a block model, J. Glaciol., 64, 247–258, https://doi.org/10.1017/jog.2018.15, 2018. a
Bamber, J., Westaway, R. M., Marzeion, B., and Wouters, B.: The land ice contribution to sea level during the satellite era, Environ. Res. Lett., 13, 099502, https://doi.org/10.1088/17489326/aac2f0, 2018. a
Church, J., Clark, P., Cazenave, A., Gregory, J., Jevrejeva, S., Levermann, A., Merrifield, M., Milne, G., Nerem, R., Nunn, P., Payne, A., Pfeffer, W., Stammer, D., and Unnikrishnan, A.: Sea Level Change, book section 13, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1137–1216, https://doi.org/10.1017/CBO9781107415324.026, 2013. a
Colton, D. and Kress, R.: Inverse acoustic and electromagnetic scattering theory., vol. 93, Berlin, SpringerVerlag, 1992. a
Eis, J., Maussion, F., and Marzeion, B.: Initialization module of OGGM, available at: https://github.com/OGGM/initialization, last access: 10 December 2019. a
Farinotti, D., Huss, M., Bauder, A., Funk, M., and Truffer, M.: A method to estimate the ice volume and icethickness distribution of alpine glaciers, J. Glaciol., 55, 422–430, https://doi.org/10.3189/002214309788816759, 2009. a
Farinotti, D., Brinkerhoff, D. J., Clarke, G. K. C., Fürst, J. J., Frey, H., Gantayat, P., GilletChaulet, F., Girard, C., Huss, M., Leclercq, P. W., Linsbauer, A., Machguth, H., Martin, C., Maussion, F., Morlighem, M., Mosbeux, C., Pandit, A., Portmann, A., Rabatel, A., Ramsankaran, R., Reerink, T. J., Sanchez, O., Stentoft, P. A., Singh Kumari, S., van Pelt, W. J. J., Anderson, B., Benham, T., Binder, D., Dowdeswell, J. A., Fischer, A., Helfricht, K., Kutuzov, S., Lavrentiev, I., McNabb, R., Gudmundsson, G. H., Li, H., and Andreassen, L. M.: How accurate are estimates of glacier ice thickness? Results from ITMIX, the Ice Thickness Models Intercomparison eXperiment, The Cryosphere, 11, 949–970, https://doi.org/10.5194/tc119492017, 2017. a
Farinotti, D., Huss, M., Fürst, J. J., Landmann, J., Machguth, H., Maussion, F., and Pandit, A.: A consensus estimate for the ice thickness distribution of all glaciers on Earth, Nat. Geosci., 12, 168–173, https://doi.org/10.1038/s4156101903003, 2019. a
Gardner, A. S., Moholdt, G., Cogley, J. G., Wouters, B., Arendt, A. A., Wahr, J., Berthier, E., Hock, R., Pfeffer, W. T., Kaser, G., Ligtenberg, S. R. M., Bolch, T., Sharp, M. J., Hagen, J. O., van den Broeke, M. R., and Paul, F.: A Reconciled Estimate of Glacier Contributions to Sea Level Rise: 2003 to 2009, Science, 340, 852–857, https://doi.org/10.1126/science.1234532, 2013. a
Giesen, R. H. and Oerlemans, J.: Calibration of a surface mass balance model for globalscale applications, The Cryosphere, 6, 1463–1481, https://doi.org/10.5194/tc614632012, 2012. a
Giesen, R. H. and Oerlemans, J.: Climatemodel induced differences in the 21st century global and regional glacier contributions to sealevel rise, Clim. Dynam., 41, 3283–3300, https://doi.org/10.1007/s0038201317437, 2013. a
Goelzer, H., Nowicki, S., Edwards, T., Beckley, M., AbeOuchi, A., Aschwanden, A., Calov, R., Gagliardini, O., GilletChaulet, F., Golledge, N. R., Gregory, J., Greve, R., Humbert, A., Huybrechts, P., Kennedy, J. H., Larour, E., Lipscomb, W. H., Le clec'h, S., Lee, V., Morlighem, M., Pattyn, F., Payne, A. J., Rodehacke, C., Rückamp, M., Saito, F., Schlegel, N., Seroussi, H., Shepherd, A., Sun, S., van de Wal, R., and Ziemen, F. A.: Design and results of the ice sheet model initialisation experiments initMIPGreenland: an ISMIP6 intercomparison, The Cryosphere, 12, 1433–1460, https://doi.org/10.5194/tc1214332018, 2018. a
Goosse, H., Barriat, P.Y., Dalaiden, Q., Klein, F., Marzeion, B., Maussion, F., Pelucchi, P., and Vlug, A.: Testing the consistency between changes in simulated climate and Alpine glacier length over the past millennium, Clim. Past, 14, 1119–1133, https://doi.org/10.5194/cp1411192018, 2018. a
Gregory, J. M., White, N. J., Church, J. A., Bierkens, M. F. P., Box, J. E., van den Broeke, M. R., Cogley, J. G., Fettweis, X., Hanna, E., Huybrechts, P., Konikow, L. F., Leclercq, P. W., Marzeion, B., Oerlemans, J., Tamisiea, M. E., Wada, Y., Wake, L. M., and van de Wal, R. S. W.: TwentiethCentury GlobalMean Sea Level Rise: Is the Whole Greater than the Sum of the Parts?, J. Climate, 26, 4476–4499, https://doi.org/10.1175/JCLID1200319.1, 2013. a, b
Harris, I. C., Jones, P. D., Osborn, T., and Lister, D.: Updated hightresolution grids of montly climatic observations the CRU TS3.10 Dataset, Int. J. Climatol., 34, 623–642, https://doi.org/10.1002/joc.3771, 2014. a
Heimbach, P. and Bugnion, V.: Greenland icesheet volume sensitivity to basal, surface and initial conditions derived from an adjoint model, Ann. Glaciol., 50, 67–80, https://doi.org/10.3189/172756409789624256, 2009. a
Henderson, L. S. and Subbarao, K.: “Inverse Crime” and Model Integrity in Lightcurve Inversion applied to unresolved Space Object Identification, J. Astronaut. Sci., 64, 399–413, https://doi.org/10.1007/s4029501601051, 2017. a
Hock, R., Bliss, A., Marzeion, B., Giesen, R. H., Hirabayashi, Y., Huss, M., Radić, V., and Slangen, A. B. A.: GlacierMIP – A model intercomparison of globalscale glacier massbalance models and projections, J. Glaciol., 65, 453–467, https://doi.org/10.1017/jog.2019.22, 2019. a
Huss, M. and Hock, R.: A new model for global glacier change and sealevel rise, Front. Earth Sci., 3, 54, https://doi.org/10.3389/feart.2015.00054, 2015. a, b, c, d
Jacob, T., Wahr, J., Pfeffer, W. T., and Swenson, S.: Recent contributions of glaciers and ice caps to sea level rise, Nature, 482, 514–518, https://doi.org/10.1038/nature10847, 2012. a
Jarvis, A., Guevara, E., Reuter, H., and Nelson, A.: Holefilled SRTM for the globe: version 4: data grid, CGIARCSI, 2008. a
Kienholz, C., Rich, J. L., Arendt, A. A., and Hock, R.: A new method for deriving glacier centerlines applied to glaciers in Alaska and northwest Canada, The Cryosphere, 8, 503–519, https://doi.org/10.5194/tc85032014, 2014. a
Leclercq, P. W., Oerlemans, J., and Cogley, J. G.: Estimating the Glacier Contribution to SeaLevel Rise for the Period 1800–2005, Surv. Geophys., 32, 519, https://doi.org/10.1007/s1071201191217, 2011. a
Lee, V., Cornford, S. L., and Payne, A. J.: Initialization of an icesheet model for presentday Greenland, Ann. Glaciol., 56, 129–140, https://doi.org/10.3189/2015AoG70A121, 2015. a
Lüthi, M. P.: Transient response of idealized glaciers to climate variations, J. Glaciol., 55, 918–930, https://doi.org/10.3189/002214309790152519, 2009. a
Marzeion, B., Jarosch, A. H., and Hofer, M.: Past and future sealevel change from the surface mass balance of glaciers, The Cryosphere, 6, 1295–1322, https://doi.org/10.5194/tc612952012, 2012. a, b, c, d
Marzeion, B., Cogley, J. G., Richter, K., and Parkes, D.: Attribution of global glacier mass loss to anthropogenic and natural causes, Science, 345, 919–921, https://doi.org/10.1126/science.1254702, 2014. a, b
Marzeion, B., Leclercq, P. W., Cogley, J. G., and Jarosch, A. H.: Brief Communication: Global reconstructions of glacier mass change during the 20th century are consistent, The Cryosphere, 9, 2399–2404, https://doi.org/10.5194/tc923992015, 2015. a, b
Marzeion, B., Kaser, G., Maussion, F., and Champollion, N.: Limited influence of climate change mitigation on shortterm glacier mass loss, Nat. Clim. Change, 8, 305–308, https://doi.org/10.1038/s4155801800931, 2018. a
Maussion, F., Butenko, A., Champollion, N., Dusch, M., Eis, J., Fourteau, K., Gregor, P., Jarosch, A. H., Landmann, J., Oesterle, F., Recinos, B., Rothenpieler, T., Vlug, A., Wild, C. T., and Marzeion, B.: The Open Global Glacier Model (OGGM) v1.1, Geosci. Model Dev., 12, 909–931, https://doi.org/10.5194/gmd129092019, 2019a. a, b, c, d, e, f, g, h
Maussion, F., Rothenspieler, T., Dusch, M., Recinos, B., Vlug, A., Marzeion, B., Landmann, J., Eis, J., Bartholomew, S. L., Champollion, N., Gregor, P., Butenko, A., Smith, S., and Oberrauch, M.: OGGM, v1.1, available at: https://github.com/OGGM/oggm, last access: 10 December 2019b. a
Maussion, F., Rothenspieler, T., Dusch, M., Recinos, B., Vlug, A., Marzeion, B., Landmann, J., Eis, J., Bartholomew, S. L., Champollion, N., Gregor, P., Butenko, A., Smith, S., and Oberrauch, M.: OGGM/oggm: v1.1, Zenodo, https://doi.org/10.5281/zenodo.2580277, 2019c. a
Mosbeux, C., GilletChaulet, F., and Gagliardini, O.: Comparison of adjoint and nudging methods to initialise ice sheet model basal conditions, Geosci. Model Dev., 9, 2549–2562, https://doi.org/10.5194/gmd925492016, 2016. a
Oerlemans, J.: Climate Sensitivity of Franz Josef Glacier, New Zealand, as Revealed by Numerical Modeling, Arct. Alpine Res., 29, 233–239, available at: https://www.tandfonline.com/doi/abs/10.1080/00040851.1997.12003238 (last access: 10 December 2019), 1997. a
Oerlemans, J.: Glaciers and Climate Change, Swets and Zeitlinger, Lisse, 2001. a
Pfeffer, W. T., Arendt, A. A., Bliss, A., Bolch, T., Cogley, J. G., Gardner, A. S., Hagen, J. O., Hock, R., G., K., Kienholz, C., Miles, E. S., Moholdt, G., Mölg, N., Paul, F., Radić, V., Rastner, P., Raup, B. H., Rich, J., Sharp, M. J., and the Randolph Consortium: The Randolph Glacier Inventory: a globally complete inventory of glaciers, J. Glaciol., 60, 537–551, https://doi.org/10.3189/2014JoG13J176, 2014. a, b, c
Radić, V. and Hock, R.: Regionally differentiated contribution of mountain glaciers and ice caps to future sealevel rise, Nature Geoscience, 4, 91–94, https://doi.org/10.1038/ngeo1052, 2011. a
Radić, V. and Hock, R.: Glaciers in the Earth's Hydrological Cycle: Assessments of Glacier Mass and Runoff Changes on Global and Regional Scales, Surv. Geophys., 35, 813–837, https://doi.org/10.1007/s107120139262y, 2014. a
Raup, B., Racoviteanu, A., Khalsa, S. J. S., Helm, C., Armstrong, R., and Arnaud, Y.: The GLIMS geospatial glacier database: A new tool for studying glacier change, Global Planet. Change, 56, 10 –110, https://doi.org/10.1016/j.gloplacha.2006.07.018, 2007. a
Slangen, A. B. A., Adloff, F., Jevrejeva, S., Leclercq, P. W., Marzeion, B., Wada, Y., and Winkelmann, R.: A Review of Recent Updates of SeaLevel Projections at Global and Regional Scales, Surv. Geophys., 38, 385–406, https://doi.org/10.1007/s1071201693742, 2017a. a
Slangen, A. B. A., Meyssignac, B., Agosta, C., Champollion, N., Church, J. A., Fettweis, X., Ligtenberg, S. R. M., Marzeion, B., Melet, A., Palmer, M. D., Richter, K., Roberts, C. D., and Spada, G.: Evaluating Model Simulations of TwentiethCentury Sea Level Rise. Part I: Global Mean Sea Level Change, J. Climate, 30, 8539–8563, https://doi.org/10.1175/JCLID170110.1, 2017b. a
van Pelt, W. J. J., Oerlemans, J., Reijmer, C. H., Pettersson, R., Pohjola, V. A., Isaksson, E., and Divine, D.: An iterative inverse method to estimate basal topography and initialize ice flow models, The Cryosphere, 7, 987–1006, https://doi.org/10.5194/tc79872013, 2013. a
WCRP Global Sea Level Budget Group: Global sealevel budget 1993–present, Earth Syst. Sci. Data, 10, 1551–1590, https://doi.org/10.5194/essd1015512018, 2018. a
WGMS: Fluctuations of Glaciers Database, World Glacier Monitoring Service, Zurich, Switzerland, https://doi.org/10.5904/wgmsfog201811, 2018. a
Wouters, B., Gardner, A. S., and Moholdt, G.: Global Glacier Mass Loss During the GRACE Satellite Mission (2002–2016), Front. Earth Sci., 7, 96, https://doi.org/10.3389/feart.2019.00096, 2019. a, b
Zekollari, H. and Huybrechts, P.: On the climate–geometry imbalance, response time and volume–area scaling of an alpine glacier: insights from a 3D flow model applied to Vadret da Morteratsch, Switzerland, Ann. Glaciol., 56, 51–62, https://doi.org/10.3189/2015AoG70A921, 2015. a
Zekollari, H., Huss, M., and Farinotti, D.: Modelling the future evolution of glaciers in the European Alps under the EUROCORDEX RCM ensemble, The Cryosphere, 13, 1125–1146, https://doi.org/10.5194/tc1311252019, 2019. a, b
Zemp, M., Haeberli, W., Hoelzle, M., and Paul, F.: Alpine glaciers to disappear within decades?, Geophys. Res. Lett., 33, L13504, https://doi.org/10.1029/2006GL026319, 2006. a, b, c, d
Zemp, M., Frey, H., GärtnerRoer, I., Nussbaumer, S. U., Hoelzle, M., Paul, F., Haeberli, W., Denzinger, F., Ahlstrøm, A. P., Anderson, B., Bajracharya, S., Baroni, C., Braun, L. N., Cáceres, B. E., Casassa, G., Cobos, G., Dávila, L. R., Delgado, G. H., Denuth, M. N., Espizua, L., Fischer, A., Fujita, K., Gadek, B., Ghazanfar, A., Hagen, J. O., Holmlund, P., Karimi, N., Li, Z., Pelto, M., Pitte, P., Popovnin, V., Portocarrero, C. A., Prinz, R., Sangewar, C., Severskiy, I., Sigurdsson, O., Soruco, A., Usubaliev, R., and Vincent, C.: Historically unprecedented global glacier decline in the early 21st century, J. Glaciol., 61, 745–762, https://doi.org/10.3189/2015JoG15J017, 2015. a
Zemp, M., Huss, M., Thibert, E., Eckert, N., McNabb, R., Huber, J., Barandun, M., Machguth, H., Nussbaumer, S. U., GärtnerRoer, I., Thomson, L., Paul, F., Maussion, F., Kutuzov, S., and Cogley, J. G.: Global glacier mass changes and their contributions to sealevel rise from 1961 to 2016, Nature, 568, 382–386, https://doi.org/10.1038/s4158601910710, 2019. a
Instead, it is also possible to choose ${s}_{t}^{\mathrm{med}}$ for the uncertainty analyses, but this would require acceptance criteria for the fitness functions J_{A} and J_{l}, which would influence the state. For simplification, we choose the state with the minimal fitness function.
We also tested a step of 0.5 K, leading to a larger sample size but no significant change to the correlation analysis and our results. Thus, we kept the 1 K step change here.
 Abstract
 Introduction
 Methods
 Test site and input data
 Results
 Hardware requirements and performance
 Discussion and conclusions
 Code availability
 Appendix A: Temperature bias for the synthetic experiments
 Appendix B: Initialization at different starting times
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Methods
 Test site and input data
 Results
 Hardware requirements and performance
 Discussion and conclusions
 Code availability
 Appendix A: Temperature bias for the synthetic experiments
 Appendix B: Initialization at different starting times
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement