The Open Global Glacier Model ( OGGM ) v 1 . 0

Despite of their importance for sea-level rise, seasonal water availability, and as source of geohazards, mountain glaciers are one of the few remaining sub-systems of the global climate system for which no globally applicable, open source, communitydriven model exists. Here we present the Open Global Glacier Model (OGGM, www.oggm.org), developed to provide a modular and open source numerical model framework for simulating past and future change of any glacier in the world. The mod5 elling chain comprises data downloading tools (glacier outlines, topography, climate, validation data), a preprocessing module, a mass-balance model, a distributed ice thickness estimation model, and an ice flow model. The monthly mass-balance is obtained from gridded climate data and a temperature index melt model. To our knowledge, OGGM is the first global model explicitly simulating glacier dynamics: the model relies on the shallow ice approximation to compute the depth-integrated flux of ice along multiple connected flowlines. In this paper, we describe and illustrate each processing step by applying the 10 model to a selection of glaciers before running global simulations under idealized climate forcings. Even without an in-depth calibration, the model shows a very realistic behaviour. We are able to reproduce earlier estimates of global glacier volume by varying the ice dynamical parameters within a range of plausible values. At the same time, the increased complexity of OGGM compared to other prevalent global glacier models comes at a reasonable computational cost: several dozens of glaciers can be simulated on a personal computer, while global simulations realized in a supercomputing environment take up to a few hours 15 per century. Thanks to the modular framework, modules of various complexity can be added to the codebase, allowing to run new kinds of model intercomparisons in a controlled environment. Future developments will add new physical processes to the model as well as tools to calibrate the model in a more comprehensive way. OGGM spans a wide range of applications, from ice-climate interaction studies at millenial time scales to estimates of the contribution of glaciers to past and future sea-level change. It has the potential to become a self-sustained, community driven model for global and regional glacier evolution. 20 1 Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2018-9 Manuscript under review for journal Geosci. Model Dev. Discussion started: 23 March 2018 c © Author(s) 2018. CC BY 4.0 License.

interaction studies at millenial time scales to estimates of the contribution of glaciers to past and future sea-level change.It has the potential to become a self-sustained community-driven model for global and regional glacier evolution.

