Robust uncertainty assessment of the spatio-temporal transferability of glacier mass and energy balance models

Energy and mass-balance modelling of glaciers is a key tool for climate impact studies of future glacier behaviour. By incorporating many of the physical processes responsible for surface accumulation and ablation, they offer more insight than simpler statistical models and are believed to suffer less from problems of stationarity when applied under changing climate conditions. However, this view is challenged by the widespread use of parameterizations for some physical processes which introduces a statistical calibration step. We argue that the reported uncertainty in modelled mass balance (and associated energy flux components) are likely to be understated in modelling studies that do not use spatio-temporal cross-validation and use a single performance measure for model optimization. To demonstrate the importance of these principles, we present a rigorous sensitivity and uncertainty assessment workflow applied to a modelling study of two glaciers in the European Alps, extending classical best guess approaches. The procedure begins with a reduction of the model parameter space using a global sensitivity assessment that identifies the parameters to which the model responds most sensitively. We find that the model sensitivity to individual parameters varies considerably in space and time, indicating that a single stated model sensitivity value is unlikely to be realistic. The model is most sensitive to parameters related to snow albedo and vertical gradients of the meteorological forcing data. We then apply a Monte Carlo multi-objective optimization based on three performance measures: model bias and mean absolute deviation in the upper and lower glacier parts, with glaciological mass balance data measured at individual stake locations used as reference. This procedure generates an ensemble of optimal parameter solutions which are equally valid. The range of parameters associated with these ensemble members are used to estimate the cross-validated uncertainty of the model output and computed energy components. The parameter values for the optimal solutions vary widely, and considering longer calibration periods does not systematically result in better constrained parameter choices. The resulting mass balance uncertainties reach up to 1300 kg m⁻², with the spatial and temporal transfer errors having the same order of magnitude. The uncertainty of surface energy flux components over the ensemble at the point scale reached up to 50 % of the computed flux. The largest absolute uncertainties originate from the short-wave radiation and the albedo parameterizations, followed by the turbulent fluxes. Our study highlights the need for due caution and realistic error quantification when applying such models to regional glacier modelling efforts, or for projections of glacier mass balance in climate settings that are substantially different from the conditions in which the model was optimized.

1 Introduction

Surface energy and mass balance models are valuable tools for estimating the response of glaciers to meteorological forcing (Oerlemans, 2011). Model results can be used to estimate regional run-off and resultant sea level rise (e.g. Hock, 2005), but additionally, and unlike results of empirical melt models, they can also be used to characterize the fundamental processes and key drivers of melt on glaciers, which is critical for understanding how they may behave under the influence of changing climate (e.g. Mölg and Hardy, 2004; Klok and Oerlemans, 2004; Hock and Holmgren, 2005; Mölg et al., 2008; Prinz et al., 2016; Willeit and Ganopolski, 2017).

All glacier surface mass and energy balance models contain a degree of parameterization of physical relationships. These parameters are either optimized to fit observations, or chosen based on previously established empirical relationships, or are a mix thereof. Uncertainty surrounding the transferability of parameterizations in both space and time poses a critical limitation on the usefulness of such models for regional upscaling of glacier behaviour or forward projections of global glacier behaviour under changing climate conditions.

Early energy balance studies typically apply models at a single point in space for which local physical relations can be readily established empirically, or direct measurements are available to tune the parameterizations (e.g. Mölg and Hardy, 2004; Greuell and Konzelmann, 1994; Bintanja and Van Den Broeke, 1995). Optimizing a model to local measurements can successfully reproduce local melt rates or surface temperature (e.g. Oerlemans and Knapp, 1998), and where this is the case, reliable simulation of glacier ablation is often taken to mean that the model also accurately reveals the relative importance of specific energy sources to ice ablation. Model optimization based on data from a single site, or from a very short time series, is, however, prone to parameter over-fitting, meaning that parameters are specifically adjusted to the study location and/or time (Beven, 1989). This can be evident in upscaling point optimizations to the glacier scale: for example, Klok and Oerlemans (2002) applied a distributed energy balance model to a mid-latitude glacier, using a combination of previously published parameter values and values estimated from local point-scale measurements, and found reasonable agreement for local energy fluxes. The albedo parameterization was identified as a potential source of uncertainty as it was based on data from a single point and 1 year of observations (Klok and Oerlemans, 2002; Oerlemans and Knapp, 1998) and may not be valid elsewhere on the glacier surface throughout all seasons (Van De Wal et al., 1992; Konzelmann and Braithwaite, 1995).

In studies of spatially distributed glacier mass balance (e.g. Klok and Oerlemans, 2004; Hock and Holmgren, 2005; Hock, 2005; Reijmer and Hock, 2007; Mölg et al., 2009; Rye et al., 2012; Gurgiser et al., 2013) optimization of free parameters to in situ measurements can be successful if the processes being parameterized are quasi-constant over the whole glacier surface, or if a dense measurement network is available for spatially distributed optimization. Brock et al. (2000) concludes that the accuracy of spatially distributed models is strongly dependent on the ability to apply multiple local optimizations, and on the importance of individual energy components. Nevertheless, most temperature index models (e.g. Hock, 2005; Pellicciotti et al., 2005; Carenzo et al., 2009; Robinson et al., 2010, 2011) and also a number of energy balance models (e.g. Mölg et al., 2009; Gurgiser et al., 2013) have been optimized towards a single best fit to the glacier-wide mass balance measurement, which requires a subjective choice of the single mass balance metric to be used. For example, optimizing for cumulative mass balance, mass balance gradient or stake measurements have been shown to be problematic as different optimal solutions are found depending on the mass balance metric chosen for optimization (Rye et al., 2012). The associated differences in the individual optimal parameter values and resultant values of the energy components have not been studied in detail, and furthermore, uncertainties of mass balance measurements (e.g. Zemp et al., 2013; Galos et al., 2017) imply that a single best-fit model simulation may not be found at all (Beven and Binley, 1992).

A way forward may be found in multi-objective optimization of glacier energy balance modelling, first applied in a glaciological context by Rye et al. (2012). They optimized a mass and energy balance model, on two Arctic glaciers in Svalbard over ∼40 years using three objectives for optimization: (i) the mass balance gradient, (ii) the mean absolute error (MAE) at the stake location, and (iii) the cumulative mass balance. This approach creates an ensemble of optimal solutions which all are equally “good” in respect to all three objectives. With this approach the authors could reconstruct the mass balance of the glaciers before direct measurements were available and also give an estimate of the model uncertainty from the parameter spread within the optimal solution set. This work demonstrated that it is likely that stated model performance based on single objective optimizations does not adequately represent model performance at a glacier scale or over longer time periods.

Figure 1 The sequential approach used in this study can be classified in three steps. First, data management and model setup in beige, the simulations (blue) first use a global sensitivity analysis to reduce the parameter space followed by a multi-objective optimization. All simulations are performed independently for three summers on two glaciers. The data analysis (green) is done independently for sensitivity, parameter and model uncertainty analyses.

Mass balance models are required to be transferable in space and time in order to estimate run-off on a larger scale or the impact of a changing climate (Oerlemans et al., 2005; De Woul and Hock, 2005; Raper and Braithwaite, 2006). To study the transferability of an enhanced temperature-index model Carenzo et al. (2009) used the optimized parameters from one particular year and glacier and compared it to the locally optimized run at different glaciers and over different time periods. They concluded that their model shows a rather good transferability in space, except during overcast conditions. Furthermore, they observed that the parameters vary depending on year and location and are correlated to each other. MacDougall and Flowers (2011) and Prinz et al. (2016) investigate transferability of full energy balance models; while MacDougall and Flowers (2011) find satisfactory temporal transferability in the Arctic over 2 years, albeit with some local parameter adjustment, Prinz et al. (2016) fail to do so in the tropics over an interval of a century. This is attributed to a substantially changed climate over the century and/or different micro-meteorological setting due to dramatic glacier shrinkage (Prinz et al., 2016), which implies the problem of transferring a calibrated model to rather different climatic settings and glaciers and raises the question about the general uncertainty and transferability of such models.

