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:

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 1The 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.

Download Print Version | Download XLSX

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

Download Print Version | Download XLSX

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 3The 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.

Download

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.1 Global sensitivity analysis

Figure 4The 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.

Download

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 3Five 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.

Download Print Version | Download XLSX

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 5Each 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.

Download

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.

Download

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 7The 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).

Download

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 8The 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.

Download

Table 4The 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.

Download Print Version | Download XLSX

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.