Introduction
Glaciers constitute natural low-pass filters of atmospheric variability.They allow people to directly perceive slow changes of the climate system, that would otherwise be masked by short-term noise in human perception.Since glaciers form prominent features of many landscapes, shrinking glaciers have become an icon of climate change.
However, impacts of glacier change -whether growth or shrinkage -go far beyond this sentimental aspect: glaciers are important regulators of water availability in many regions of the world (Kaser et al., 2010;Huss, 2011;Immerzeel et al., 2012), and retreating glaciers can lead to increased geohazards (see Richardson and Reynolds, 2000, for an overview).Even though the ice mass stored in glaciers is small compared to the Greenland and Antarctic ice sheets (< 1%), glacier melt has contributed significantly to past sea-level rise (SLR; e.g.Cogley, 2009;Leclercq et al., 2011;Marzeion et al., 2012b;Gardner et al., 2013).
They probably have been the biggest single source of observed SLR since 1900 and they will continue to be a major source of SLR in the 21st century (e.g.Church et al., 2013;Gregory et al., 2013).
It is therefore a pressing task to improve the knowledge of how glaciers change when subjected to climate change, both natural and anthropogenic (Marzeion et al., 2014a).The main obstacle to achieve progress is a severe undersampling problem: direct glaciological measurements of mass balances have been performed on ∼300 glaciers world wide (≈ 0.1% of all glaciers on Earth).The number of glaciers on which these types of measurements have been carried out for time periods longer than 30 years, i.e. over periods that potentially allow for the detection of a climate change signal, is one order of magnitude smaller (Zemp et al., 2009).Length variations of glaciers have been observed for substantially longer periods of time (Oerlemans, 1994(Oerlemans, , 2005)).These variations are, however, much more difficult to understand, as large glacier length fluctuations may arise from intrinsic climate variability (Roe and O'Neal, 2009;Roe, 2011).Data obtained by remote sensing allow for gravimetric assessments of ice mass change or volume change estimates obtained by differencing digital elevation models.Unfortunately, they are only available for the past decade (e.g.Gardner et al., 2013).
During the past few years, great progress has been made in methods to model glaciers globally (Radić andHock, 2011, 2014;Giesen andOerlemans, 2012, 2013;Marzeion et al., 2012aMarzeion et al., , b, 2014a, b;, b;Huss and Hock, 2015).While these approaches yield consistent results at the global scale, all of them suffer from greater uncertainties at the regional and local scales.These stem from the great level of abstraction of the key processes (Marzeion et al., 2012b(Marzeion et al., , 2014b)), from the need to spatially interpolate model parameters (Radić andHock, 2011, 2014;Giesen andOerlemans, 2012, 2013), and from uncertainties of the boundary and initial conditions.All models lack ice dynamics, most lack frontal ablation (with the exception of Huss and Hock, 2015), and all lack modulation of the surface mass balance by debris cover and snow redistribution (wind and avalanches).Only one model (Marzeion et al., 2012b) was able to provide estimates of past glacier volume changes for the 20th century.None of these models is open-source.
Mountain glaciers are one of the few remaining subsystems of the global climate system for which no globally applicable, open source, community-driven model exists.The ice sheet modelling community shows a better example, with models such as the Parallel Ice Sheet Model (Winkelmann et al., 2011), Elmer/Ice (http://elmerice.elmerfem.org/),Glimmer-CISM (https: //csdms.colorado.edu/wiki/Model:Glimmer-CISM),or ISSM (https://issm.jpl.nasa.gov/).These models have been applied to mountain glaciers as well, but cannot be applied globally out-of-the-box.While the atmospheric modelling community has a long tradition of sharing models (e.g. the Weather Research and Forecasting model, or WRF) or comparing them (e.g. the Coupled Model Intercomparison Project or CMIP), recent initiatives originating from the glaciological community show a new willingness to better coordinate global research efforts following the CMIP example (e.g. the Glacier Model Intercomparison Project1 or the Glacier Ice Thickness Estimation Working Group2 ).
In the recent past, great advances have been made in the global availability of data and methods relevant for glacier modelling, spanning glacier outlines (Pfeffer et al., 2014), automatized glacier centerline identification (e.g., Kienholz et al., 2014), bedrock inversion methods (e.g.Huss and Farinotti, 2012), and global topographic data sets (e.g.Farr et al., 2007).Taken together, these advances now allow the ice dynamics of glaciers to be simulated at the global scale, provided that adequate modelling platforms are available.In this paper, we present the Open Global Glacier Model (OGGM), developed to provide a modular and open source numerical model framework for consistently simulating past and future global scale glacier change.
Global not only in the sense of leading to meaningful results for all glaciers combined, but also for any small ensemble of glaciers, e.g. at the headwater catchment scale.Modular to allow different approaches to the representation of ice flow and surface mass balance to be combined and compared against each other.Open source so that the code can be read and used by anyone and so that new modules can be added and discussed by the community, following the principles of open governance.
Consistent between past and future in order to provide uncertainty measures at all realisable scales.This paper describes the basic structure and primordial assumptions of the model (as of version 1.1).We present the results of a series of single glacier and global simulations demonstrating the model's usage and potential.This will be followed by a description of the software requirements and the testing framework.Finally, we will discuss the potential for future developments that could be conducted by any interested research team.

Fundamental principles
The starting point of OGGM is the Randolph Glacier Inventory (RGI; RGI Consortium, 2017;Pfeffer et al., 2014): our goal is to simulate the past and future evolution of all of the 216,502 inventoried glaciers worldwide (as of RGI V6).This "glacier centric" approach is the one followed by most global and regional models to date; its advantages and disadvantages will be discussed in Sect.3.6.4.Provided with the glacier outlines, topographical and climate data at reasonable resolution and accuracy, the model should be able to (i) provide a local map of the glacier including topography and hypsometry, (ii) estimate the glacier's total ice volume and compute a map of the bedrock topography, (iii) compute the surface climatic mass balance and (if applicable) at its front via frontal ablation, (iv) simulate the glacier's dynamical evolution under various climate forcings, and (v) provide an estimate of the uncertainties associated with the modelling chain.
For each of these steps, several choices are possible regarding the input data to be used, the numerical solver or the parameterisations to be applied.Any given choice is driven by subjective considerations about data availability, the estimated accuracy of boundary conditions (such as topography), and by technical considerations such as the available computational resources.
In this paper we present one way to realize these steps using OGGM, which to date is in our opinion the best compromise between model complexity, data availability and computational effort.The OGGM software, however, is built in such a way that future improvements and new approaches can be implemented, tested, and applied at minimal cost by ourselves or a larger community.

Example workflow
We illustrate with an example how the OGGM workflow is applied to Tasman glacier, New Zealand (Fig. 1).Here we describe shortly the purpose of each processing step, and more details will be provided in Sect.3: Preprocessing The glacier outlines extracted from the RGI are projected onto a local gridded map of the glacier (Fig. 1a).
Depending on the glacier's location, a suitable source for the topographical data is downloaded automatically (here SRTM) and interpolated to the local grid.The map's spatial resolution depends on the size of the glacier (here, 150 m).
Flowlines The glacier centerlines are computed using a geometrical routing algorithm (adapted from Kienholz et al., 2014, Fig. 1b), filtered and slightly modified to become glacier flowlines with a fixed grid spacing.

Catchment areas and widths
The geometrical widths along the flowlines are obtained by intersecting the normals at each grid point with the glacier outlines and the tributaries' catchment areas (Fig. 1c).Each tributary and the main flowline has a catchment area, which is then used to correct the geometrical widths so that the flowline representation of the glacier is in close accordance with the actual altitude-area distribution of the glacier (Fig. 1d, note that the normals are now corrected and centred).
Climate data and mass balance Gridded climate data (monthly temperature and precipitation) are interpolated to the glacier location and temperature is corrected for altitude using a linear gradient.These climate time series are used to compute the glacier mass balance at each flowline's grid point for any month in the past.
Ice thickness inversion Using the mass balance data computed above and relying on mass-conservation considerations, an estimate of the ice flux along each glacier cross-section can be computed.By making assumptions about the shape of the cross-section (parabolic or rectangular) and using the physics of ice flow, the model computes the thickness of the glacier along the flowlines and the total volume of the glacier (Fig. 1e).
Glacier evolution A dynamical flowline model is used to simulate the advance and retreat of the glacier as a response of the surface mass-balance forcing.Here (Fig. 1f), a 100-yr long random climate sequence leads to a glacier advance.

Model structure
The OGGM model is built around the notion of tasks, which have to be applied sequentially to single or a set of glaciers.There are two types of tasks: Entity tasks are tasks which are applied on single glaciers individually and do not require information from other glaciers (this encompasses the majority of OGGM's tasks).Most often they need to be applied sequentially (for example, it is not possible to compute the centerlines without having read the topographical data first).
Global tasks are tasks which are run on a set of glaciers.This encompasses the calibration and validation routines, which need to gather data across a number of reference glaciers.
This model structure has several advantages: the same entity task can be run in parallel on several glaciers at the same time, and they allow a modular workflow.Indeed, a task can seamlessly be replaced by another similar one, as long as the required input and output formats are agreed upon beforehand.The output of each task is made persistent by storage on disk, allowing a later use by a subsequent task, even in a separate run or on another machine.For example, the preprocessing tasks store the topography data in a netCDF file, which is then read by the centerlines task, which itself writes it's output in a vector file format.
In this paper we will refrain from naming the tasks by their function name in the code, as these are likely to change in the future and are sometimes organised in a non-trivial way as a result of implementation details.The next section therefore is called "Modules", where each module can be seen as a collection of tasks developed towards a certain goal.

Modules
The modules are described in the order in which they are applied for a model run.When we provide a specific value for a model parameter in the text, we refer to the model's default parameter value: it can be changed by the user at runtime.

Preprocessing
The objective of the preprocessing module is to set up the geographical input data for each glacier (the glacier outlines and the local topography).First, a Cartesian local map projection is defined: we use a local Transverse Mercator projection centred on the glacier.Then, a suitable topographical data source is chosen automatically, depending on the glacier's location.Currently we use: the Shuttle Radar Topography Mission (SRTM) 90m Digital Elevation Database v4.1 (Jarvis et al., 2008) for all locations in the [60 • S; 60 • N] range (data acquisition: 2000) the Greenland Mapping Project (GIMP) Digital Elevation Model (Howat et al., 2014)  for global glacier modelling will be the topic of a follow-up study3 .
The spatial resolution of the target local grid depends on the size of the glacier: the default is to use a square relation to the glacier size (dx = aS 1 2 with a = 14 and S the area of the glacier in km 2 ) clipped to a predefined minimum (10 m) and maximum (200 m) value.After the interpolation to the target grid, the topography is smoothed with a Gaussian filter of 250 m radius (this value does not change with the local glacier map resolution because it is meant to be applied to the original DEM, not the interpolated one).This smoothing is driven by practical considerations, since the model becomes unstable if the boundary conditions are too noisy (see also Bahr et al., 2014, for a discussion about the unavoidable trade-off between resolution and accuracy).

Flowlines and catchments
The glacier centerlines are computed following an algorithm developed by Kienholz et al. (2014) and adapted for our purposes.
This algorithm was chosen because it allows one to compute multiple centerlines and to define a main branch fed by any number of tributaries.In general we found the method to be very robust, although some glaciers obviously won't have the optimal number of centerlines, with either too many (frequent in the case of large cirque glaciers) or not enough (some tributary branches have no centerlines).These errors however are assumed to play a relatively minor role compared to other uncertainties in the model chain.
In the model semantics, the original "centerlines" are then converted to "flowlines": the points defining the line geometries are interpolated to be equidistant from each other (the default spacing along the line is twice that of the local glacier map, i.e. varying between 20 m and 400 m depending on glacier size), and the tail of the tributaries are cut before reaching their descendant (see the differences between Figs. 1b and c, or between Figs. 2a and b).Each grid point's elevation is obtained from the underlying topography.By construction, upslope trajectories or sinks along the flowline are rare: this can still occur when the glacier outlines are poorly defined or because of errors in the gridded topography.In these cases, we interpolate the heights (in the case of a deepening) or cut the first grid points of the line (in case of an upslope starting from the flowline's head) until only positive slopes larger than 1.5 • remain.This is necessary because sinks along a flowline are incompatible with the forward dynamical model, which will fill them with ice and create undesirable spin-up issues.
The flowlines are then sorted according to their Strahler number (a measure of branching complexity defined by Strahler, 1952, and commonly used in hydrological applications), from the lowest (line without tributaries but with possible descendants) to the highest (the main -and longest -centerline).This ordering is important for the mass flow routing: indeed, each flowline contains a reference to its descendant, and this reference is used by the inversion and dynamical models to transfer mass from the tributaries towards the main flowline.
The width of each grid point along the flowline is computed in four steps.First, the catchment area of each flowline is computed using a routing algorithm similar to that used to compute the centerlines (Fig. 2a).Then the geometrical widths are computed by intersecting the flowline's normal to the boundaries of either the individual catchments or the glacier itself (Fig. 2b).These geometrical widths are then corrected by a factor specific for each altitudinal bin (Fig. 2c), so that the true altitude area distribution of the glacier is approximately preserved (Fig. 2d).Finally, these widths are multiplied by a single factor ensuring that the total area of the glacier is the exact same as the one provided by the RGI, ensuring consistency with future model intercomparisons.
At this stage, it is important to note that the map representation of the flowline glacier presented in Fig. 2c is purely artificial.
The fact that the glacier cross-sections are overlapping is irrelevant: the role of the flowlines is to represent the actual flow of ice as accurately as possible while conserving the fundamental aspects of the real glacier: slope, altitude, area, geometry.The flowline approximation is going to work better for valley glaciers (like Tasman glacier shown above) than for cirque glaciers (like the Upper Grindelwald).For ice-caps, the flowline representation is likely to work poorly, as discussed in Sect.3.6.2.
From Fig. 2c one can see that future improvements of the mass balance model based on e.g.topographical shading or snow redistribution are made possible by the knowledge about the flowlines' location.

Climate data and mass balance
The mass balance model implemented in OGGM is an extended version of the temperature index melt model presented by Marzeion et al. (2012b).The monthly mass balance m i at an elevation z is computed as: where P Solid i is the monthly solid precipitation, p f a global precipitation correction factor (defaults to 2.5, see Appendix A), µ * the glacier's temperature sensitivity, T i the monthly air temperature, T M elt the monthly air temperature above which ice melt is assumed to occur (default: -1 • C, chosen because melting days can occur even if the monthly average temperature is below 0 • C), and ε a residual (or bias correction) term.Solid precipitation is computed as a fraction of the total precipitation: 100% , and linearly interpolated in between.The parameter 5 µ * indicates the temperature sensitivity of the glacier and needs to be calibrated.For this paper, the temperature and precipitation time series  are obtained from gridded observations (CRU TS4.01; Harris et al., 2014, see Appendix A).The temperature lapse-rate is set to a constant value (default: 6.5 K km −1 ) or it can be time-dependant and computed from a linear fit of the 9 surrounding grid-points.
For the calibration of the temperature sensitivity parameter µ * we use the method described by Marzeion et al. (2012b) and successfully applied many times since then (e.g.Marzeion et al., 2014aMarzeion et al., , 2015Marzeion et al., , 2018)).Although the general procedure didn't change, its peculiarity justifies describing it here.We will start by noting that µ * depends on many factors, most of them being glacier-specific (e.g.avalanches, topographical shading, cloudiness), and others being related to systematic biases in the input data (e.g.climate, topography).As a result, µ * can vary greatly between neighbouring glaciers without obvious physical reasons.The calibration procedure implemented in OGGM makes use of these apparent handicaps by turning them into assets.
The procedure begins with glaciers for which we have direct observations of specific mass balance (N = 254, see Appendix B).For each of these glaciers, annual sensitivities µ(t) are computed from Eq. 1 by requiring that the glacier specific mass balance m(t) is equal to zero4 .m(t) is the glacier integrated mass balance computed for a 31-yr period centred around the year t and for a constant glacier geometry fixed at the RGI outline's date (e.g.2003 in the Alps).The process is illustrated in Fig. 3c (blue line): around 1920 the climate was cold and wet (Figs.3a and b), and as a consequence the hypothetical temperature sensitivity required to maintain the 2003 glacier geometry needs to be high.Inversely, the more recent climate is warmer and the temperature sensitivity needs to become smaller for the glacier to remain stable.
These hypothetical, time-dependent µ(t) are called "candidates": it is likely (but not certain) that at least one of them is the correct µ * .To determine which of the candidates is suitable, we then compute the mass balance time series for each of the µ(t) and compute their bias ε with respect to observations (red line in Fig. 3c).Note that the period over which the observations are taken is not relevant for the bias computation: each µ candidate can produce a mass balance for any year, as per Eq. 1.In comparison to observations, µ(t = 2000) is too low and produces mass balances with a positive bias.Inversely, µ(t = 1920) is too high and leads to a negative bias.For three years, the bias is close to or crossing the zero line and µ(t) is therefore very close to the ideal µ * .These dates represent the center of a 31-yr long climate period where today's glacier would be in equilibrium and maintain its current geometry.From these three candidates, we pick the date with the smallest bias and call it t * .This t * is an actual date but is mostly an abstract concept: we are going to make use of it for the next step.
For the vast majority of the glaciers, µ * and t * are unknown.For these we could interpolate the µ * (probably the most obvious solution), or we could interpolate t * : indeed, the procedure above can be reversed and t * can be used to retrieve µ * , again by requiring that m(t * ) is equal to zero (Eq.1).We interpolate t * to all glaciers without observations using inverse distance interpolation from the 10 closest locations (which can be quite far away, see Appendix B and C).The residual bias ε for glaciers with observations can be close to zero (the case for Hintereisferner in Fig. 3, where the bias curve crosses the zero line) but can also be higher (indicating that no 31-yr period in the last century would sustain the current glacier geometry).
When no perfect t * is found, the date with the smallest absolute bias is chosen.This residual ε is also interpolated between locations and added to the modelled mass balance.This residual may be significant at certain locations (up to 1.5 m yr −1 , median of 6 cm yr −1 ) and would benefit from further calibration e.g. with regional geodetic mass-balance estimates.The benefit of this approach is best shown by cross-validation (Fig. 4), where one can see that the error increases considerably when using µ * interpolation instead of the proposed method.This is due to several factors: the equilibrium constraint applied on µ(t) implies that the sensitivity cannot vary much during the last century.In fact, µ(t) at one glacier often varies less in one century than between neighbouring glaciers, because of the local driving factors mentioned earlier.In particular, it will vary comparatively little around a given year t: errors in t * (even large) will result in relatively small errors in µ * .
the equilibrium constraint will also imply that systematic biases in temperature and precipitation (no matter how large) will automatically be compensated by all µ(t), and therefore also by µ * .In that sense, the calibration procedure can be seen as an empirically driven downscaling strategy: if a glacier is located there, then the local climate (or the glacier temperature sensitivity) must allow a glacier to be there.For example, the effect of avalanches or a negative bias in precipitation input will have the same impact on calibration: the value of µ * should be lowered to take these effects into account, even though they are not resolved by the mass balance model.The most important drawback of this calibration method is that it assumes that two neighbouring glaciers should have a similar t * .This is not necessarily the case, as other factors than climate (such as the glacier size) will influence t * too.Our results (and the arguments listed above) show however that this is an approximation we can cope with.
In a final note, it is important to mention that µ * and t * should not be over-interpreted in terms of real temperature sensitivity or response time of the glacier.This procedure is primarily a calibration method, and as such it can be statistically scrutinized (for example with cross-validation).It can also be noted that the mass balance observations play a relatively minor role in the calibration: they could be entirely avoided by fixing a t * for all glaciers in a region (or even worldwide), without much performance loss.The observations, however, play a major role for the assessment of model uncertainty (Fig. 4).For more information about the climate data and the calibration procedure, refer to Appendix A.

Ice thickness
Measuring ice thickness is a labour-intensive and complex task, therefore only a fraction of the world's glaciers is monitored and direct measurements are sparse.A physical or statistical approach is necessary for modelling glacier evolution at the global scale.For a recent review of available techniques for ice thickness modelling, see Farinotti et al. (2017).OGGM implements a new ice thickness inversion procedure, physically consistent with the flowline representation of glaciers and taking advantage of the mass balance calibration procedure presented in the previous section.It is a mass-conservation approach largely inspired by Farinotti et al. (2009), but with distinct characteristics.
The principle is quite simple.The flux of ice q [m 3 s −1 ] through a glacier flux-gate (cross-section) of area S [m 2 ] reads: with u the average velocity [m s −1 ].Using an estimate for u and q obtained from the physics of ice flow and the mass balance field, S and the local ice thickness h [m] can be computed relying on some assumptions about the bed geometry.We compute the depth-integrated ice velocity using the well known shallow-ice approximation (Hutter, 1981(Hutter, , 1983)): with A the ice creep parameter [s −1 Pa −3 ], n the exponent of Glen's flow law (n=3), and τ the basal shear stress, computed as: with ρ the ice density (900 kg m −3 ), g the gravitational acceleration (9.81 m s −2 ) and α the surface slope computed numerically along the flowline.Optionally, a sliding velocity u s can be added to the deformation velocity to account for basal sliding.We use the same parameterisation as Oerlemans (1997), who relied on Budd et al. (1979): with f s a sliding parameter (default: 5.7 × 10 −20 s −1 Pa −3 ).If we consider a point on the flowline and the catchment area Ω upstream of this point, mass conservation implies: with ṁ the mass balance [kg m −2 s −1 ], and m = ṁ − ρ∂h/∂t the "apparent mass balance" after Farinotti et al. (2009).If the glacier is in steady state, the apparent mass balance is equivalent to the actual (and observable) mass balance.In the non-steady state case, ∂h/∂t is unknown, and neither is the time integrated (and delayed) mass balance Ω ṁ responsible for the flux of ice through a section of the glacier at a certain time.Farinotti et al. (2009) and Huss and Farinotti (2012) deal with the issue by prescribing an apparent mass balance profile as a parameterized linear gradient which is, arguably, more a semantic than a physical way to deal with the transience of the problem.
Like Huss and Farinotti (2012), OGGM cannot deal with the transient problem yet: we deliberately assume steady state and therefore set m = ṁ.This has the strong advantage that we can make direct use of the equilibrium mass balance m(t * ) computed earlier, which satisfies m = 0 by construction.q is then obtained by integrating the equilibrium mass balance m along the flowline(s).The tributaries will have a positive flux at their last grid point: this mass surplus is then transferred to the downstream line, normally distributed around the 9 grid points centred at the flowlines' junction.By construction, q starts at zero and increases along the major flowline, reaches its maximum at the equilibrium line altitude (ELA) and decreases towards zero at the tongue (for glaciers without frontal ablation).
Equation 2 turns out to be a polynomial of degree 5 in h with only one root in R + , easily computable for each grid point.
Singularities due to flat areas are avoided since the constructed flowlines are not allowed to have a local slope α below a certain threshold (default: 1.5  m −1 , ELA altitude of 2600 m a.s.l.).For (d), the glacier is first grown to equilibrium then shrunk for 60 years after an ELA shift of +200 m. the glacier width) or rectangular (S = hw) bed shape.The default in OGGM is to use a parabolic bed shape, unless the section touches a neighbouring catchment (see Fig. 2c), neighbouring glacier (ice divides, computed from the RGI), or at the terminus of a calving glacier.In these cases the bed shape is rectangular.Optionally, OGGM can also compute the effect of lateral bed stresses (Cuffey and Paterson, 2010) following a parameterization and tabular correction factors developed by Adhikari and Marshall (2012b).
Figure 5 displays some examples taken from the OGGM test suite, where the automated inversion procedure is applied on idealized glaciers generated with OGGM's flowline model (see Sect. 3.5).In the equilibrium cases (Fig. 5a to c), the inverted topography is nearly perfect.Differences arise at strong surface gradients, mostly because of numerical differences (the inversion method uses a second order central difference which tends to smooth the slope).The transient case (Fig. 5d) illustrates the consequences of the steady-state assumption: although the glacier is retreating, the constraint m = 0 leads to 10 a lowered ELA and, even with a perfectly known mass balance gradient, leads to an overestimated ice thickness (in this case, 25%).This effect is visible everywhere, but is strongest at the tongue.The importance of the steady state assumption on ice thickness estimates has been studied with numerical experiments (e.g.Adhikari and Marshall, 2012a) and is often compensated by calibration in real world applications.Bahr et al., 1997Bahr et al., , 2015) ) and based on point observations (Fischer and Kuhn, 2013).For (a), additional sensitivities are computed with an additional sliding velocity (Oerlemans, 1997) and his sliding parameter fs.For (b), additional sensitivities are computed with a varying creep parameter A.
The sensitivity of the inversion procedure to various parameters is illustrated with the example of the Hintereisferner glacier (Fig. 6).The total volume (and the local thickness) is very sensitive to the choice of the creep parameter A, varied from a factor 1/10 to 10 times the default value of 2.4 × 10 −24 s −1 Pa −3 (Cuffey and Paterson, 2010).With a smaller A, the ice is stiffer and the glacier gets thicker (A is expected to get smaller by one or more orders of magnitude with colder ice temperatures).
Inversely, softer ice leads to a thinner glacier.The shape of the curve is proportional to the fifth root of the fraction 1/A, explaining why the volume gets very sensitive to small values of A. Adding sliding reduces the original thickness significantly for the same reasons as an increasing A, since both sliding and ice rheology (A) have a strong influence on the computed ice flux q.Inversely, adding lateral bed stresses reduces ice velocity and increases the computed ice volume.Changing from a rectangular to a parabolic bed shape yields a volume loss of exactly 1/3: this is expected from geometrical considerations.The mixed parabolic/rectangular bed shape model implemented by default therefore lies in between.
The total precipitation amount, by acting on the mass balance gradient and therefore on the ice flux q will also play a nonnegligible role for the ice thickness (Fig. 6b).The effect is small in comparison to the influence of A, but it is noticeable: glaciers located in maritime climates (with high values of accumulation) will be thicker on average than similar glaciers in drier conditions.
This example shows that one can always find an optimum (and non-unique) set of parameters leading to the correct total volume.In practice, however, calibrating the model for accurate global glacier volume estimates is a major challenge for global glaciological models and will be the topic of a separate study.The IACS Working Group on Glacier Ice Thickness Estimation 5 is working towards this goal: OGGM participated in the first Ice Thickness Models Intercomparison eXperiment (ITMIX, Farinotti et al., 2017), ranking amongst the best models with limited data requirements.