It can be expected that models with more parameters have greater variation in their solutions. Reduction of free parameters for optimization based on a sensitivity analysis is therefore a helpful tool to reduce both the effect of parameter correlation and computational expense (Spear and Hornberger, 1980; Saltelli et al., 2000; van Griensven et al., 2006). Gurgiser et al. (2013) applied such a parameter reduction procedure on a tropical glacier to reduce the free parameters prior to assessing model transferability.

Figure 2 Model simulations are performed at the stake locations shown as points; points marked in black are only used in the optimization, while green points indicate the seven stakes on each glacier that were also used in the sensitivity analysis. Detailed maps are available in the Supplement (Figs. S1–S2).

Model sensitivity and model uncertainty are often evaluated together, and assessments of varying robustness have been presented in the literature. For example, Mölg et al. (2012) used an arbitrarily chosen spread of the most positive and negative deviation simulations around their single best fit in respect to root mean square deviation (RMSD) of cumulative mass balance to estimate uncertainty. This gives only a rough estimate, as only two particular runs determine the uncertainty estimate. Anslow et al. (2008) first optimize their model and then vary the optimized parameters within certain bounds (5 %) and perturb the meteorological input to quantify the impact on the mass balance. This provides the sensitivity of the model output towards the parameter values and inputs, but the created range is also used as model error estimate. Machguth et al. (2008) perform a similar assessment but base their perturbation ranges on probability density functions, whereby model uncertainty is assessed by applying random and systematic errors and uncertainties to the meteorological input data as well as to the mean value of parameters. The reported uncertainty, of 700 kg m⁻² for a 400-day simulation at a single point (roughly 10 % of the total melt), is related to the standard deviation of the probability density function. Rye et al. (2012) used multi-objective optimization to better constrain their model parameters but do not evaluate their model on independent observations (i.e. observations not used for calibration).

In this study, we present a model calibration and uncertainty assessment workflow built upon a combination of these ideas. Our aim is to bring awareness that uncertainty estimates of physically-based models with many free parameters are likely to be under-estimated when applied in different settings (geographical and or temporal) than those for which the model was calibrated. Using an established distributed energy and mass balance model (Mölg and Hardy, 2004; Mölg et al., 2008, 2009), we simulate 3 years of summer mass balances on two mid-latitude glaciers (Fig. 1). We start by applying a global sensitivity analysis to reduce the parameter space extending the local sensitivity analysis used by Gurgiser et al. (2013) to a global variance-based method (Saltelli et al., 2006), a procedure which has recently been applied in snow pack modelling (Sauter and Obleitner, 2015). Subsequently, we use the multi-objective optimization applied by Rye et al. (2012) to calibrate our model based on a set of three quality measures. The parameter uncertainty and resulting uncertainty of the energy components are evaluated based on this calibration procedure. Finally, the temporal and spatial transfer of such a model ensemble is assessed with cross-validation.

In this paper, we will use the term “model uncertainty” to describe the difference between any modelled quantity and its counterpart in reality (“the truth”). An uncertainty value is a measure of how much trust can be given to a modelled quantity: in practice, model uncertainty can be estimated based on observations, and in any modelling activity which includes parameter calibration model uncertainty must be estimated separately from the calibration procedure (cross-validation). For quantities without equivalent in reality (e.g. model parameters), we use the term “uncertainty” to refer to the fact that their true value is really unknown, and that this uncertainty in the parameters is also conveyed in the model uncertainty. When we speak from “model sensitivity”, we mean the variance of the model output as function of the variance of an input quantity (e.g. forcing data, model parameters). A model sensitivity analysis does not require observations. In our paper, we restrict our sensitivity analysis to the internal model parameters, not to the input meteorological variables.

2 Study sites and model input data

Two glaciers in the eastern European Alps were selected as test sites in this study (Fig. 2). Hintereisferner (HEF; 46.80^∘ N, 10.75^∘ E) is a sizeable valley glacier in the Austrian Ötztal Alps that spanned 2454 to 3720 m a.s.l. in 2013, when the glacier area was ca. 6.7 km². Langenferner (Vedretta Lunga) (LGF; 46.46^∘ N, 10.61^∘ E) is a smaller valley glacier in the Italian Ortler Alps that spanned 3370 to 2711 m a.s.l. in 2013. These glaciers were chosen since the model used here requires topographic and meteorological input data, and measurements of surface mass balance for evaluation. For both these glaciers (i) topographic data are available in the form of high-resolution digital elevation models (DEMs) derived from airborne laser-scanning data acquired in autumn 2013 (Galos et al., 2015); (ii) meteorological data are available from automatic weather stations (AWSs) in the vicinity of the glaciers for the period 2012 to 2014 and (iii) intense glaciological observations, including measurements of seasonal mass balance (e.g. Klug et al., 2017; Galos et al., 2017), are available.

At HEF the AWS is located on a small plateau within a rock slope north of the upper tongue area of the glacier at an altitude of 3025 m a.s.l. The horizontal distance of this AWS from the glacier is about 300 m and it provides all meteorological data required for the model except for precipitation. Precipitation data were taken from the gauge operated by the Bavarian Academy of Sciences at Vernagtbrücke, 3.5 km east of HEF at an elevation of 2600 m a.s.l., and were scaled to the elevation of the AWS on the basis of precipitation gradients derived from 11 totalizing rain gauges in the vicinity of the glacier (Strasser et al., 2017). At LGF the AWS data came from the station of the Hydrological Service of the province of Bozen operated at Sulden Madritsch, 2.5 km north of the glacier at an altitude of 2825 m a.s.l. (Galos et al., 2017).

3 Model and methods

3.1 Energy balance model

The energy and mass balance model used in this study is a process-based model that has been applied in a range of glacier environments (Mölg and Hardy, 2004; Mölg et al., 2008, 2009, 2012; Gurgiser et al., 2013; Prinz et al., 2016; Galos et al., 2017). The model was run in hourly time-steps for three summer periods over each glacier. The model is a distributed mass and energy balance model, but in this study simulations were limited to 18 stake locations on each glacier to reduce computational expense. The model tracks the accumulation of solid precipitation and uses the surface energy balance to calculate the ablation at the glacier surface:

$$\begin{array}{}\text{(1)}& Q_{\text{M}}+Q_{\text{ice}}={\text{SW}}_{\text{net}}+{\text{LW}}_{\text{net}}+Q_{\text{S}}+Q_{\text{L}}+Q_{\text{G}}+Q_{\text{P}},\end{array}$$

where LW_net and SW_net are the net radiation balances for long-wave (thermal) and short-wave (solar) radiation and the other energy fluxes are the sensible (Q_S), latent (Q_L), ground (Q_G) and precipitation (Q_P) heat flux. The available energy is used to raise the glacier surface temperature (Q_ice) if below freezing point or for melting (Q_M) if the glacier surface is at the melting point. Mass losses of the glacier are represented via melt (Q_M) and sublimation (Q_L). Refreezing of liquid precipitation and resublimation lead to additional mass accumulation at the surface. We use the model in a similar configuration to Prinz et al. (2016). The only difference is given by a change in the short-wave radiation scheme which is explained in the detailed model description in the Appendix (Sects. A1–A6).

3.2 Methods

3.2.1 Global sensitivity analysis (GSA)

Variance-based sensitivity testing methods work in a probabilistic framework judging sensitivity by relative variances of model input and output (van Griensven et al., 2006; Saltelli et al., 2000, 2006, 2010). This is a global method that is independent of model calibration, i.e. independent of a local optimal run, and is hereafter referred to as global sensitivity analysis (GSA). The method treats the model as a simple function f with

$$\begin{array}{}\text{(2)}& y=f\left(X\right)X=X_{\mathrm{1}},X_{\mathrm{2}},\mathrm{\dots },X_n,\end{array}$$

where y is the single model result (in this case mass balance) and $X_{\mathrm{1},\mathrm{\dots },n}$ are the individual input parameters.

The influence of an individual parameter can be examined by the main effect (V_i) of X_i on Y.

$$\begin{array}{}\text{(3)}& V_i=V_{X_i}\left(E_{X-i}\right(Y\left|X_i\right))\end{array}$$

X_−i is the whole parameter space except any variation in X_i (a fixed X_i), E is the expectation value and V the variance. $E_{X-i}\left(Y\right|X_i$) is the mean model output with whole parameter variation except in X_i. The variance over all values for X_i yield the variance attributed to parameter X_i. The sensitivity of the model towards single parameters is evaluated by normalizing by the total variance of the output.

