Initialization of a global glacier model based on present-day glacier geometry and past climate information: an ensemble approach

Abstract. 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 present-day 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 ill-posed 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 present-day 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 present-day 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 real-world applications.



Introduction
Glaciers contributed significantly to past sea-level rise (SLR; e.g. Gregory et al., 2013;Slangen et al., 2017a;Cazenave and et al., 2018) and they will continue to be a major contributor in the 21 th century (e.g. Church et al., 2013;Slangen et al., 2017b).
A large fraction of this contribution will be caused by the ongoing adjustment to the past climate (Marzeion et al., 2014. 20 Reconstructions of past glacier mass change are therefore not only necessary to determine the budget of past sea-level change  and to increase the confidence in projections (by allowing to quantify the agreement with observations, 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 modelled velocities or surface elevations.
In this study, we aim to bring elements of response to the following questions: to which degree does the past evolution of a glacier constrain its present day geometry? 5 how much information does the present day glacier geometry contain about its past states?
is it possible to reconstruct past glacier states from the partial information available to us?
To this aim, we present a new method estimating past glacier states and apply it to synthetic numerical experiments. After introducing the relevant features of the Open Global Glacier Model (OGGM; Maussion et al., 2019) in Sect. 2.1, we describe the design of the experiments in Sect. 2.3. The synthetic framework serves to test the skill of our approach in a surrogate 10 model world where everything is known, and allows 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 3 steps: generation of plausible glacier states, identification of glacier state candidates, and their evaluation based on the misfit between the modelled and the observed surface elevation at the year of the observation. We applied our approach to 2621 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 15 length, area or geometry) is shown in Sect. 4.2. Finally, we summarize the results and discuss the method's limitations and its applicability to real-case studies, as well as needed and possible future developments in Sect. 6.

The Open Global Glacier Model
The Open Global Glacier Model (OGGM; Maussion et al., 2019) is an open source numerical framework that allows the 20 modelling 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 also large glacier states can be generated. The resolution of the map topography dx depends on the size of the glacier and is constrained to 10 m ≤ dx ≤ 200 m. After the preprocessing, 25 glacier centerlines are computed using a geometrical routing algorithm (adapted from Kienholz et al., 2014). They are then considered as glacier flowlines, 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 flowline coordinates are then obtained from the underlying topography file and glacier section widths are computed by intersecting the flowline's normal 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 30 with a mass-conservation approach (Farinotti et al., 2009(Farinotti et al., , 2017. Information on bed topography at each grid point results from the calculated ice thickness and the surface elevation. From this, the glacier's 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 model then computes the glacier mass balance 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 for more specialized applications.

5
To force the mass balance model with a randomized climate time series, a window size h and a center year y 0 need to be set first. All years ∈ [y 0 − h−1 2 , y 0 + h−1 2 ] are then shuffled in the next step. Additionally, it is possible to set a temperature bias β, which shift all values of the temperature series. The objective of forcing the mass balance model with randomized climate is to easily produce a great number of realistic climate settings, representative of a given time period. The dynamical flowline model can then be used to determine the evolution of the glacier under certain climate forcings by solving the shallow ice 10 approximation.
For more details concerning the glacier model, please refer to Maussion et al. (2019) and http://docs.oggm.org.

Problem description
Here, we define a glacier state (hereinafter referred as state) as follows: Definition 1. Let m ∈ N be the total number of grid points of all flowlines of a glacier. Then s t = (z t , w t , b) is a glacier state 15 at time t, with surface elevation z t ∈ R m + , widths w t ∈ R m + , and bed topography b ∈ R m + . The set S ti = {s t |t = t i } 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:

20
Given an initial state s t0 ∈ S t0 at time t 0 , the state s t ∈ S t at time t > t 0 can be computed by: where G past : S t0 → S t is an operator representing the equations of OGGM, using known climate time series as boundary condition. 25 For inverse problems, the solution is known by direct observations: s te = s obs te , whereas the desired initial state s t0 is unknown. The inverse problem consists of finding the initial state s t0 ∈ S t0 , such that the forward modelled solution at time t e fits the observations from the same year t e : Unfortunately, we do not have an explicit formulation for G −1 past in our case. A backwards reconstruction is impeded by the 30 non-linear 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 modelled state is as close as possible to the observation: This function calculates the averaged difference in surface elevation and width between the observed and forward modelled glacier state. 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's 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(s t0 ). As uncertainties of the model can safely be assumed to be 10 larger than the differences between those states at time t e , it is impossible to identify the global minimum of j(s t0 ). I.e., the solution of our inverse problem is non-unique.
The objective of our approach is therefore to identify the set S t0 of all states which correspond to the observed state s obs te within a given uncertainty after being modelled forward. We call this condition acceptance criterion: The function J(s t0 ) is called in in the following fitness function. Assuming a vertical error of 5 m in x and an horizontal error of 10 m in w, we propose to set = (5m) 2 + (10m) 2 = 125m 2 . These numbers can be changed easily, and in a real-world application should be based on the vertical uncertainty of the used DEM and the horizontal uncertainty of the used outline.
All states s t0 ∈ S 125 t0 that have a fitness value smaller than one are called acceptable states. The first expectation would be that the glacier with the smallest fitness value, is be the best solution. However, due to uncertainties that derive from the model 20 itself, this is not always the case. As an alternative, we determine the 5th quantile of all states in S t0 . This set contains the best solutions of all acceptable states referred to the fitness value. 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.

Synthetic experiments
In order to be able to evaluate the results of the initialization method, we use synthetic experiments. To this end, we use OGGM 25 to initialize a model run. We apply a random climate scenario (window size h = 31 years, and center year y 0 = t 0 ) and run the model 400 years forward. The temperature bias is set to β = −1K to ensure that a relatively big glacier state is created, because this the case for most real glaciers in 1850. The resulting state is defined to be the state in year t 0 . We model this state forward, applying the past climate time series from t 0 until t e (here: 2000) and obtain the observed state of the synthetic experiment.
Through this procedure, a time series of glacier states is created, which in the following we will try to recover applying the 30 initialization method described below, but only using information about the observed state from the synthetic experiment in t e . Thanks to the initial temperature bias of β = −1K, these synthetic states in t e are close to the real observed states in 2000 (total area difference of about 1%, but individual glaciers can vary). In the following we call the state derived from the synthetic experiment s exp t . Figure 1 shows the surface elevations in 1850 and 2000 of the synthetic experiment for the Guslarferner as an example. Black line indicates the bed rock and the dashed line the ice surface of the synthetic experiment.

5
We will discuss the limitations of the usage of the synthetic experiments in the context of real-world applications in Sect. 6.

Reconstruction of initial glacier states
Our initialization method consists of three main steps: generation of a set of physically plausible glacier states S t0 , identification of glacier candidates s t0 ∈ S t0 , and their evaluation based on the fitness function J(s t0 ) (see Sec. 2.2).

Generation of glacier states 10
In a first step, we generate a set of different, physically plausible states. 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 (see Sect. 2.3). This procedure has the advantage that a realistic climate representative of a given time period can be created, while interannual variability is uncorrelated to that of the period. Hence, all generated climate models differ from each other, but represent the climate conditions around t 0 15 at the same time. Additionally, we vary the temperature bias β to create further variations. We start with temperature biases  β ∈ [−2, 2] K. If β = 2 K is not large enough to result in a glacier with zero ice thickness, higher values will be used. If β = −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's boundary), smaller values will be used. We use these climate conditions to force a 400-year glacier simulation. 400 years turned out to be a sufficient time to generate a large number of different states. This model run is initialized with the only known state, which is the actual observed glacier topography, taken from the RGI (see Fig. 2a).  (Jóhannesson et al., 1989). Glacier volume increases with decreasing temperature biases. Thus, we choose the upper ten volume trajectories, corresponding to the lowest temperature biases, to determine t stag . To this end, we smooth each of the ten curves 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 fluctuate around an equilibrium value is a potential glacier state candidate (in the following referred as candidate). This holds true for all states s t with t > t stag . In order to avoid testing 5 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.

Evaluation
The last step evaluates all previously selected candidates. Each candidate is used as 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 10 for the initial condition and exactly the same climate time series. Afterwards, we compare the resulting modelled states s te with the observed state s obs te by applying the fitness function J(s t0 ) (Eq. 5). This function calculates the averaged difference between the glacier geometries at the grid points, more specifically between the surface elevations z te and the widths w te , of the observed and the modelled glacier. In Fig. 2c the candidates are colored by their fitness value. The resulting output of this method is the temporal evolution of all candidates from t 0 until t e sorted by their fitness value. 15

Test site and Input Data
We tested our approach on alpine glaciers. The glacier outlines are derived from the Randolph Glacier Inventory (RGI v6.0, region 11; Pfeffer et al., 2014). For each of the glacier outlines a transverse Mercator projection, which is centered on the glacier, is defined in order to preserve map projection properties (e.g. area, distances, angles). Next, the topographical data is Here we show (i) the results for two example glaciers in 1850, as well as an error analysis for all tested glaciers, and (ii) the influence of the choice of the fitness function on the quality of our results.

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   5 show the results of the Guslarferner, as an example glacier. A second case, the Hintereisferner, is shown in Fig. 5 and 6.
Especially the result of the Guslarferner shows clearly that the determination of past states is not unique (see Fig. 3 The results of Hintereisferner are different to the results from Guslarferner (see Fig. 5). The fitness values of most candidates 25 are large compared to the ones of the Guslarferner. Only a few candidates have extremely small fitness values. In the case of the Hintereisferner, 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 modelled states and the observed one in 2000 are larger. The fitness values of all candidates range between 1.2×10 −3 and 47. 68 candidates fulfill the acceptance criterion (Eq. 5). Figure 6 shows that the acceptance criterion in this case confines the where v exp (  Both states fit well the synthetic experiment. In many cases, s med t is equal to s min t , but for some glaciers either s min t or s med t have a clearly better performance. In all cases, the uncertainties quickly reduce after around 1880 to 1900. The error analysis also shows that our method has a slight tendency to overestimate the glacier volume, and that the error distribution is strongly skewed (see Figure 9a).

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 reconstructabilty of past glacier states. Thus, we additionally evaluate the candidates by only using information about the glacier area or length.
For the glacier area based evaluation, we used the following fitness function: where A te 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 in our test site, we evaluate the 200 candidates also 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 5 volume error to the synthetic experiment. Figure 10 shows the relative errors of all three evaluation methods. Figure 9a shows the distribution of the relative errors of the 5th-95th percentile in 1850, whereas the the evolution from 1850-2000 of the median and the 5th-95th percentile of the relative errors are shown in Fig. 9b. The more information is taken into account for the evaluation, the smaller are the errors.
The greatest uncertainties are associated with using the glacier length based fitness function (Eq. 9), whereas the differences error in 1850 of the glacier length evaluation has with 137% the worst performance. This holds also true regarding the values of the 95th precentile; 95% of the tested cases have in 1850 a relative volume error smaller than 1008%, if the length based fitness function is used for the evaluation. For the area based evaluation this value is 329% and for the geometry based fitness function 95% of the tested glaciers have an error smaller than 286%. This shows that the advantage using the geometry instead of the glacier area to evaluate the candidates is not very high; both evaluations shows a very good performance. Especially if 5 the states are modelled forward (e.g. to 1900), both approaches perform well. However, it is not advisable to use the glacier length based evaluation.

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). Our method requires to run hundreds of dynamical model runs for each single glacier and as described

Discussion and conclusions
In this study, a new method to initialize past glacier states is presented and applied in 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 differentiation between the different error sources in a real-world application is difficult and thus an evaluation of the skill of the method would be difficult, too. For real world applications, we 5 will have to validate the reconstructions against past outlines derived from moraines, historical maps or remote sensing (e.g. as provided by GLIMS, Raup et al., 2007).
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 modelled initial glacier states in 1850 do not correspond to reality either. The past states determined in this study only can serve to verify the functionality of the developed method.

10
The results in Sect. 4.1 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 don't address here. In the context of model initialization, this non-uniqueness is impeding the reconstruction. We evaluated the candidates with a fitness function based on averaged geometry differences between the forward modelled and the observed state in t e = 2000. The threshold value = 125 m 2 was derived by assuming a typical error of 5 15 m in surface elevations and 10 m in glacier width, but how these values should be chosen depends on the situation. 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. Otherwise, an that is too strict could lead to none of the candidates fulfilling the criterion.
Due to uncertainties that derive from the model itself, the glacier state with the minimal fitness value is not always close to the synthetic experiment. As a more robust alternative, we propose to use s med t , 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. On the one hand, the median In Sect. 4.2 we compare different fitness functions for the candidate evaluation. We show that using limited information only (glacier area and glacier length) lead to an increase of the errors in 1850. The differences between the geometry based eval-10 uation and the area based evaluation are small, but the differences to the length based evaluation are significant. This reflects that glacier states with the same length could differ strongly in volume and area. But this effect is also influenced by the spatial resolution of the model grid: a higher resolution of the grid would lead 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 to use either the geometry or the area based fitness function for the evaluation. In this study, 15 we only take the observation of the year t e into account. Multi-temporal outlines would 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. A correlation analysis of the estimated reconstruction spread and error with glacier characteristics (e.g. size, slope, ice thickness distribution per altitude, 20 climatic setting, etc.) may lead to improved understanding of the factors that influence the reconstructability of a glacier's past states, and will be the topic of a separate study.
Future work will also include the application of the method on real-world cases, which will come with additional challenges.
For example, we will have to consider the merging of neighboring glaciers when growing (a work already under way, Dusch, 2018). Importantly, the effect of uncertainties in the boundary conditions (in particular the glacier bed, its outlines and un-25 certainties in the climate forcing) will have to be quantified. Here again, the "surrogate world" framework will be useful by allowing data denial and data alteration experiments. To ensure the robust reconstruction of real-world glacier states additional changes and model developments are necessary, but our study is a first important step in this direction.
Code availability. 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), the doc-30 umentation is hosted on ReadTheDocs (http://oggm.readthedocs.io), and the project webpage for communication and dissemination can be found at http://oggm.org. 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).