Ice dynamics
At this stage of the processing workflow, the ice-dynamics module is straightforward to implement.Provided with the mass balance, slope, width w and bed topography along the flowline, we solve the equation: numerically with a forward finite difference approximation scheme on a staggered grid.Numerical stability is ensured by the use of an adaptive time stepping scheme following the Courant-Friedrichs-Lewy (CFL) condition ∆t = γ ∆x max(u) with γ as the dimensionless Courant number chosen between zero and one.Unlike many solvers of the shallow-ice equation, we do not transform Eq. 7 to become a diffusivity equation in h, but solve it as it is formulated here.This has the advantage that the numerical solver is the same regardless of the shape of the bed (parabolic, trapezoidal or rectangular).The new section S at time t + ∆t allows to compute h(t + ∆t) according to the local bed geometry.Therefore, it is possible to have changing bed geometries along a single flowline using the same numerical solver.The drawback of our approach is that we cannot take advantage of the diffusivity equation solvers already available elsewhere.We tested our solution against the robust and and mass-conservative solver presented by Jarosch et al. (2013): our model yields accurate (and faster) results in most cases, but fails to ensure mass-conservation for very steep slopes like most other solvers to date.While a flowline version of the solver presented by Jarosch et al. (2013) is available in OGGM, it is not used operationally since it cannot yet handle varying bed shapes and multiple flowlines -it will become the default solver when these elements are implemented.
At a junction between a tributary and its downstream line, an artificial grid point is added to the tributary line.This grid point has the same section area S and thickness h as the previous one, but the surface slope is computed from the difference in elevation between tributary and descendant flowline.This is necessary to ensure a dynamical connection between the two lines: when the main flowline is at a higher elevation than its tributary, no mass exchange occurs and the tributary will build up mass until enough ice is available.At a junction point, Eq. 7 therefore contains an additional mass flux term from the tributary.
Before the actual run, a final task merges the output of all preprocessing steps and initialises the flowline glacier for the model.For the glaciers to be allowed to grow, a downstream flowline is computed using a least cost routing algorithm leading the glacier towards the domain boundaries (this algorithm is similar to the algorithm used to compute the glacier centerlines).
The bed geometries along the downstream line are computed by fitting a parabola to the actual topography profile.In case of bad fit, the values are interpolated or a default parabola is used.Along the glacier, where the bed geometries are unknown before the inversion, the bed geometries are either rectangular (ice divides and junctions) or parabolic.Very flat parabolic shapes can happen occasionally, for wide sections with a shallow ice thickness.These geometries are unrealistically sensitive to changes in h.They create a strong positive feedback (the thickening of ice leading to a highly widening glacier) and are b:  c:   therefore prevented: when the parabola parameter falls below a certain threshold, the geometry is assumed to be trapezoidal instead.
The coupling between the mass-balance and ice dynamics modules is a user choice.The spatially distributed mass-balance used by the dynamical model can be updated: (i) at each time-step of the dynamical model's computation, (ii) each month, (iii) each mass-balance year (the default), or (iv) only once (for testing and feedback sensitivity investigations).In practice, this does not make much difference for the yearly averages of glacier change (except for option iv), and the choice of a yearly update is mostly driven by performance considerations.Note that the mass-balance model can compute the mass-balance at shorter time intervals if required by the physical parameterizations: the interface between the model elements simply requires the mass-balance model to integrate the mass-balance over a year before giving it to the dynamical model.
The results of two idealized simulations with an advancing and a shrinking scenario are shown in Fig. 7.When put under the cold and wet climate of the beginning of the 20th century, Hintereisferner would grow about 2/3 larger than it is today.
Inversely, the glacier is in strong disequilibrium with today's climate: it would lose about 2/3 of it's volume if the climate remained as it was during the past 31 years.The response time of the glacier is approximately twice as fast in the shrinking case, and the natural random variations of the glacier are much smaller than for a large glacier with more inertia and a longer response time.The previous results were obtained with the default set-up of OGGM.In Fig. 8 we assess the sensitivity of the dynamical model to changes in the creep parameters A and to the addition of lateral drag and basal sliding velocity.As expected, these dynamical parameters affect the equilibrium volume and the response time of the glacier (faster ice leading to a thinner glacier, and inversely).Because of the mass-balance elevation feedback, the stiffer and therefore thicker glacier is also larger and longer, but its response to climate variability is smaller in amplitude than that of the fast moving sliding glacier.
A and f s depend on many factors such as ice temperature or basal characteristics and they cannot be assumed to be globally constant.They are considered as calibration parameters in OGGM, and will be tuned towards observations of ice thickness or glacier length changes.In this study we only calibrate the mass-balance model while the ice dynamics parameters are set to their default values (A = 2.4 × 10 −24 s −1 Pa −3 , f s = 0, no lateral drag).Nevertheless, we discuss the model sensitivity to these dynamical parameters for individual glaciers (Fig. 8) or global runs (Fig. 10).