$$\begin{array}{}\text{(4)}& S_{X_i}={\displaystyle \frac{V_{X_i}\left(E_{X-i}\right(Y\left|X_i\right))}{V_y}}\end{array}$$

$S_{X_i}$ is the first-order sensitivity index. The total sensitivity index ($S_{{\text{T}}_i}$) is the effect of X_i with all its interactions on the model variance:

$$\begin{array}{}\text{(5)}& S_{{\text{T}}_i}={\displaystyle \frac{E_{X-i}\left(V_{X_i}\right(Y\left|X_{-i}\right))}{V_y}}.\end{array}$$

This can be related to the sensitivity obtained from local sensitivity analysis. The model sensitivity (variance) to X_i is tested ($V_{X_i}\left(Y\right|X_{-i}$)) at every point of the parameter space (X−i fixed). To clarify, consider the example of a simple non-additive model $Y=X_{\mathrm{1}}\cdot X_{\mathrm{2}}+X_{\mathrm{3}}$ with the variables X_i as input parameters with a given variance/uncertainty. Assuming unified distribution within the intervals

$$X_{\mathrm{1}}\in [\mathrm{1},\mathrm{3}],X_{\mathrm{2}}\in [\mathrm{0.1},\mathrm{0.3}],X_{\mathrm{3}}\in [\mathrm{0.5},\mathrm{1}],$$

leads to a model output range of $Y\in [\mathrm{0.6},\mathrm{1.9}]$. The variance-based method yields the results for $S_{X_i}$, the first-order sensitivity index and $S_{{\text{T}}_i}$, the total sensitivity index for an ensemble of 10 000 runs as shown in Table 1. The first-order effect of X₃ is the largest, while the other two are similar if computational uncertainty is neglected. The most variance is caused by the last parameter. X₃ has no interactions, so its total index is the same as the first-order one, while interaction between X₁ and X₂ creates additional variance, so their total index is higher. In the example X₁ and X₂ contribute to ≈60 % of the total variance and X₃≈40 %, as $X_{\mathrm{1}}\cdot X_{\mathrm{2}}\in [\mathrm{0.1},\mathrm{0.9}]$ and $X_{\mathrm{3}}\in [\mathrm{0.5},\mathrm{1}]$.

Table 1 The sensitivity indices for the simple model $Y=X_{\mathrm{1}}\cdot X_{\mathrm{2}}+X_{\mathrm{3}}$. The indices for X₁ and X₂ are similar as they both have the same normalized variance. X₁⋅X₂ creates additional variance by the interaction of the two parameters yielding higher total indices.

Table 2 The ranges for the 22 different parameters used in the sensitivity study. Most parameterizations are explained in the Appendix A.

The estimation of the sensitivity indices follows the algorithm from Saltelli et al. (2010). The model used here has 22 free parameters. A base sample of 12 000 parameter settings was created with a quasi-random Sobol sequence. The random numbers are linearly transformed onto the parameter intervals. The distribution is always treated as uniform and the limits for every parameter are given in Table 2. The indices are estimated with $N\cdot (k+\mathrm{2})$ runs, where k is the number of parameters and N the base sample size. The GSA consists of a total ensemble size of 300 000 simulations per year and glacier, fulfilling the convergence criteria for the algorithm ($S_{{\text{T}}_i}\ge S_{X_i}$, $\sum S_{X_i}\le \mathrm{1}$, $S_{X_i}\ge \mathrm{0}$). Note that we did not investigate if fewer solutions could already fulfill the convergence criteria. To reduce computational expenses the GSA model was limited to seven stake locations on each glacier (Fig. 2).

The parameter sensitivity results from the GSA are also used as a tool to reduce the number of free parameters in the model by identifying those parameters which have only a marginal influence on the model output (Spear and Hornberger, 1980; Saltelli et al., 2000; van Griensven et al., 2006). The model is considered insensitive to parameters with a total sensitivity index ($S_{{\text{T}}_i}$) of <0.05, and these parameters were fixed at the median value of the range shown in Table 2 in subsequent model simulations.

3.2.2 Multi-objective optimization and uncertainty quantification

A multi-objective optimization allows for more than one optimal solution in the calibration procedure, and offers a way to assess a range of plausible parameter sets that we will use later on for model predictions. The multi-objective optimization used here follows previous approaches in hydrology and glaciology (Yapo et al., 1998; Rye et al., 2012). Where the model is given n objectives, with f_n to be minimized in respect to the model parameter input X, the optimization approach can be written as follows:

$$\begin{array}{}\text{(6)}& \text{minimize}\left(f_{\mathrm{1}}\right(X),f_{\mathrm{1}}(X),\mathrm{\dots },f_n(X\left)\right).\end{array}$$

The result of Eq. (6) is an ensemble of optimal solutions that represent trade-offs between the objectives and no single one can be deemed superior to the other optimal solutions. Therefore, they are called the non-dominated set of optimal solutions, or Pareto set (Pareto, 1971). As an illustration, consider an optimization with two objectives (f₁,f₂): The concept of a Pareto optimal set is shown in Fig. 3 in which the (classic) single objective solutions are the points $f_{\mathrm{1}}^{min}$ and $f_{\mathrm{2}}^{min}$ for the two objectives. A solution at the utopian point is desirable as all functions would be at their minimum, but the models generally cannot optimize the different objectives simultaneously. There are only compromise solutions between the objectives. The members of the set of optimal solutions defining the Pareto front are superior to the other solutions, but are all equal to each other without subjective ranking by the modeler. The variation of the parameters of the optimal solution set defines the minimum parameter uncertainty (Vrugt et al., 2007). This uncertainty is a result of shortcomings in the model and/or the variations of parameters, such as spatial or temporal change in the true parameter value over the simulation period (Oerlemans and Greuell, 1986; Marshall and Warren, 1987). If a single simulation must be chosen to be the optimal model set up, the compromise solution, defined as the point with the lowest Euclidean distance to the utopian point, is a common choice.

Figure 3 The figure displays a 2-dimensional Pareto space which comprises a 2-dimensional Pareto front. The solutions on this front (black solid line) are referred to as the non-dominated set of solutions. In comparison, all other solutions within the solution space are inferior in at least one objective relative to the Pareto front. Classic single objective optimization yields the points $f_{\mathrm{1}}^{min}$ and $f_{\mathrm{2}}^{min}$, which represent the minimum of those objectives that the model can achieve. The utopian point (black) is the point ($f_{\mathrm{1}}^{min},f_{\mathrm{2}}^{min}$) where both objectives are at their minimal value. Commonly, the compromise solution (red) of the Pareto set is considered an objective choice for a single solution as it has the minimum Euclidean distance of the optimal solution towards the utopian point.

In this study the multi-objective optimization is based on a Monte Carlo simulation. The non-sensitive parameters from the GSA were fixed to their median value from the range used in the GSA (Table 2). Then 20 000 model simulations with random value combinations of the remaining parameters were created and the mass and energy balance were simulated for 18 stake locations. This approach was chosen in favour of an evolutionary algorithm so that different objective function spaces and all single objectives could be investigated with the same set of simulations. Various objective functions were initially explored, including root mean square deviation (RMSD) and mean absolute deviation (MAD) over all simulation points, but finally three objective functions that captured the main patterns of behaviour were applied: (i) the BIAS over all simulated stakes, (ii) the mean absolute deviation (MAD) of the lower nine stakes (MAD_low9) and (iii) the MAD over the upper nine stakes (MAD_top9). The BIAS is used as a proxy for the cumulative mass balance with the avoidance of interpolation errors. The RMSD is a commonly used measure for optimization in glaciological modelling (e.g. Gurgiser et al., 2013; MacDougall and Flowers, 2011). By using the MAD here we want to reduce the effect of individual stakes which could be influenced by processes which are not captured by the model (snow redistribution through wind or avalanches, dust and debris cover and related changes in radiation, etc.), but the general features of those two statistical functions are similar. Previous studies (e.g. Klok and Oerlemans, 2004; Hock, 2005; Sauter and Obleitner, 2015) have focused on the accumulation and ablation area separately or exclusively, but without a distinct mathematical comparison. Therefore, the approach of the split MAD was chosen. The Pareto front was identified, and a second ensemble including solutions within a certain range (100 kg m⁻²) from the Pareto front was additionally identified to account for errors in the field measurements of mass balance at each simulation point. However, results of this second ensemble will only be mentioned briefly throughout the discussion. The spread of the parameter settings of all optimal solutions of the Pareto and near-Pareto Sets are used to indicate the parameter uncertainty for each case, and the calculated surface energy balance components of these optimal sets are also used to estimate the uncertainty of the energy components on the point scale, as well as on the glacier scale.