Special cases and model limitations
The previous experiments demonstrate that the OGGM model is capable of simulating the dynamics of glaciers in a fully automated manner.In this section we describe the implications of the flowline approximation in the special cases of waterterminating glaciers and ice caps, and discuss some examples of glaciers with a less trivial geometry.

Water-terminating glaciers
Glaciers are defined as "water-terminating" in OGGM when their RGI terminus attribute is either flagged as marine-terminating or lake-terminating.The major difference between a water-terminating glacier and a valley glacier is the additional mass-loss that occurs at the glacier front (frontal ablation).This has implications for the bed thickness inversion which currently assumes that the mass-flux at the front is zero (by setting m = 0, see Sect.3.4) and for the dynamics of the glacier.The current treatment of water-terminating glaciers in OGGM is very simple but explicit: we do not take frontal ablation into account for the bed inversion (i.e. the original glacier front has a thickness of zero), but we do have a basic parametrisation in the ice dynamics module.We add a grid point behind the glacier front which is reset to zero ice thickness at each time step: the ice mass suppressed this way is the frontal ablation flux, that we store.This parametrisation has the advantage of preventing waterterminating glaciers from advancing while still allowing them to retreat (in which case they stop calving).We are working on a more advanced frontal ablation parametrisation for both the ice dynamics and the ice thickness inversion (Recinos et al., 2018).

Ice caps and ice fields
Ice caps and ice fields in the RGI are divided in single dynamical entities separated by their ice divide (Fig. 9).However, the entities that belong to an ice cap are classified as such in the RGI: currently, the only special treatment for these entities in OGGM is that only one major flowline is computed (without tributaries).Indeed, the geometry of ice caps is often non trivial and it is not clear whether tributaries would really improve the model results.An example of an ice cap is shown in Fig. 9: while the general behaviour of the ice cap is reasonably simulated by the flowline model (e.g. at the outlet glaciers), other features appear to be unrealistic (e.g.close to the ice-divides).Moreover, the mass-conservation inversion method is probably underestimating the real ice-thickness at the location of the ice-divide, where other processes related to the past history of the ice cap are at play.A possible way forward would be to run a distributed shallow-ice model instead of the flowline representation, and it is part of our long-terms plans to do so.

Glacier complexes
Single glaciers can be defined as the smallest dynamically independent entity, i.e. the boundaries between two glaciers should approximately follow the ice divides or hydrological basin boundaries.The flowline assumption strongly relies on this condition being true, and indeed most of the RGI glaciers are properly outlined.Unfortunately there are notable exceptions, for three main reasons: human decision: some well known glaciers have historical boundaries that the inventory provider wanted to keep, although the glacier is now disintegrated in smaller entities.A good example is the Hintereisferner glacier (Fig. 7), which should have three outlines instead of one.
uncertainties in the topography: the inventories are often generated using both automated processes and manual editing.
There is no guarantee that we use the same DEM as the original inventory, and therefore OGGM and RGI might disagree on the ideal position of an ice divide.
unavailable data: some remote glaciers and ice caps are outlined in the RGI, but not divided at all.These are the most problematic cases, and should be a matter of concern for all RGI users.For example, the largest glacier in RGI (an ice  cap in north-eastern Greenland with id RGI60-05.10315and area 7537 km 2 ) is wrongly outlined and should be separated in at least a dozen of smaller entities.
Most of the small errors are filtered out by OGGM with algorithms based on surface slope thresholds (see Sect. 3.2), but the latter group of glaciers should be handled upstream.We have developed an open source tool to automatically compute glacier divides (https://github.com/OGGM/partitioning,based on Kienholz et al., 2013), but do not use it here.This issue is a large 5 source of uncertainty for ice thickness estimates and dynamical modelling of glaciers in general, and could be the subject of a dedicated study.

Glacier centric modelling
Like most global glacier models, OGGM simulates each glacier individually.This has evident practical advantages, and is also a strong asset for our mass balance model calibration algorithm.However, this has two major drawbacks: (i) neighbouring 10 glaciers won't merge although they grow together, and (ii) we can only simulate glaciers which are already inventoried, while uncharted glaciers are simply ignored.Both errors are a source of uncertainty for long or past simulations but less for short term projections in a warming world.The most obvious way to deal with this issue is to use distributed models (e.g.Clarke et al., 2015), with their own drawbacks (mostly: computational costs and the need for distributed mass balance fields).Another way would be to allow the dynamical merging of neighbour flowline glaciers at run time.While both are viable options for the OGGM workflow, they represent a considerable increase in complexity and are not available yet.Like other fundamental issues described in this paper (such as missing topographical data or wrongly outlined glaciers), this problem will affect other glacier models as well.We hope that some of the tools we introduce here will help to solve some of these issues upstream, and that the community will soon be able to put pressure on commercial data providers for better data availability.

Global simulations
Thanks to its automated workflow, OGGM is able to apply all processes described in the previous section to all glaciers of the globe with the exception of Antarctica, where no CRU data is available (see Appendix C for an overview of the RGI regions).
No special model setup is needed, we use all model default settings without any calibration (this is not strictly true for the µ * calibration, which is an automated process and cannot be tuned or turned off).In the following analyses the focus is put on the model behaviour and not on the quantitative results.However, as we are going to see our results are close to expectations even without calibration, indicating a realistic model behaviour.

Hardware requirements and performance
Thanks to the computational efficiency of the flowline model, OGGM runs quickly enough to be used on a personal computer for up to a hundred glaciers.At the global scale a high performance computing environment is required.For these global simulations we used a small-sized cluster comprising two nodes with 16 quad-core processors each, resulting in 128 parallel threads.With this configuration, the model preprocessing chain (including the ice thickness inversion) takes about 7 hours to complete (without data download).The total size of the (compressed) preprocessed output is 122G, an amount which can be reduced by deleting intermediate computing steps.The amount of required storage increases with each dynamical run; here again it is possible to reduce the data amounts by storing only diagnostic variables such as volume, area, length, ELA instead of the full model output.The dynamical runs are the most expensive computations: running five 300-yr long global runs takes about 24 hours on our small cluster, a very satisfying performance.It is interesting to note that because of the adaptive time-step, glacier shrinkage scenarios run faster than growing ones.

Invalid glaciers
Due to uncertainties in the input data (topography, outlines, climate), a certain number of glaciers fail to be modelled by OGGM.The statistics of these invalid glaciers are summarized in Table 1.The largest amount of errors (2.6% of the total area) is due to invalid climate series.Errors occur mostly when the climate is too cold for melt to happen or, inversely, too warm or too dry for accumulation to happen.While some of these errors are directly due to incorrect climate data, some can also be attributed to missing processes in the OGGM mass balance model: sublimation and frontal ablation, both leading to mass-loss even at cold temperatures.The least problematic source of error (0.2% of the total area) is due to failures during the actual dynamical run.The large majority of dynamical failures (751 out of 772) happen because the glacier exceeded the domain boundaries at run time.Some of these errors could be mitigated by increasing the domain size (at the cost of computational efficiency).Only 21 glaciers fail due to numerical instabilities.Finally, there is a number of other errors (0.3% of the total area) happening at other stages of the model chain.Examples include errors in processing the geometries or failures in computing certain topographical properties due to invalid DEMs.Altogether, 7084 glaciers (3.1% of the total area) cannot be modelled by OGGM.There are strong regional differences, remote low and high latitude regions accounting for most of the errors.

Volume inversion
A summary of the volume inversion results is presented in Fig. 10.As expected from theory (Bahr et al., 1997(Bahr et al., , 2015)), our glacier volume estimates approximately follow a power law relationship with the glacier area (V = c S γ ).The coefficients obtained by a linear fit in log space are close, but not equal to the coefficients computed by Bahr et al. (1997).In particular, the OGGM fit is slightly flatter than the theoretical value (Fig. 10, upper panel), in accordance with empirical coefficients (e.g.Bahr et al., 2015;Grinsted, 2013).This is an encouraging result, especially because it was reached with the OGGM default settings and without calibration.
The global volume estimates are particularly sensitive to the choice of the ice dynamics parameters, as shown in Fig. 10 (lower panel).As for individual glaciers, the total volume follows an inverse polynomial curve as expected from the equations of ice flow.Changing from a rectangular to a parabolic bed shape yields a volume loss of exactly 1/3 (see Sect. 3.4).Adding lateral drag yields a volume very close to the rectangular case: although this is fortuitous (individual glaciers can show different results, see Fig. 6), it matches nicely the original purpose of the parametrization, which is to compute a more realistic ice flow for parabolic bed shapes.The three independent estimates plotted as straight dotted lines (VAS; Huss and Farinotti, 2012;Grinsted, 2013) illustrate that A is a relatively straightforward parameter to act upon in order to fit the model to observations.The effect of A, however, is going to be the same on all glaciers and therefore will be a poor measure of performance (see also Bahr et al., 2015, Sect. 8.11).In fact, the added value of OGGM is more likely to be found in the deviations from the scaling law (Fig. 10, lower panel).The deviations are the result of a range of possible factors such as slope, total accumulation, or altitude area distribution.With accurate boundary conditions, OGGM should be able to provide more accurate estimates, within the limits of the assumptions and simplifications behind the model equations.The calibration and validation of the OGGM inversion model will be the topic of a subsequent study.

Dynamical runs
We test the model behaviour by running several 300-yr long global simulations under various climate "scenarios".In the first simulations (Fig. 11), we run the model under the climate of the past 31 years.In order to keep the forcing realistic, we create a pseudo-random climate by shuffling the years infinitely.We also run two additional simulations with a 0.5 • C positive and negative bias.The unbiased simulation illustrates the committed glacier mass-loss, i.e. the ice mass which is not sustainable under the current climate.Figure 11 shows that all regions will continue to lose ice even if the climate remains constant.The regions with the largest committed mass loss relative to the initial volume are Western Canada and US (02), Svalbard (07), and the three High Mountain Asia regions (13,14,15).Inversely, the regions Arctic Canada South (04), Greenland (05) and Iceland (06) are least affected.The reasons for these regional differences are complex: they are due to the climate itself of course, but also to glacier properties such as size, slope, and continentality.The regions that are far from equilibrium also tend to be less sensitive to the temperature bias experiments, although this should not be over-interpreted (indeed, the range of the y-axes can hide differences which appear small in comparison to the large regional glacier loss).In general, the model behaviour looks reasonable and the regional differences are in qualitative agreement with other global studies (e.g.Huss and Hock, 2015, where the regions with a stronger response to 21st century climate change are the same as the ones listed above).Also our global estimate of the committed mass loss (approx.33% at the end of the 300-yr simulation, probably more at equilibrium) is in agreement with other studies (27±5%, 38±16%, and 36±8% for Bahr et al., 2009;Mernild et al., 2013;Marzeion et al., 2018, respectively).
A further model test is presented in Fig. 12.Here, we apply a new climate scenario: the climate at t * which, for each glacier individually, represents a theoretical equilibrium climate.In addition to the global response to these scenarios, we separate between the majority group of smaller glaciers and the much smaller group of very large glaciers.Both groups are selected so that they sum up to one quarter of the total glacier volume.A striking feature of the runs is that the glaciers tend to grow under the artificial t * climate: the growth is slow at first and accelerates with time, hinting towards a positive feedback.This feedback is driven by two factors: first, a higher surface elevation leads to a positive change in mass balance (mass balance-elevation feedback); second, because of the parabolic and trapezoidal bed shapes, a larger ice thickness leads to a wider accumulation area above the ELA and to a wider ablation area below the ELA.It appears that the positive width-accumulation feedback is stronger than the negative width-ablation feedback.This can be explained by the larger accumulation area of glaciers in an equilibrium climate: the average accumulation area ratio at t * in OGGM is 51%.In order to test which of these feedbacks is stronger, we run a simulation with rectangular bed shapes exclusively (dotted light purple line in Fig. 12), therefore eliminating the width-accumulation but keeping the mass balance-elevation feedback.The results show that for the vast majority of glaciers the feedback almost disappears, while the very large glaciers still show a weak and delayed altitude feedback.
It is unclear whether this is a bug or a feature.On the one hand, this behaviour is not really desirable since one would expect glaciers to remain constant under a theoretical equilibrium climate.On the other hand, t * is just a vehicle to calibrate the model and was not supposed to yield a particular insight (for example, many glaciers can only have an equilibrium t * climate after the application of a bias to the operational mass balance model).There are many reasons why small initial perturbations such as numerical noise or the differences between the bed inversion and forward model numerical schemes might lead to a different equilibrium.It must also be noted that this feedback is slow to appear, and will only have a notable influence on the largest glaciers for long term simulations in a cooler climate (the global volume change after 100 years due to the feedback is 2.4% for the default and 1% for the all rectangular cases).The very simple definition of an "equilibrium climate" for these very large glaciers is problematic anyway: large glaciers have a very slow but potentially large response to the smallest changes in climate.At the global scale, most of the 300-yr volume loss is due to the small glaciers, which respond faster and stronger than larger ones.

Conclusions
We present a new model of global glacier evolution, the Open Global Glacier Model (OGGM, v1.1).The panoply of tools available to compute past and future glacier change range from simple box models (e.g.Harrison, 2013) to more complex, geometry aware models (Huss and Hock, 2015, to cite the most recent in date).OGGM undoubtedly belongs to the complex  Note the units of the y-axes (10 3 km 3 ) and the marked regional differences.side of this scale.Different model complexities are justified by different problem settings, taking into account the modelspecific merits and drawbacks.Instead of endorsing one approach over the other, OGGM aims to provide a framework which allows one to switch between models and allows objective intercomparisons.In fact, the ice dynamics module represents only a small fraction of the OGGM codebase: a huge amount of work has been invested to provide a series of tools which will help others in their own modelling endeavours.Any interested person can download, install, and run these tools at no cost.This includes the automated download of topographic and climate data for any location on the globe, the collation of glacier attributes, the automated computation of glacier centerlines, or the delineation of glacier dynamical entities.While some of these tools have been described elsewhere, the added value of OGGM is that they are now centralised, documented, and available for public review via the open-source model.
In the future, we will continue to encourage external contributions in several ways.First, it must be as easy as possible for a new user to detect where and how a contribution can be implemented: here, documentation is key.Then, the model must be able to cope with different ways to simulate a considered process: every single task in the OGGM workflow can be replaced or enhanced, as long as the format of the input and output files is agreed beforehand.Perfect modularity will be hard to achieve, but the recent implementation of alternative numerical solvers show that modularity is possible.Finally, we need to ensure attribution to the original contribution (e.g. a scientific publication) in order to engage the wider community.For this purpose, we developed a template repository for external OGGM modules: https://github.com/OGGM/oggmcontrib.This development model will ensure that users importing OGGM extensions will be aware of the source of each module they are using and will be able to refer to the original contribution appropriately.We hope that this development model will foster new collaborations.
We cannot (and do not want to) demonstrate that OGGM will provide more accurate estimates of future sea-level rise than earlier attempts.However, OGGM allows new studies which weren't possible before.The dynamical representation of glacier advances and retreat enables studies of glacier evolution at long (paleo-) time scales, where ice dynamics and geometrical attributes such as the accumulation area ratio play an important role (e.g.Mackintosh et al., 2017).First OGGM simulations over the last millennium show very promising results (Goosse et al., 2018).The modular framework allows one to compare the performance of various parametrisations such as the mass balance and downscaling algorithms.It may be argued that the amount of available data is not enough to constrain modelling studies such as ours at the global scale: OGGM can now be used to test this argument by allowing simpler modules to be added to the codebase and test the added value of increased complexity.
Planned and envisioned future developments for the model follow the general guidelines of modularity and extendability.
While some of the authors are working on adding even more complexity to the model (for example by improving the frontal ablation and mass balance parametrisations or by implementing a distributed ice dynamics module), it is part of our plans to implement simpler approaches such as the original Marzeion et al. (2012b) model or the Huss and Farinotti (2012) approach to ice thickness estimation.A considerable amount of work will be needed to correctly assess the uncertainties associated with the model chain: here, Monte Carlo and Bayesian approaches might be the way to follow.
The non-linear dynamical behaviour of glaciers raises a wide range of very interesting inverse problems.For example, how to deal with the transient climate issue in the ice thickness inversion algorithm?How much information about past climate can be extracted from moraine proxies and today's glacier extent?What are the uncertainties associated with global sea-level rise estimates, and where do they originate?How much complexity is just right?These are all questions the authors hope will be easier to address through the publication of OGGM.
6 Code availability, testing, and software requirements The OGGM software is coded in the Python language and licensed under the GPLV3 free software license.The latest version of the code is available on Github (https://github.com/OGGM/oggm), the documentation is hosted on ReadTheDocs (http: //docs.oggm.org),and the project webpage for communication and dissemination can be found at http://oggm.org.Past and intermediate versions are available in a permanent DOI repository (https://zenodo.org/badge/latestdoi/43965645).The software ships with an extensive test suite which can be used by the users to test their configuration.The tests are triggered automatically at each new code addition, reducing the risk of introducing new bugs (https://travis-ci.org/OGGM/oggm).The suite contains unit tests (for example for the numerical core) and integration tests based on sets of real glaciers.At the time of writing, 85% of all relevant lines of code are covered by the tests (i.e.called at least once by the test suite).The remaining 15% are challenging to monitor because they mostly concern the automated downloading tools which are used in production and cannot be tested automatically.

Figure 1 .
Figure 1.Example of the OGGM workflow applied to Tasman glacier, New Zealand; a: topographical data preprocessing; b: computation of the flowlines; c: geometrical glacier widths determination (the colors indicate the different flowlines); d: width correction according to catchment areas and altitude-area distribution (see Fig. 2 and main text for details); e: ice thickness inversion; f: random 100-yr long glacier evolution run leading to a glacier advance.See Sect.2.1 for details.

Figure 2 .
Figure 2. Example of the flowlines' width determination algorithm applied to the Upper Grindelwald glacier, Switzerland; a: determination of each flowline's catchment area; b: geometrical widths; c: widths corrected for the altitude-area distribution, the bold lines representing the grid points where the cross-section touches a neighbouring catchment; d: frequency distribution of the glacier area per altitude bin, as represented by OGGM and by the SRTM topography.

Figure 3 .
Figure 3. Calibration procedure for µ * applied to the Hintereisferner glacier, Austria.a and b: annual and 31-yr average of temperature and precipitation obtained from the nearest CRU grid point (altitude 2700 m a.s.l.).c: time series of the µ candidates (mm yr −1 K −1 ) and their associated mass balance bias (mm w.e.yr −1 , right-axis) in comparison to observations.The vertical dashed line marks the time where the bias is closest to zero (t * ).

Figure 4 .
Figure 4. Benefit of spatially interpolating t * instead of µ * as shown by leave-one-out cross-validation (N = 254).Left: error distribution of the computed mass balance if determined by the interpolated t * .Right: error distribution of the mass balance if determined by interpolation of µ * .The vertical lines indicate the mean, median, 5% and 95% percentiles.See https://cluster.klima.uni-bremen.de/~github/crossvalfor an online visualisation of model performance for each glacier.

Figure 5 .
Figure 5. Idealized inversion experiments: we compute the bed topography out of the surface elevation obtained from a flowline model applied to a predefined bed topography.(a), (b) and (c): glacier grown to equilibrium with different bed topographies (flat, cliff, random).(d): transient experiment with a shrinking glacier.The same mass balance profile is used for all experiments (linear gradient of 3 mm w.e.

Figure 6 .
Figure 6.Total volume of the Hintereisferner glacier computed with (a) varying factors for the default creep parameter A, and (b) varying precipitation factors.The dotted and dashed black lines display the total volume estimated with volume-area scaling (VAS,Bahr et al., 1997,

Figure 7 .
Figure 7. Evolution of the Hintereisferner under two random forcing scenarios and for the default parameter set.For each scenario, the "climate years" during a 31-yr period are shuffled randomly, therefore creating a realistic climate representative for a given period.(a): the glacier volume evolution for each scenario (the black line marks the initial computed glacier volume).(b) and (c): the glacier shape at the end of the 800 yrs simulation for each case.

Figure 8 .
Figure 8. Evolution of volume (a) and length (b) of the Hintereisferner glacier under a random climate forcing generated by shuffling the "climate years" representative for the 31-yr period centred on t * .The glacier is reset to zero for each simulation, and the bed topography is obtained with the default parameters.The sensitivity to the addition of a sliding velocity or to a halving of the creep parameter A are also shown.The noisy patterns of the length time series are due to the fact that the length of a glacier on a discrete grid is sensitive to small interannual variations.

Figure 9 .
Figure 9. OGGM inversion workflow applied to the RGI entities of the Eyjafjallajökull ice cap, Iceland.Upper panel: outlines and topography.Lower panel: glacier thickness.

Figure 10 .
Figure 10.Global glacier volume modelling.Upper panel: binned scatter plot of volume versus area for all valid glaciers (N=207,438) with the default OGGM setup.Color shading indicates the number of glaciers in each bin.Note the logarithmic scale of the axes and the irregular color scale levels.The dashed lines indicates the volume area scaling relationship with either the theoretical parameters from Bahr et al. (1997) (V = 0.034 S 1.375 ) or fitted on our own data (V = 0.042 S 1.313 ).Lower panel: global volume estimates as a function of the multiplication factor applied to the ice creep parameter A, with five different set-ups: defaults, with sliding velocity, with lateral drag, and with rectangular and parabolic bed shapes only (instead of the default mixed parabolic/rectangular).In addition, we plotted the estimates from standard volume area scaling (VAS, V = 0.034 S 1.375 ), Huss and Farinotti (2012) (HF2012) and Grinsted (2013) (G2013).The two latter estimates are provided for indication only since they are based on a different glacier inventory.

Figure 12
Figure 12. a: global glacier volume change under various climate scenarios (1985-2015 climate with three temperature biases and climate at t * which, for each glacier individually, represents a theoretical equilibrium climate) and model configurations (rectangular bed instead of the mixed default), plotted as a fraction of the initial volume.b and c: volume changes of all glaciers making up for the first and last quartile of the sorted cumulative total volume.

Figure C1 .
Figure C1.Upper panel: map of the RGI regions; the red dots indicate the glacier locations and the blue circles the location of the 254reference WGMS glaciers used by the OGGM calibration.Lower panel: region names and basic statistics of the database (number of glaciers per region, regional contribution to the global area in percent, and percentage of the regional area which cannot be modelled by OGGM).

Table 1 .
Statistics of the model errors for each RGI region.The column names indicate which processing step produces an error, the value is the number of invalid glaciers and (in parentheses) the percentage of regional area they represent.