4 Results and discussion

4.1 Global sensitivity analysis

Figure 4 The amount of sensitive stakes per year for (a) HEF and (b) LGF. The sensitivity analysis was performed at seven stakes on each glacier, though the vertical gradients can only be tested at six stakes as one is located at the same altitude as the reference weather station. Every parameter with a sensitivity index higher than 0.05 got a score of 1, giving a maximum count of seven per year (meaning the model is sensitive to this parameter at all stakes). Parameters involved in the parameterization of surface albedo are dominating, with snow related parameters in the upper section of the glacier and ice related ones at the lower stakes. Hintereisferner shows a total of 11 sensitive parameters and Langenferner 6.

The focus of this GSA is not on the absolute sensitivity towards single parameters, but rather to reduce the dimension of the parameter space. Therefore, the following discussion is limited to two classes: parameters to which the model is sensitive ($S_{{\text{T}}_i}>\mathrm{0.05}$) and non-sensitive ($S_{{\text{T}}_i}<\mathrm{0.05}$). On each glacier the mass and energy balance at seven stake locations over 3 years was simulated for the GSA, so the maximum count of sensitivity for a parameter would be 21, meaning that the model is always sensitive to that parameter at every point of the glacier.

At Hintereisferner, 11 out of 22 parameters are identified as sensitive (Fig. 4a), and these sensitive parameters are classified in two general categories. Firstly, all but the lowest stake location are sensitive to parameters related to surface albedo, particularly of snow and firn, and secondly, for stakes with high elevation differences compared to the AWS, the model is also sensitive to the vertical temperature gradient.

The sensitivities show spatial and temporal variability, which can be explained by the varying mass balance conditions of the respective year (mean specific summer/annual mass balance with 2012: −2643/−1560, 2013: −1841/−510 and 2014: −1494/−122 kg m⁻²). For example, sensitivity to the ice-related parameters is most evident in 2012, which was the driest (in terms of precipitation, not air humidity) and most negative mass balance year, with large parts of the glacier surface free of snow and firn for most of the ablation season. By contrast, the roughness length of fresh snow is only influential at the upper stakes in 2014, where snowfall was frequent during the ablation season resulting in the least negative mass balance of the three study years. Sensitivity towards the elevational precipitation gradient is only relevant at the lowest stakes (500 m below the weather station) in the wet years.

On the smaller Langenferner, 6 of the 22 parameters were identified as sensitive (Fig. 4b). Similar to HEF, the model shows consistent sensitivity to surface albedo and the vertical temperature gradient. As LGF is smaller than HEF, the sensitivity shows less variability in space and time, though the annual mass balances during the three study years range from about −1500 to +400 kg m⁻², and, as the tongue of LGF does not extend to such a low elevation as the one of HEF, it is less sensitive to ice-related parameters. Variations in the ice albedo within the bounds of 0.15 and 0.25 hardly influence the mass balance model results on the smaller glacier, even though ice is exposed for the majority of the summer at the lowest stake. This low sensitivity to the ice albedo compared to the snow albedo parameters is explained by the fact that, as the removal of snow cover is accompanied by a large drop in albedo (0.4–0.65 to 0.15–0.25), the time of exposure is more crucial than the final ice albedo, and this time of ice exposure is itself influenced by the snow albedo via its dominant control on the short-wave radiation budget. Within the chosen parameter ranges, the net short-wave radiation varies by 50 % in the case of fresh snow (10 %–20 % absorbed) and only by 12 % over ice.

4.2 Calibration

Table 3 Five objective functions are used to analyse the model performance. The minimum value for every function and each year are given in kg m⁻². While the BIAS is low in all cases, absolute errors and RMSD are much higher, and highest in 2012. Note that the minimum MAD does not refer to the same run over the whole glacier and its upper and lower parts.

First we consider the best model performance with respect to each individual objective function tested (Table 3), before presenting the multi-objective optimization based on the first, fourth and fifth objective in Fig. 5.

In all cases a model simulation with very low bias (<1 kg m⁻²) with respect to the stake mass balance can be found. This illustrates that apparently a good optimization on the single value of cumulative mass balance over the stakes is relatively easy to achieve (Table 3). In comparison, the deviation of all other objective functions is much higher, ranging from 81 to 537 kg m⁻². The deviations in these objectives are all largest in 2012 on both glaciers. RMSD and MAD vary in a similar way between the years at each glacier, with the higher RMSD values indicating a non-uniform deviation from the measurements over the stakes. With the exception of 2014, the glacier-averaged MAD is larger than the MAD calculated for either the upper or lower section of the glaciers. This is to be expected as the stakes within each section of the glacier experience more similar climate conditions, resulting in a lower MAD. The fact that MAD in the lower glacier section is larger than in the upper section in 2012 and 2013 is probably related to the incapability of the model in its current configuration to correctly reproduce the date of ice exposure. In 2014 the upper glacier sections show a slightly higher MAD, associated with above average accumulation in the previous winter and the frequent summer snowfall in this season.

The multi-objective optimization, using BIAS, MAD_top9 and MAD_low9, yields an ensemble of solutions. The non-dominated set for each of the 3 years has 27, 17, 69 members for HEF and 58, 61, 14 members for LGF, respectively (Fig. S4). The fewest solutions are found in years with the lowest total MAD (HEF, 2013; LGF, 2014). Figure 5 shows the Pareto front of optimal solutions for HEF 2012 and the corresponding parameter settings. A low bias is easily achieved by the model if no other objectives are considered because it is a single value (the sum of the mass balance at all stakes) and, for example, deviations in the ablation and accumulation area may compensate each other. The projections onto the BIAS planes are less curved (the distance between the utopian and compromise point is smaller) and the performance in respect to the MADs can be drastically increased with only a small cost in the BIAS. The 2-dimensional projections of the Pareto space (Fig. 5a and b) illustrates, for example, that allowing for a model bias of 25 kg m⁻² can improve the MAD by 200 and 300 kg m⁻² in the lower and upper glacier sections, respectively. The MADs plane (Fig. 5c) is more curved (larger distance between the utopian and compromise point), indicating that the two objectives cannot be optimized by the model at the same time, such that some parameter sets leading to good results for the ablation zone of the glacier may not sufficiently reproduce the relevant processes in the accumulation zone.

Figure 5 Each individual member of the Pareto set for HEF in 2012 is displayed with a different colour and the compromise solution highlighted (red triangle or red line). The different panels are the 2-dimensional projections of the Pareto space onto the (a) BIAS and MAD_top9; (b) BIAS and MAD_low9; and (c) MAD_low9 and MAD_top9 planes. Panel (d) shows the normalized parameter values for each case in the same colours as in the Pareto space plots. The parameter settings of the optimal solutions are quite diverse and span over most of the parameter space. The firn albedo and albedo timescale are the only parameters showing some confinement to a narrower range. The chaotic nature of the parameter settings furthermore shows that a single solution is not representative in its parameter settings for the ensemble of optimal solutions.

The parameter values of those optimal solutions span the entire allowed space apart for some of those relating to snow albedo which span (almost) the whole parameter space in all years for both glaciers, and show no obvious tendencies towards a certain albedo range (Sect. 3). For HEF in 2012, snow albedo values cluster in the higher range (0.52–0.6) for firn and (0.86–0.9) for fresh snow (Fig. 5d), while on LGF lower firn and fresh snow albedo values (<0.5/<0.84) are optimal. Similar behaviour is observed for the albedo time scale (see Appendix A3), which tends towards higher values for HEF in 2012/13 and towards lower ones for LGF in 2013/14. The confinement of snow albedo is mainly a result of the highest model sensitivity towards this parameter, nevertheless it still varies and the converse argument, of less sensitive parameters showing greater span is not valid: for example, the roughness length over fresh snow is generally at the lower margin of the allowed parameter range (0.1–$\mathrm{0.14}\times {\mathrm{10}}^{-\mathrm{3}}$ m) in 2012, even though the model is considered insensitive ($S_{\text{T}{z}_{{\text{h}}_{\text{fs}}}}<\mathrm{0.05}$) to this parameter in that particular year. These results highlight that the parameter settings of multiple optimal solutions for this type of mass and energy balance models can vary drastically. There are no clear correlations between two individual parameters, instead all parameters interact simultaneously to some degree. Without the a priori reduction of model parameters by GSA, even less information could be extracted from the optimization. Compared to Rye et al. (2012) our parameters span a wider range of the normalized parameter space, which is due to a wider initial parameter range in our study. Despite the relatively narrow ranges of values reported in the literature, our study clearly reveals that many of the parameters could take almost any value in the optimization process. Changes to the parameter ranges accounting for potentially unrealistic values may quantitatively change the results, but within the range no change in the sensitive parameters is expected. For example, Rye et al. (2012) applied values for fresh snow albedo in the range of 0.65 to 0.95, while we restricted the initial range to values between 0.8 and 0.9, as reported in the literature (e.g. Cuffey and Paterson, 2010). We also used fresh snow densities which are relatively low compared to those reported in recent studies Helfricht et al. (e.g. 2018). However, the used values are based on previous studies (e.g. Mölg et al., 2008; Gurgiser et al., 2013; Prinz et al., 2016) and the choice of those does not significantly influence the results.

4.3 Transferability studies

To investigate the transferability of the optimized mass balance model settings, all the optimal solutions of the Pareto set of one glacier summer mass balance case were applied to the five other summer and glacier cases. While each Pareto set was identified based on the multi-objective optimization, the transferability study uses only the Euclidean distance towards the utopian point as a quantification tool. The individual optimal parameter settings for HEF 2012 for example yield quite varying performances for the other summers (Fig. S4a). While the performance on the same glacier (HEF) is reasonably good for 2012 (200–800 kg m⁻² compared to 152–600 kg  m⁻² in the optimization period of summer 2013) and slightly worse for 2014, the optimal solutions do not perform so well for LGF, resulting in Euclidean distances of up to 3500 kg m⁻². Analogous analysis of the ensemble behaviour of other summers shows that the optimal solution for 2012 also performs well in 2013 and vice-versa, and shows acceptable performance in 2014. The deviation of the 2012 and 2013 optimal values of HEF yields errors greater than 2000 kg m⁻² on LGF. The 2014 HEF ensemble performs on average better on HEF, but two simulations perform better on LGF in 2012/13 and around 20 are within the same error as for HEF. On LGF also 2012 and 2013 agree better, and the ensembles produce reasonable results for both glaciers in 2014. The ensemble of 2014 on LGF yields similar errors (250–800 kg m⁻²) for LGF 12/13 and HEF 14. All ensembles of LGF produce larger errors on HEF in 2012 and 2013.

Figure 6 (a) Performance of the single compromise solution for each season and glacier (GGG_yy), with HEF in solid blue colours and LGF in red. The simulations which perform best over a 3-year (RMSD_HEF and RMSD_LGF) and 6-year period (RMSD_all) are given with dashed lines, following the same colour scheme. (b) The corresponding parameter setting of the optimal solutions to the left. The colour scheme is equivalent. The compromise solutions for the individual years show different parameter settings and also varying performance out of the calibration period. Only the snow-albedo related parameters show a trend as they take rather large values on HEF and small on LGF. No clear trend is visible for the other parameters.

The cross-validation (Fig. 6) focuses on the transferability of the single compromise solution to other season and glacier cases. This can be considered as a classical best guess solution. The features follow the structure of the ensemble behaviour discussed above with HEF 2012 and 2013 seeming to be distinct from the other four cases. The compromise solutions for HEF 2012 and 2013 are similar in performance and parameter value and, while they perform adequately for HEF in 2014, within the estimated model uncertainty of 1300 kg m⁻², the error is greater than 1500 kg m⁻² when either of these compromise solutions is applied on LGF, no matter for which year. Similarly, the compromise solution for the 3-year period for HEF (RMSD_HEF in Fig. 6), which is dominated by the characteristics of 2012 and 2013, also performs poorly when applied to LGF (errors of up to 3500 kg m⁻²). The compromise solution of HEF 2014, however, generally performs better on LGF than for other years at HEF, and reciprocally, the compromise solution over the whole period at LGF performs best at HEF in 2014, and the maximum error (up to 2500 kg m⁻²) is lower than for cases of HEF compromise solutions being applied to LGF. This is probably due to the domination of more negative mass balances in 2012 and 2013 at HEF, where good model performance is linked to capturing the large extent of the ablation area, whereas the shorter glacier tongue at LGF has smaller impact on the mass balance of this glacier. The compromise solution (RMSD_all) for all six cases also highlights that within this set of six the cases HEF 2012 and HEF 2013 are more distinct from the other cases as the overall compromise solution performs worst in these two cases. For most parameters no clear separation between the two glaciers is evident, except for fresh snow albedo and the albedo timescale which are both larger at HEF and smaller at LGF. Inspection of the optimal parameter values reveals that runs with a longer calibration period (RMSD_xxx) do not necessarily take trade-off values between the individual years. For example, in this case the solution that performs best over both glaciers and the whole time period (RMSD_all) takes larger values of fresh snow density and ice-albedo than any other compromise solution (Fig. 6b). This further highlights the model complexity and is suggestive of the effects of physical shortcomings (such as parameter values that are constant in space and time) compensating each other.

Figure 7 The energy balance components for 8 out of 18 selected stake locations close to the central flow line, displayed in different colours for HEF 2012. Solid bars represent the fluxes of short-wave radiation (SW), long-wave radiation (LW), turbulent heat fluxes (Q_L and Q_S), penetrating short-wave (Q_PS), precipitation heat flux (Q_P),conductive heat flux (Q_C) and the resultant available heat for melting (Q_M, here plotted as a positive flux).

4.4 Energy balance components

Analysis of the energy balance components associated with Pareto set solutions offers a qualitative means of verifying that the identified optimal parameter settings are in line with expected physical processes at the glacier surface. The energy balance components calculated by the model are expected to vary depending on the parameter settings of an optimal ensemble, which have been demonstrated to span almost the whole parameter space. This variation in energy balance components is indicative of the uncertainty in the modelled energy fluxes (we say “indicative”, because the true uncertainty can only be assessed using observations, which are not available here). Figure 7 illustrates such variations in the energy balance components for the case of HEF 2012 based on the our model, not accounting for uncertainties in the meteorological input itself. In this case, the most uncertain energy balance component is the short-wave radiation, which at the same time is the largest energy source for the surface. Total energy flux from short-wave radiation decreases with altitude, while the associated uncertainty increases. The sensible and latent heat flux provide a net energy source to the surface and their value and uncertainty also decrease with altitude. The long-wave radiation budget is a net energy loss from the surface in summer, and its value increases and uncertainty decreases with elevation. As a result of these elevational patterns in uncertainty, the uncertainty in melt energy is also largest at low elevations.

The variation of the averaged energy components over the stakes for HEF 2012 are given in Fig. 8. The uncertainties are generally lower than on a stake basis. The short-wave, conductive ground heat flux and sensible heat flux supply a net heating to the surface on both glaciers. The precipitation heat flux is also a minor energy source. The penetration of short-wave radiation and the long-wave budget remove energy from the glacier surface. Latent heat is the only energy flux that has either a positive or negative effect on the surface energy balance depending on stake location, glacier and year. On both glaciers, lower elevation locations tend to show more positive energy fluxes from latent heat. At HEF this flux is mostly an energy addition to the glacier surface while on LGF it mostly serves to remove energy from the surface. In the beginning of the summer, sublimation during the day and condensation or re-sublimation during the night is dominant on HEF, and the general trend over the summer is to progressively more condensation. LGF shows less condensation during mid-summer, which is mainly attributed to less windy conditions than at HEF.

Figure 8 The energy balance components averaged over all stakes has less uncertainty than on the point scale for HEF 2012. The objective functions are all integrated over the whole glacier and thus the uncertainty is lower. Glacier wide the short-wave radiation is the largest component with the largest absolute uncertainty as well, followed by the turbulent fluxes. The long-wave balance and the penetrating short-wave radiation provide a net cooling effect for the surface.

Table 4 The energy balance components are averaged over all stake locations. The uncertainty is given in respect to the minimum and maximum of the ensemble. The short-wave radiation (SW_net) has the largest impact, decreasing in importance from 2012 to 2014, with a less negative mass balance (Q_M). The penetrating short-wave radiation (Q_PS) follows the same pattern with the opposite effect. The long-wave budget (LW_net) is lower for LGF. The turbulent fluxes are greatest in 2012 and larger on HEF. The precipitation (Q_P) and convective (Q_C) heat flux are of minor importance.

The total contribution of the energy balance components averaged over the glacier are listed in Table 4. The relative uncertainties of the energy balance components are up to 50 % of their contribution on single stake basis and 30 % averaged over HEF; slightly lower (30 % and 25 %, respectively) for LGF. This leads to a variation in the available heat for melting and the mass balance of about 30 % on a point scale. The absolute uncertainty of the seasonally averaged available energy for melting can reach up to 35 W m⁻² at the tongue area of HEF. This corresponds to a daily melt uncertainty of 9 kg m⁻² and seasonal uncertainty of up to 1.3 kg m⁻². The glacier averaged available heat for melting is much less uncertain over all stakes. This is a result of the calibration process. The sum of total available melt energy is directly linked to the bias as objective function, which shows the largest value among the optimal solutions on HEF 2012 with 600 kg m⁻². In comparison the MADs which are more influenced by the mass balance at the individual stake reach values up to 1000 kg m⁻².

The largest uncertainties in our study are associated with the short-wave radiation as a result of the albedo parameterization, which relies on five model parameters. Alternative albedo parameterizations are also known to be a source of substantial uncertainty (e.g. Klok and Oerlemans, 2004; Willeit and Ganopolski, 2017). The greatest uncertainty is commonly found in the accumulation area and around the equilibrium line altitude. This is because (i) the parameterization for snow albedo has more variation and free model parameters than albedo over ice, and (ii) around the ELA the variation of the ice exposure date increases the uncertainty of short-wave radiation flux. Point scale albedo measurements combined with localized optimization schemes may solve this issue, but for distributed models a more detailed model may be necessary to better capture the full complexity of the processes governing initial snow albedo and its change through time (e.g. Flanner and Zender, 2006).

The long-wave radiation shows a lower uncertainty in this study than in Sauter and Obleitner (2015) and its uncertainty is mainly due to the air temperature, the related temperature gradient parameter, and the surface temperature. It is important to note that we cannot state that the general uncertainty of energy balance models associated with incoming long-wave radiation is low, since in this study the parameterization was optimized prior to the sensitivity analysis as direct measurements are available at the weather station. Consequently, long-wave radiation is considered a meteorological forcing here and therefore it was decided to do this prior optimization. The parameterization gives no bias for the station but the hourly RMSD was up to 30 W m⁻², which is in the range of the net long-wave budget. This therefore also mainly influences short-term differences in the long-wave budget rather than the seasonal energy flux. Nevertheless, as with albedo, it remains unclear whether long-wave radiation modules based on air-temperature, cloudiness and sky-view factor are sufficient to model spatio-temporal variation over a glacier.

The turbulent fluxes are associated with the second largest uncertainty in this study, which is in agreement with other studies that found larger uncertainties in the radiative forcing (Willis et al., 2002). Turbulent fluxes are important for determining short-term variations of melt rates due to, for example, changes in the stability regimes (Lang, 1981). However, the uncertainties in our model are due to differences in roughness lengths and the temperature gradient. Roughness lengths over ice and snow vary substantially (e.g. Braithwaite, 1995) in space and time (Greuell and Konzelmann, 1994; Calanca, 2001). The appropriateness of using constants for these values in glacier modelling is also questionable, and stability corrections may differ from the glacier margins to the interior, for example. It is therefore also questionable how appropriate constant roughness lengths and stability corrections for ice and snow in space and time are. Furthermore, recent studies (Sauter and Galos, 2016) showed that the application of the bulk approach in complex mountain terrain can generally be problematic.

Heat supply from rain is negligible in our study, which is in agreement with other studies on alpine glaciers (e.g. Hock, 2005).

4.5 Implications of this study

The larger glacier, Hintereisferner, has more sensitive parameters and the variation over the stakes is larger than at Langenferner, as a result of more distinct climate regions on the longer tongue of the larger glacier. This is also true for the uncertainty of energy balance components, with the exception of the net solar radiation, which is comparable on both glaciers. Short-wave radiation is the most uncertain of the energy balance components, due to the albedo parameterization, which accounts for the change in albedo over time, but does not account for any possible spatial variation in temperature or grainsize-dependent albedo decay rates. We have shown that the model has difficulties optimizing the upper and lower part of the glacier simultaneously, as a result of the variable parameter values of physical quantities like albedo. The large spread of our ensemble is a result of trade-off solutions between the real albedo at any time and any location and the temporally and spatially averaged parameterization applied. Other parameterizations that are assumed as constant in space and/or time, or only indirectly affected by temperature and altitude dependencies, are also subject to similar trade-off effects. Although the physical relations may not be the same at all times and the lower tongue area may be quite different from the upper glacier, this does not mean that the model performance is worse on the larger glacier (HEF) with more variation in a quantitative matter (Table 3), but rather that the solutions of the Pareto front show more variation in the parameter settings. This analysis clearly identifies the issue of governing parameters and parameterizations not being constant in space and time as the main problem of distributed energy balance modelling. The most readily appreciable example in this regard is ice albedo, which is often lower near the terminus due to debris and dust accumulation and water saturation of the glacier surface.

To improve this we suggest two potential approaches: (1) although optimizing all key parameters serves a purpose for a broad range of applications, fixing low sensitivity parameters to common values, which are not optimized, results in a type of a simplification of the model that reduces over-fitting and potentially increases the stability and comparability of the energy balance model over short-timescales. The overall performance of such a model will be lower because the tuning possibilities have been restricted, but better estimates of the model uncertainties for out-of-sample periods can be generated. (2) Parameters or parameterizations could be allowed to vary in space and/or time. This could be achieved either by increasing the measurements and data availability or increasing the model complexity. More complex albedo schemes are available, for example, for snowpack models like Crocus (Vionnet et al., 2012) or SNOWPACK (Lehning et al., 2002). However, if new parameterizations are introduced they require sufficient field data to constrain the physical process and should not be just added as additional model free parameters to optimize.

The approaches in this study are helpful tools to combine these suggestions. A clear understanding of the model sensitivity, independent of the optimization of the model is necessary to decide on the importance of certain parameters. It gives the option of fixing parameters and focusing on the key processes. We have shown that multi-objective optimization is a valid tool to assess uncertainties in the model. The objectives used are all based on the same data (i.e. stake data). This allowed us to show the uncertainty that is just associated with treating the available data in a different way without requiring additional measurements. The model can readily be optimized to minimise bias or meet any single value objective; therefore, model performance based on single best-fit approaches should be treated with caution. Furthermore, a single solution may suffer significantly from parameter over-fitting and is not representative in its parameter settings for other plausible solutions. The chosen objectives show that there is inter-annual variation in the performances of the upper and lower section of the glacier in our cases. The curved nature of the Pareto front highlights that simultaneous optimization of both areas is difficult for the model. Parameters are just not constant in either space or time, so the model uncertainty increases when the model is applied to other time periods or on another glacier. The model uncertainty is in the range of 1000 kg m⁻² per summer season for each glacier. It is larger when transferring the calibrated model to another alpine glacier, but still of the same order of magnitude. Our results reveal larger model uncertainties related to spatial transfer than found in previous studies (MacDougall and Flowers, 2011). This can be explained by the relatively large inter-annual variability of mass balance, as well as the comparably large distance between the glaciers in our study. Together with an uncertainty estimation of the energy balance components the key parameterizations, which need further improvement, can be identified. Furthermore, within the multi-objective framework it is possible to focus on processes individually: for example if the albedo is measured on the point scale, the difference to its model value could be used as an objective, instead of an a priori calibration of the albedo parameterization itself.

Neither meteorological forcing on the point scale nor mass balance measurements are free of errors, and the related model uncertainties were not formally disentangled from other uncertainties in this study. Zemp et al. (2013) have estimated an annual measurement uncertainty of 140 kg m⁻² on point scale glaciological mass balance measurements, while Galos et al. (2017) report somewhat lower values for Langenferner. More information about the propagation of those errors is needed to quantitatively include them in the optimization. However, if an uncertainty of 50 kg m⁻² in the MAD and BIAS is included, the Pareto sets increase by 1 order of magnitude. This complicates further interpretations and increases the total model uncertainty.

The analysis presented here indicates that while mass and energy balance models help us to understand the physical processes on the glacier, the necessity for parameterizations within these models introduces considerable, variable uncertainty to the model output. Calibration of surface mass balance models is complex and uncertainty studies are helpful to understand those models, and it is not advisable to draw general conclusions from such modelling efforts without first fully understanding the inherent model sensitivity and the properties of the uncertainty of the calculated mass balance and associated energy fluxes in detail.

5 Conclusions

Based on a well developed mass and energy balance model, applied to two well-studied glaciers in the European Alps, this study gives a robust estimate of the model uncertainty and discusses the advantages of parameter space reduction and multi-objective optimization in glaciological modelling.

Using a variance-based global sensitivity method, model sensitivity to the free model parameters was identified, independent of the calibration data. Model sensitivity to specific parameters is both site- and time- specific, and this should be acknowledged in wider applications of such models. By separating the parameters into two sensitivity categories, the model parameters to be optimized can be reduced. Those that the model output is sensitive to were subject to a multi-objective optimization, while non-sensitive parameters were fixed to literature values.

The multi-objective optimization was based on three objectives related to stake mass balance data measured using the glaciological method. We used the model bias over all stakes and the mean absolute deviation over the upper and lower part of the glaciers. It proved difficult to optimize model performance in the upper and lower section of the glacier simultaneously. The bias over all stakes, which was used as a proxy for the cumulative mass balance, can be minimized easily, and this should be considered when optimizing for a single best fit against single values. The ensemble of optimal solutions shows a wide spread of parameter settings within the physically reasonable range. This implies that the common approach of a single best optimized parameter set is subject to over-fitting and may significantly differ from other equally plausible solutions, meaning that they are not representative by default. Furthermore, our results show that the constraint of plausible parameters is only marginally linked to the sensitivity, with very sensitive parameters also taking multiple optimal values. This implies that keeping these parameters constant in space and time increases the model uncertainty. The overall model uncertainty (not accounting for uncertainties related to meteorological forcing data) is in the range of 1000 kg m⁻² per summer season on the same glacier, and increases when applied to the other glacier. The model performance is worse when applied to another glacier, but is of the same order of magnitude as for the temporal transfer, suggesting the model can be applied, within its uncertainty, to other glaciers with similar climatic settings.

Parameter uncertainty is connected with uncertainty in the energy balance components, which, in the cases studied here, reached 30 % averaged over the glacier and 50 % at individual stake locations. In our study the most uncertain energy balance components are the net short-wave radiation and the turbulent fluxes, reasserting the findings of other studies (e.g. Van De Wal et al., 1992; Klok and Oerlemans, 2002) that indicate the snow and ice albedo representation is the most crucial parameter on mid-latitude glaciers for the summer mass balance.

Overall, the findings of this study highlight that understanding the sensitivity and uncertainty of surface energy and mass balance models is complex, and simplistic assessments, in particular single best guess approaches, of model performance are likely to overstate model capabilities. Further studies such as this, incorporating more models, glaciers and years would help to constrain the degree to which results from such models can be considered reliable for regional applications and for projections of glacier mass balance.

Appendix A: Model description

The mass and energy balance model used here consists of coupled surface and subsurface components. The model computes mass balance as the sum of solid precipitation, surface deposition, internal accumulation (refreezing of liquid water in snow), change in englacial liquid water storage, subsurface and surface melt and sublimation. This approach is based on the surface energy balance of a glacier in the following form:

$$\begin{array}{}\text{(A1)}& Q_{\text{M}}+Q_{\text{ice}}={\text{SW}}_{\text{net}}+{\text{LW}}_{\text{net}}+Q_{\text{S}}+Q_{\text{L}}+Q_{\text{G}}+Q_{\text{P}},\end{array}$$

where SW_net is net short-wave radiation, LW_net is the sum of incoming and outgoing long-wave radiation at the glacier surface; Q_S and Q_L are the turbulent fluxes of sensible and latent heat, respectively; Q_G is the subsurface energy flux comprised of Q_C; the conductive heat flux in the subsurface; and Q_PS is the energy flux from short-wave radiation penetrating into the subsurface, and finally, Q_P is the heat flux from precipitation. The sum of these fluxes yields a residual flux F which, if the glacier surface temperature (TS) reaches 273.15 K, represents the latent energy for melting. If TS is below 273.15 K, energy conservation is achieved by solving TS to balance the fluxes (e.q. Mölg et al., 2009). The model is fully described in the previously mentioned publications and briefly below.

A1 Long-wave radiation

The calculation of the incoming long-wave radiation is based on Stefan–Boltzmann law (Mölg et al., 2009; Klok and Oerlemans, 2002; Konzelmann et al., 1994):

$$\begin{array}{}\text{(A2)}& {\text{LW}}_{\text{in}}=\mathit{\sigma }\mathit{ϵ}{T}_{\text{a}}^{\mathrm{4}}\end{array}$$

with σ being the Stefan–Boltzmann constant and ϵ the emissivity:

$$\begin{array}{}\text{(A3)}& \mathit{ϵ}={\mathit{ϵ}}_{\text{cs}}(\mathrm{1}-n^p)+{\mathit{ϵ}}_{\text{cl}}{n}^p,\end{array}$$

where “cs” and “cl” are the clear-sky and cloud emissivity, respectively, n is the cloud cover fraction calculated in the solar module as n_eff and p is an exponent related to the importance of cloud emissivity (Greuell et al., 1997). The cloud emissivity is computed using

$$\begin{array}{}\text{(A4)}& {\mathit{ϵ}}_{\text{cl}}=\mathrm{0.23}+b{\left({\displaystyle \frac{e_{\text{a}}}{T_{\text{a}}}}\right)}^{\mathrm{1}/\mathrm{8}},\end{array}$$

with e_a as the atmospheric vapour pressure. The three parameters ϵ_cs, p and b were optimized (using a 5000 member Monte Carlo) to reproduce the measured long-wave radiation. First the runs within 10 % of the best run, in respect to a weighted average of BIAS and RMSD, between the simulated and the measured incoming long-wave radiation at the HEF Station were determined. The run of this ensemble with the lowest RMSD/BIAS on LGF was taken as the best compromise solution. The parameters are fixed within the model for the whole study period and are based on three summers of data at HEF and one and a half at LGF (therefore a larger impact of the longer data at HEF on the optimization). The trade-off values are taken to be applicable on both glaciers with the final values of b=0.515, n=1.95 and ϵ_cs=0.994. These setting results in an hourly RMSD below 31/37 W m⁻² for HEF/LGF and no bias, not far from the optimal setting for either glacier with 30/36 W m⁻².

The outgoing long-wave radiation follows the Stefan–Boltzmann law as in Eq. (A2), with T the glacier surface temperature and the emissivity of ice ϵ_i is assumed 1.

A2 Convective fluxes

The latent heat flux (Q_L) and the sensible (Q_S) are computed similar to Mölg and Hardy (2004). The calculations are based on the Monin–Obhukov similarity theory (Garratt, 1992).

$$\begin{array}{}\text{(A5)}& Q_{\text{L}}=\mathrm{0.623}{L}_{\text{v}}{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{\mathrm{1}}{p_{\mathrm{0}}}}{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}\mathit{\nu }(e_{\text{a}}-E_{\text{s}})}{\ln \frac{z_{\text{m}}}{z_{{\mathrm{0}}_{\text{m}}}}\ln \frac{z_{\mathit{\nu }}}{z_{\mathrm{0}\mathit{\nu }}}}},\end{array}$$

with L_v being the enthalpy of vaporization (2.514 MJ kg⁻¹), ρ₀ the air density at mean sea level (1.29 kg m⁻³), p₀ is 1013 hPa, κ the van Karman constant (0.4), e_a is the water vapour pressure in air and E_s the surface value, respectively. $z_{{\mathrm{0}}_{\text{m}}}$ and z_0ν are the momentum and scalar roughness length of water vapour. z_m and z_v is the height above ground where the wind speed and the water vapour (e_a) is measured and calculated. The sensible heat flux

$$\begin{array}{}\text{(A6)}& Q_{\text{S}}=c_{\text{p}}{\mathit{\rho }}_{\mathrm{0}}{\displaystyle \frac{p}{p_{\mathrm{0}}}}{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}\mathit{\nu }(T_{\text{a}}-T_{\text{s}})}{\ln \frac{z_{\text{m}}}{z_{{\mathrm{0}}_{\text{m}}}}\ln \frac{z_{\text{h}}}{z_{{\mathrm{0}}_{\text{h}}}}}},\end{array}$$

is computed with c_p, the specific heat of air at constant pressure, T_a, T_s the air and surface temperature and z_h the scalar roughness length for temperature. The roughness lengths (z_j) are a model free parameters in this study. The model distinguishes three different roughness lengths depending on the glacier surface: fresh snow, firn and ice. For a stable stratified atmosphere, a stability correction based on Phi functions is applied (Mölg and Hardy, 2004).

A3 Surface albedo and the albedo module

The albedo parameterization is based on Oerlemans and Knapp (1998). It computes the broadband albedo for each grid cell, based on the ice and snow albedo and the depth of the snow pack:

$$\begin{array}{}\text{(A7)}& \mathit{\alpha }={\mathit{\alpha }}_{\text{snow}}+({\mathit{\alpha }}_{\text{ice}}-{\mathit{\alpha }}_{\text{snow}})\cdot \exp \left({\displaystyle \frac{-d}{d^{*}}}\right),\end{array}$$

α_ice is a model free parameter d is the snow depth, and d^* is the characteristic scale for the snow depth and a free parameter (Oerlemans and Knapp, 1998). The relation for the snow albedo (α_snow) is

$$\begin{array}{}\text{(A8)}& {\mathit{\alpha }}_{\text{snow}}={\mathit{\alpha }}_{\text{firn}}+({\mathit{\alpha }}_{\text{freshsnow}}-{\mathit{\alpha }}_{\text{firn}})\cdot \exp \left({\displaystyle \frac{-t}{t^{*}}}\right),\end{array}$$

with α_firn, α_freshsnow and t^* as model free parameters subject to optimization. The albedo module (t^*) is a characteristic time scale in days (Klok and Oerlemans, 2002) and t the time since the last snowfall event (>0.1 cm fresh snow).

A4 Surface temperature and ground energy flux

The conductive heat flux (Q_C) and the energy flux from penetrating short-wave radiation (Q_PS) determine the ground heat flux (Q_G) of the energy balance (Eq. 1). The model solves the thermodynamic energy equation for a multi-layer grid with a fixed bottom temperature (15 Layers, 0.1 m steps in the first metre, gradually increasing to a total depth of 7 m). The bottom temperature is a model free parameter. Q_C is computed from the temperature difference between the surface and the first layer.

The calculation of the penetration of short-wave radiation is based on Bintanja and Van Den Broeke (1995). A constant fraction (1−ζ_i) of the net-short-wave radiation is penetrating the surface and the intensity is exponentially decreasing with depth. The optimization and sensitivity analysis in this study uses four parameters with the extinction coefficient and the absorbed fraction (ζ_i) for snow and ice.

A5 Surface accumulation and precipitation

The surface accumulation is directly related to the precipitation. The model has two threshold values for all liquid and all solid precipitation (Mölg et al., 2012). In between these thresholds, the portion increases linearly. The temperature threshold and the density of solid precipitation are subject of the sensitivity analysis and optimization.

A6 Solar module and solar module sensitivity

The parameterization of the short-wave radiation is based on the calculation of the cloudiness, in the form of the effective cloud cover fraction n_eff:

$$\begin{array}{}\text{(A9)}& n_{\text{eff}}={\displaystyle \frac{\mathrm{1}-{\text{SW}}_{\text{mea}}/(D_{\text{cs}}+S_{\text{cs}})}{k}},\end{array}$$

with SW_mea being the measured short-wave and D_cs, S_cs the calculated diffuse and direct radiation under clear sky conditions. The parameter k determines at which fraction of the clear sky value full cloudiness is achieved; i.e. all incoming radiation is diffuse. It is important to note, we allow n_eff>1 if such low radiation was measured. The influence of k on the model output was investigated (Appendix A6). The calculation of the clear sky values is described in Mölg et al. (2009). The diffusion portion of radiation under clear sky conditions was determined using a manual selection of clear sky days. The values varied between the snow free (K_dif=0.51) and snow covered days (K_dif>0.65). For the calculations an averaged value of 0.6 was used. As K_dif is a fixed glacier wide value, while snow cover might vary, a modulation depending on the conditions at the weather station is not possible. The applicability of K_dif as a single value might need to be reevaluated for other models, applications and research questions.

The calculation of the incoming short-wave radiation on every point of the glacier is based on the assumption of homogeneous cloudiness (n_eff). It is a reversing of Eq. (A9)

$$\begin{array}{ll}\text{(A10)}& {\displaystyle }& {\displaystyle }{\text{SW}}_{\text{diff}}={\displaystyle }& {\displaystyle }(D_{\text{cs}}+S_{\text{cs}})\cdot (\mathrm{1}-n_{\text{eff}}\cdot k)\left(\right(\mathrm{1}-p_{\text{diff}})\cdot n_{\text{eff}}+p_{\text{diff}}),\end{array}$$

with SW_diff being the calculated diffuse radiation and p_diff the portion of diffuse radiation under clear sky conditions. p_diff is calculated as the ratio of the clear sky diffuse and total radiation. It was 0.084 and 0.085 for the two glaciers and set to 8.5 % (for both to have a common value). Compared to previous works using the solar module, we changed the increase of diffuse radiation. Instead of a linear increase of diffuse radiation, the portion of diffuse radiation is linearly increasing with increased cloudiness. This is a basic parameterization and reproduces the measured radiation fully. Via n_eff, k is determining the ratio of direct and diffuse radiation. This could alter the energy balance. The direct radiation is calculated analogously and corrected for slope and aspect.

The calculation of solar radiation incorporates the free parameter k, which determines at which fraction of the total possible global radiation everything is considered as diffuse radiation. The parameter k varies with latitude (Hastenrath, 1984) and is not constant in time either; therefore the effective cloud cover incorporates some of its variability and is not exactly the cloudiness (Mölg et al., 2009). With the new parameterization (Eq. A10), the global solar radiation at the weather station can be fully reproduced so k cannot be optimized. But it determines the portions of direct and diffuse radiation, which may have a significant influence on the energy and mass balance. Therefore, an additional GSA was performed with the parameter k as the 24th model free parameter. Based on the values for the tropics 0.65 (Mölg et al., 2009) and the Arctic with ≈0.85 (Hock and Holmgren, 2005) it was varied in this range for the sensitivity analysis. Its maximum sensitivity index over all seven investigated stakes in the GSA was $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$, which is 1 order of magnitude lower than the threshold for our sensitive parameters. Therefore, the choice of k within the given range is not influential on the simulation of the mass balance on the glacier. The model albedo does not vary between direct and diffuse radiation, so it only influences the total amount of radiation at less or more shaded areas than the weather station.

Furthermore, the change in the calculation of direct and diffuse components from linear with cloudiness to a linear increase of the fraction are better suited to represent the site radiation. This is in agreement with the measured radiation by Hock and Holmgren (2005) on the Arctic glacier, Storglaciären (Fig. S5). The slightly higher starting value (p_diff) is due to larger portion of diffuse radiation under clear sky conditions in the Arctic than in the mid-latitudes, and the higher final value is due to a smaller k in this study (0.8) compared to around 0.85 in the Arctic. The influence of this change in parameterization is probably also rather small, as the model is not sensitive to changes in the relative fractions of diffuse and direct radiation on the chosen glaciers and stake locations.