Sensitivity of the Greenland mass and energy balance to uncertainties in key model parameters

. We investigate the sensitivity of a glacier surface mass and energy balance model using a variance based analysis, for two distinct periods of the last glacial cycle: present day climate and the Last Glacial Maximum (LGM). The results can be summarized in three major ﬁndings: The sensitivity towards individual model parameters and parameterizations is neither invariant in space nor in time. The model is most sensitive to uncertainty related to down-welling long-wave radiation. Turbulent latent heat ﬂux has a sizable contribution to the surface mass balance uncertainty in central Greenland today and over the entire 5 ice sheet during the cold climate of the LGM, in spite of its low impact on the overall surface mass balance of the Greenland ice sheet in modern climate. We conclude that quantifying the model sensitivity is very helpful for the subsequent tuning of free model parameters, because it clariﬁes the relative importance of individual parameters and highlights interactions between them that need to be considered. radiation is about 10% higher for the LGM, though the seasonal cycle deviates drastically from PD. The surface mass balance model uses a static topography, which is based on ETOPO (Amante and Eakins, 2009) for PD and on ICE-6G (Peltier et al., 2015) for the LGM (ﬁg. 1). There is an interpolation artifact in the PD simulation in the far north, but as no ice is present in this region of Greenland it does not inﬂuence the analysis. BESSI was run for 500 years with the same forcing data back and forth to account for the long response time of the ﬁrn cover. 400 years was shown to be sufﬁcient 25 to develop a dynamically and thermodynamically stable ﬁrn cover, even in the regions of very low accumulation. The analysis shown throughout this study is entirely based on an assessment of the last 100 years of every simulation. The model has nine free parameters which are listed together with their plausible range in Table 1. The initial ensemble was generated with using a Sobol sequence which consisted of 2000 ⇥ 9 members for PD and 1000 ⇥ 9 for the LGM. This sequence spans a 9-dimensional unit hypercube. It is split into two subsets A and B each consisting of one half of the initial sequence 30 ( 1000 ⇥ 9 ). For computing both sensitivity indices the estimator from Sobol et al. (2007) was used. It creates an additional set of matrices B i A , which are based on the matrix B where the values for parameter X i are replaced with those from subset


Introduction
Of the many challenges to accurately simulate past variations in the volume of the Greenland ice sheet (GrIS) and to project its future contribution to sea level rise, recent studies agree that the uncertainty associated with surface mass balance (SMB) is among the most important (Plach et al., 2019;Aschwanden et al., 2019).
Models to calculate SMB spread a whole range of complexities from empirical index models that only account for air tem- 15 perature (Ohmura, 2001;Zemp et al., 2019), or temperature and solar radiation (Robinson et al., 2011), to coupled atmospheresnow models that simulate the snowpack in multiple layers and a full representation of the atmospheric circulation, based on physical first principles (Lehning et al., 2002;Fettweis, 2007;Noël et al., 2018). On this spectrum, the empirical models perform well for the observational period and when the temperature sensitivity of the SMB is well known, but are hard to constrain for temporal and spatial variations and become unreliable for conditions outside their relatively narrow tuning interval (van de 20 Berg et al., 2011;Plach et al., 2019). This is unfortunate because while their low computational requirements make these models attractive for the long integration times that are needed to simulate continental ice sheets, their shortcomings severely limit Variations of what? I think there is something missing here.
Rephrase to avoid starting with "this", clearly state subject instead. the usefulness of these results. On the other hand, detailed snow models and especially those coupled with regional atmosphere models are too expensive to run for long periods of time. This situation motivated the development of models that balance the defensible representation of the relevant physical processes with computational efficiency (Krapp et al., 2017;Krebs-Kanzow et al., 2018;Born et al., 2019).
In this study, we use the BErgen Snow SImulator (BESSI), a model that is designed to include all relevant physical mech-5 anisms with reasonable detail but specifically prepared for long integration times by reducing its computational requirements and by strictly conserving mass and energy . Adding to the original model version, we now include parameterizations for snow aging and turbulent latent heat flux while limiting the model domain to only Greenland at a resolution of 10 km. Three different snow aging parametrizations are used based on Oerlemans and Knapp (1998), Aoki et al. (2003), and Bougamont et al. (2005). 10 Multiple studies have investigated the impact of the turbulent latent heat flux on the surface mass balance (Box and Steffen, 2001;Box et al., 2004;Van Den Broeke et al., 2008;Cullen et al., 2014;Noël et al., 2018). They find a relatively small impact of the vapor fluxes on the Greenland ice sheet total mass balance (⇡ 5 Gta 1 , (Cullen et al., 2014)), but its local impact can be up to 20% of the annual accumulation (Box et al., 2004). It is difficult to assess the importance of this component for different climatic settings, as Box and Steffen (2001) have shown that the choice of the calculation method impacts the results greatly 15 in regions of low flux like the dry interior zone of Greenland. This is exacerbated by the fact that turbulent latent heat fluxes are mostly negative in winter under the present climate conditions, and positive in summer so that the sign of the net flux may change with changing climate. During the colder climate of the glacial a much larger impact of the turbulent latent heat flux can be assumed, similarly to the much greater importance it currently has in Antarctica (e.g., Gallet et al., 2014;Van Wessem et al., 2018). A parametrization based on the bulk-method by Rolstad and Oerlemans (2005) has been added to BESSI to simulate 20 the turbulent latent heat flux.
To assess the sensitivity of the new parameterizations and that of BESSI overall, we employ a variance based approach (Saltelli et al., 2000(Saltelli et al., , 2006(Saltelli et al., , 2010Sauter and Obleitner, 2015), that has previously been used to quantify the sensitivity of a glacier model (Zolles et al., 2019). Following our model's design goal to be used over time scales of glacial cycles and accounting for potentially different sensitivities under different climate boundary conditions, we analyze two large ensembles 25 with a total of 16,500 simulations, for present day (PD) climate and that of the Last Glacial Maximum (LGM). The result is rich information on what parameters and parameterizations have the largest impact on the model's performance, how this sensitivity varies in different regions of Greenland and over time. This knowledge enables a balanced calibration of the model as it prevents parameters to be specifically adjusted to the study location and time (Beven, 1989).
The revised model is described in section 2. The sensitivity of multiple model output variables to model parameters, including 30 those for the new parameterizations for turbulent latent heat flux and snow aging, is presented in section 3. After that, we discuss our findings in section 4 and conclude in section 5.

Model description and study setup
The study uses the efficient mass and energy balance model BESSI, which is designed to simulate the mass balance over long time scales . The energy exchange of the snow with the atmosphere is altered in the model version used here, while the subsurface and internal processes are unchanged from the published version. It was enhanced to include the turbulent latent heat flux and multiple more complex snow albedo schemes were added. The model description given in the following 5 will focus entirely on the snow surface interaction with the atmosphere. For the numerical description and subsurface processes like firnification and heat conduction see Born et al. (2019).
The model setup used here has a 10 km horizontal grid over the domain of Greenland. The vertical is discretized using a mass following grid with up to 15 layers. The mass of each layer ranges 100 -500 kgm 2 . Simulations require daily input of air temperature, total precipitation, solar radiation and its reference height. Humidity is an optional input which is required 10 if the turbulent latent heat flux is computed. The air temperature is the only meteorological input which is downscaled to the actual model topography. The output written by the model can be adjusted by the user ranging from daily over monthly to annual values. Output variables include: surface mass balance, melt of snow, melt of ice, runoff, refreezing, albedo, turbulent latent heat flux, a mask containing snow, land, ice and water as well as the 3-dimensional grid values for snow mass, snow density, snow temperature and liquid water mass. 15

Surface energy fluxes
The energy exchange between the surface and the atmosphere contains five different processes of which two also imply a change in mass: The short-wave radiation (Q SW ), the long-wave/thermal radiation (Q LW ), the turbulent sensible heat flux (Q SH ), the turbulent latent heat flux (Q L ) and the heat flux associated with precipitation (Q P ).
The total surface flux can be expressed as: 20 c i m s,1 @T @t surf ace where the left hand side denotes the resulting temperature change of the snow/ice (Q i ) and the available energy for melting (Q M ) if the melting point is reached. Due to the implicit scheme the model uses, no melt is calculated at first, but only energy fluxes and temperatures (even above 273 K). The actual melt is calculated explicitly each time step for the excess heat above the melting point. The mass flux of the water vapor is also calculated explicitly. The energy input to the surface from solar radiation is calculated by using a broadband albedo value: where ↵ denotes the surface albedo, either ↵ s for snow or ↵ i for ice and SW in the incoming short wave radiation at surface height. The albedo value is assumed constant with respect to the solar incidence angle, but undergoes temporal and spatial 30 variations depending on surface properties. We implement four albedo parametrizations of different complexity to simulate the snow albedo. They all have a common maximum albedo value for fresh snow (↵ fs ), minimum value for aged snow (firn ↵ fi ), and ice albedo (↵ i ), but vary in how they calculate the aging.
Constant: This simple parametrization only uses constant values for dry snow (↵ s = ↵ fs T S < 273K), wet snow (↵ s = ↵ fi T s = 273K), and ice (↵ i ). This parameterization has been used before in BESSI . 5 Oerlemans and Knapp (1998): This parametrization assumes an exponential decay with time of the fresh snow albedo to a final value of old snow albedo (Oerlemans and Knapp, 1998): where t fs denotes the last day of snowfall, t the current day (time step) with t ⇤ as the characteristic time in days. This or similar parametrizations are usually optimized for the decay rate (t ⇤ ) using observations of albedo or mass balance (e.g., 10 Oerlemans and Knapp, 1998; Klok and Oerlemans, 2004;Bougamont et al., 2005). The very fast decay at the melting point was chosen to account for our very large upper grid box which would likely already be wet if the real surface was resolved.
Equation 3 does not consider shallow snow packs, where underlying ice or dirty firn albedo may reduce the albedo. Bougamont et al. (2005) modified the parametrization by Oerlemans and Knapp (1998) specifically for the Greenland ice sheet by introducing a snow temperature-dependent decay rate: This results in the same t ⇤ of 30 d as Oerlemans and Knapp (1998) for 273 K. This parameterization furthermore introduces a wetness-dependent albedo decay, which assumes a thin layer of water at the surface according to Here w surf denotes the thickness of the water layer and w ⇤ is a characteristic water layer thickness. Since BESSI does not 20 explicitly simulate water at the surface, this was adapted for this model by using a simple linear parametrization. The decay rate increases depending on the liquid water content ⇣ (see 2.2 for details about the liquid water content): Aoki et al. (2003): The final albedo parametrization available in the model uses both a fixed and a temperature dependent decay rate, rather than keeping it constant or relating it to the time since the last snowfall.

25
↵ s (t) = min{↵ s (t 1) ((T S 273 K) · k + c),↵ fi ,↵ s (t 1)} where t and t 1 are the current and the last time step, ↵ s the snow albedo and k = 1.35 · 10 3 K 1 and c = 0.0278 two empirically based constants. The values are based on averaged values from Aoki et al. (2003) for different spectral bands. To account for a faster decay in a wet snow pack it is linearly decreased based on the liquid water content of the topmost layer: If the layer is fully saturated with water the snow albedo instantly drops to its minimum value. In this parametrization 5 the albedo increases if new snowfall occurs, but instead of resetting it to the fresh snow value, it is incrementally increased depending on the amount of fresh snow: with d as the snow depth. The characteristic snow depth d ⇤ is at 3 cm (Oerlemans and Knapp, 1998).

10
The long-wave radiation is a simple parametrization based on the Stefan-Boltzmann relation: where is the Stefan-Boltzmann constant, T air and T s are the 2 m air and snow surface temperature, respectively. The emissivity of snow/ice ✏ s is constant at 0.98. The incoming long-wave radiation is only depending on the actual air temperature and the atmospheric emissivity ✏ atm , as the only free model parameter. The lack of confidence in cloud cover of climate models 15 in particular during the last glacial cycle, lead to this decision. Though more complex empirical relations exist (e.g., Listion and Elder, 2006), their applicability for other time scales is questionable. We therefore refrain from using these empirical relationships despite the importance of cloud cover, moisture and aerosols. ✏ atm varied over a broad range of 0.6-0.9 following the previous configurations of this or similar models (Greuell and Konzelmann, 1994;Busetto et al., 2013;Born et al., 2019).
Values spanning this entire range may in reality occur simultaneously in different regions of the Greenland ice sheet, but in the 20 current configuration there is only a single atmospheric emissivity value over the entire ice-sheet.

Turbulent sensible heat flux
The calculation of the turbulent latent heat flux is based on a bulk method (Braithwaite, 2009) which was applied previously on Greenland as residual method (Rolstad and Oerlemans, 2005). This method assumes a constant turbulent exchange coefficient (C h ) for sensible heat over time and space. The only dependency in the previously published parametrization is on the local 25 wind speed u and air temperature T air : previous where ⇢ air is the density of air and c p the heat capacity of air. Since BESSI does not use wind speed as input field, we simplify the equation with a single turbulent heat exchange coefficient D SH . This is a free parameter and subject of the sensitivity analysis. The values given in Table 1 assume an average wind speed of 5 ms 1 if compared to the reported values by Braithwaite (2009). The variation in parameter D SH therefore accounts for both the variability in average wind speed and the efficiency of the exchange C h . The previous version of BESSI did not include turbulent latent heat flux. Here it is an optional part of the setup. The implementation is analog to the turbulent sensible heat flux (eq. 11): where D LH is the turbulent latent heat exchange coefficient, and e is the water vapor pressure. The parametrization is based 10 on the bulk formulation of Rolstad and Oerlemans (2005) with the latent heat of vaporization L v and the air pressure p. The latter is calculated from the standard pressure at sea level for each grid point. While Rolstad and Oerlemans (2005) assume the same exchange coefficient C h for vapor and sensible heat, it has been shown previously that the roughness lengths and the exchange coefficient for momentum and vapor are not necessarily equal (Greuell and Smeets, 2001). Nevertheless, the parameters D LH and D SH are inherently connected by the surface structure (snow/ice) and the wind speed. To account for 15 these points the fesetup uses two free parameters determining D LH , the turbulent exchange coefficient for sensible heat D s h and r lh/sh which are defined by In the setup of this study there are three parameters determining the turbulent latent heat flux. The ratio (r lh/sh ) accounts for different exchange rates for water vapor and sensible heat. D SH the absolute exchange strength (roughness, wind, stability).

20
The parameter ( QL ) switches the subroutine on and off.

Precipitation heat flux
The heat supplied by precipitation depends on the temperature. An atmospheric temperature of 273 K is the limit of solid precipitation. In the case of snow the solid mass of the topmost grid cell increases by the amount of snow and the heat added to this box reads: while rain is added as liquid water mass to the same cell where P is the amount of precipitation, ⇢ w and c w the density and heat capacity of water. If the rain freezes or is percolating further down, this is calculated during the balance calculation, and as outlined in the next section.

Subsurface percolation and refreezing
This is the only part of the subsurface routine which is explained here. Full details are available in Born et al. (2019). The water holding capacity of each layer determines the percolation. The maximum liquid water content ⇣ max is the parameter subject to 5 the sensitivity analysis: where m w , m s and ⇢ w , ⇢ s , ⇢ i are the mass and density of water, snow and ice, respectively. If the liquid water content ⇣ exceeds ⇣ max the liquid water mass percolates to the next cell or is treated as runoff when it leaves the bottom cell. The global sensitivity analysis is variance based method that allows for a full assessment of the model sensitivity over the entire parameter space. In comparison to other methods, all parameters are changing at the same time. The method has been applied previously to snow (Sauter and Obleitner, 2015) and recently to alpine glaciers (Zolles et al., 2019). It is based on algorithms developed by Saltelli et al. (2000Saltelli et al. ( , 2006Saltelli et al. ( , 2010, utilizing the setup of the ensemble hypercube using (Sobol 15 et al., 2007). The probabilistic framework provides an estimate of the sensitivity of the model output to the individual inputs, including parameters and data. It is independent of model calibration and tuning. The model output Y is a function of the input parameters X i : Y = f (X 1 ,X 2 ,..,X n ). There are two normalized values that quantify the model sensitivity for each input parameter X i : the first or main order sensitivity index S Xi and the total sensitivity index S T i of parameter X i . The first order index denotes the sensitivity of the model towards the parameter X i only, while the latter includes all the interactions of X i : Are "full" and "entire" redundant?

Global sensitivity analysis
Awkward sentence, not sure what it means. Rephrase or delete.

Run-away sentence
where E is the expectation value of a given observable such as the SMB. V Y is the total variance of the given variable and V Xi the variance that only depends on the input paramter X i . X i denotes the whole parameter space excluding any variation 5 in X i . The first order index calculates the mean model output (E X i (Y |X i )) for each representation of X i and then assesses the sensitivity by calculating the variance for all values of X i .
The total index can be compared to the local sensitivity index that is often determined around the optimal model setting.
V Xi (Y ) varies the parameter X i along its dimension, but it is computed for all possible points of the parameter space instead of the optimal one. For the detailed algorithm refer to Saltelli et al. (2010). As under-sampling is assumed by the relatively low PD. The annual mean solar radiation is about 10% higher for the LGM, though the seasonal cycle deviates drastically from PD.
The surface mass balance model uses a static topography, which is based on ETOPO (Amante and Eakins, 2009) for PD and on ICE-6G (Peltier et al., 2015) for the LGM ( fig. 1). There is an interpolation artifact in the PD simulation in the far north, but as no ice is present in this region of Greenland it does not influence the analysis. BESSI was run for 500 years with the same forcing data back and forth to account for the long response time of the firn cover. 400 years was shown to be sufficient 25 to develop a dynamically and thermodynamically stable firn cover, even in the regions of very low accumulation. The analysis shown throughout this study is entirely based on an assessment of the last 100 years of every simulation.
The model has nine free parameters which are listed together with their plausible range in Table 1. The initial ensemble was generated with using a Sobol sequence which consisted of 2000⇥9 members for PD and 1000⇥9 for the LGM. This sequence spans a 9-dimensional unit hypercube. It is split into two subsets A and B each consisting of one half of the initial sequence I think I know what you mean, but it should be rephrased. Do you mean the firn pack is in equilibrium with the climate (or something like that)? Not a good lead-in sentence. Your model has more than 9 free parameters (e.g. $c_w$), but you select the nine parameters that you think are the most relevant. Elaborate on our choices instead of just referring to the table. A. The matrices A, B and B i A are then used to estimate the model sensitivity. The whole ensemble consists of N ⇤ (2 + k) members, with N being the base sample (1000 in the case of PD) and k the amount of parameters. A detailed description of the algorithm can be found in (Sobol et al., 2007) and (Saltelli et al., 2010).
The full ensemble for present day climate has 11000 members, that for LGM climate 5500. The normalized parameter dimensions are linearly transformed to the intervals given in table 1, with the exception of the latent heat flux switch and the 5 albedo module. These two parameters have two, respectively four discrete values, and the parameter space is split equally between them. The model simulations are carried out with the generated parameter matrix. We are using bootstrapping to estimate the sensitivity indices. For each bootstrap the equations 17 and 18 are evaluated. Finally, we report the mean sensitivity indices and their standard deviation. The results were checked for consistency ( P S Xi  1, S Xi  S T i ). The ensemble size used during the LGM is on the absolute lower limit of applicability for GSA, as the standard deviation of the sensitivity 10 indices is large (not shown). It works fine for the output quantities analyzed, some with increased uncertainty associated to the absolute values, but fails for the 10 m firn temperature for example. Due to the larger ensemble the confidence in PD ensemble is higher, but not by a large enough margin to justify the additional computation time relative to the 5500 members of the LGM ensemble. BESSI is a complex model with all parameters interacting and S Xi provides very little information. Therefore, the results mainly focus on the total sensitivity index. 15 The GSA was computed for five different outputs: albedo, vapor flux, snow melt, surface mass balance and the average surface temperature. These are based on average yearly sums or averages over 100 years. Surface temperature results are rather uncertain and only tendencies can be extracted as it is largely influenced by the annual cycle and fresh snow reseeding on an ice surface.
The global sensitivity analysis of the SMB provides the main order effect (circle) and the total effect (triangle). The two indices are displayed for all 9 parameters over the different elevation bands ranging from 0-1000, 1000-2000, 2000-2500, 2500-3000 and above 3000 m. The symbol represents the mean value of the sensitivity index with the bars as ± 1 . The elevation bands are displayed over Greenland on the right, the analysis is only done for cells where ice is present.

Global sensitivity analysis -GSA
The main focus of the results will be on the surface mass balance and discussion is limited to the total sensitivity index (S T i ), due to the low information that can be extracted from S Xi in complex models. The sensitivity of the SMB changes for different elevation classes for the present day ice sheet is shown in Fig. 2. In the region of lowest elevation the largest sensitivity is Being more specific makes the language flow better, e.g.: At elevations below 1000m, the largest...
Remove total effects and consider collapsing into one, colored by elevation band.

for present day
The GSA was also calculated for each grid cell individually. The general trends are similar to the elevation average, but the structure is much more complex. The global sensitivity maps for all parameters are displayed in Figure 3. The model is sensitive to ✏ atm over the entire ice sheet. Its impact is only reduced at the west coast glaciers where D SH is more sensitive and the north west where the choice of albedo module is important. The model is sensitive to QL in the interior of Greenland, and can be considered a sensitive parameter for most of the ice sheet apart from the regions of very high melt at low elevation.

5
In the interior of very high elevation the most sensitive of the three snow albedo related parameters is the one for fresh snow ↵ fs , while at lower elevation firn snow albedo ↵ fi becomes more important and also the chosen albedo parameterizations ↵ becomes more important, in particular in the north-east. The SMB is sensitive to D SH at very low elevations in the west and on the ice-sheet above 1500 m, only at the top of the ice sheet its influence is reduced. The ice albedo ↵ i plays a minor role in the north, but is generally of very low impact. Additionally, ⇣ max as well as r lh/sh are of minor importance.

10
The dominance of ✏ atm is a result of the change in total heat flux associated with its parameter uncertainty. At an annual average of -10°C, a change of ✏ atm from 0.6 to 0.9 increases the heat flux by 80 Wm 2 , while a change of albedo from 0.9 to 0.6 does only increase the energy input by 40 Wm 2 for typical values of solar radiation. Similarly, the sensible heat flux is smaller than the other heat fluxes over most of the ice sheet (monthly averages from -20 to +50 Wm 2 with D SH = 12 Wm 2 K 1 ) so even a doubling will not be larger than the change of atmospheric emissivity on an annual basis. The relatively low impact 15 of the ice albedo ↵ i is due to the small exposure time and the small parameter range. An ice albedo change of 0.3 to 0.4 will at most result in an annual average energy change of 10-20 Wm 2 , but the ice is never exposed for the whole season. Rather the date of ice exposure, which is a result of the energy fluxes prior to its exposure, is much more important as this is associated with an surface albedo change from 0.55-0.9 to 0.3-0.4. The SMB in regions above 1000 m is clearly sensitive to the snow albedo, which in turn is a function of snowfall, snow temperature and the chosen albedo parametrization. Each albedo model 20 treats these processes separately. While the basic one only distinguishes between dry and wet snow, the others account for snow aging ranging from time over time-and-temperature to time-temperature-wetness dependency. The choice of albedo model is not important in the interior of Greenland, as temperatures are low, albedo decay is slow and the snow does not get wet. At slightly lower elevations there is an interplay of the fresh snow ↵ fs and firn albedo ↵ fi and the chosen decay parameterization ↵ , with a bit larger impact of the fresh snow albedo, as snowfall is not frequent and decay rates at medium temperatures are 25 of the order of a few weeks (Sec. 2.1.1). An exception is found in the north east, which is characterized by the driest climate on Greenland. The less frequent precipitation and therefore albedo reseeding explains a larger dependency on the decay rate and the albedo module. This is furthermore an area higher up on the Greenland ice sheet where QL is less important. Due to the low amount of precipitation the SMB is quite sensitive to changes in the surface energy balance. Extreme ensemble members (low albedo, high emissivity) lead to a melting state, negative ensemble members result in very large changes in SMB relative 30 to the small changes associated with condensation and sublimation. In this dry region a very large atmospheric emissivity is not realistic. The large impact of D SH at the western coast is due to the large air-surface temperature difference, while in the interior Greenland the atmosphere is more in balance with the snow surface.
The model" is the subject in the first part of the sentence thus you need to introduce the new subject $\Chi_{QL}$ Be specific. Like "at elevations < 500m" or "at very low (<500m) elevations" air temperature resetting Be consistent within a sentence, don't use the math symbol in one instance and the descriptive words in another. Rephrase sentence and avoid "negative ensemble member". Explain what is shown in the figure before delving into the discussion.

Clarify
heat flux switch QL is as important as the atmospheric emissivity ✏ atm already above 1000 m and the fresh snow albedo ↵ fs is almost as important above 2000 m. The ice sheet integrated SMB shows sensitivity to atmospheric emissivity, the turbulent heat flux coefficient, the latent heat flux switch and the snow albedo. The liquid water content ⇣ max , the ice albedo ↵ i and the ratio of latent and sensible heat exchange coefficient r lh/sh do not impact the SMB in either of two climate states. On the local scale QL followed by ✏ atm and ↵ fs are the main sensitivity components (fig. 4). The sensible heat flux exchange 5 coefficient D SH is important along the margin with the largest impact in the south-east. The firn albedo and the albedo module are sensitive parameters along the margin, with a bit higher impact of the module along the western side.
The increased importance of QL is a result of a much colder and dryer climate. The SMB is positive over most of Greenland during the LGM, in the absence of melt, even for extreme parameter combinations, so that the only change in mass balance is due to sublimation and hoar formation. The dry climate also favors increased sublimation. While during the PD the sign 10 of the vapor flux also undergoes an intra-annual variability, it is mainly a net heat loss for the surface during the LGM. The other sensitive parameters still determine the strength of the sublimation via the surface temperature, because the exponential impact of the surface temperature via Clausius-Clapeyron relation impacts the latent heat flux more than the actual exchange coefficient D SH (eq. 13). The large impact of D SH in the south-east is due to the large temperature difference between the surface and the air-temperature. This region is dominated by intense precipitation and rather warm air masses even during 15 the LGM. There is a precipitation gradient from the western coast to central Greenland, the fewest precipitation amounts are found in the west of Greenland and the south of Elsmere Island due to the presence of the Laurentine ice sheet. Therefore, the albedo module is more important on the western than the eastern margin. The reduced model sensitivity towards ✏ atm is mainly a result of the reduced atmospheric temperature with annual averages being around 10 K lower, resulting in less long-wave radiation and a lower absolute impact of the emissivity.

20
The sensitivity of the surface albedo, turbulent latent heat flux, snow melt and surface temperature will only be briefly discussed and without additional figures. The global sensitivity for the annual average albedo during PD period is mainly influenced by the fresh snow albedo, with only minor importance of the final firn albedo, ice albedo, the choice of the albedo module, and the atmospheric emissivity at the ice sheet margins. However, this is an interesting result for the model tuning as the atmospheric emissivity should be treated similarly to firn and ice albedo for tuning of the albedo routine.

25
Besides the switch which disables the turbulent latent heat flux Q L completely, it is most sensitive to the atmospheric emissivity ✏ atm followed by the turbulent heat flux exchange coefficient D SH , mainly around the margins. In these areas the ratio of turbulent sensible and latent exchange coefficients r lh/sh is a parameter of low sensitivity, as are the snow albedo related ones in the north. The turbulent latent heat flux Q L is not sensitive to the ice albedo and r lh/sh in the interior of the ice sheet or the maximum liquid water content ⇣ max globally. This shows that using similar exchange coefficients for moment and 30 water vapor are justified within the framework. During the LGM Q L shows an increased sensitivity to ↵ fs , while the D SH is less important around the margin, due to lower atmospheric temperatures and slower albedo decay. The albedo module and the firn albedo play almost no role in either case.
The average snow temperature is mainly influenced by ✏ atm , ↵ fs , D SH , QL , though uncertainties of the sensitivity are rather large, due to temperature reseeding in the event of snowfall on ice or shallow snowpacks.  . The analysis was done on a regional bases. During the present day there are three regions in the west and north from 0-1000-2000-3000 m. The southeast region is precipitation driven and the change in SMB with altitude is less developed, therefore 1000-3000 m is one region. There is one additional region in the center which is at elevations above 3000 m.
Snow melt is more less linked to the SMB, and shows now significant deviation from the sensitivities reported for the SMB.
Just as expected the impact of the latent heat flux switch is much lower.

Ensemble statistics
The GSA clearly highlights that the mass balance is most sensitive to the atmospheric emissivity and the latent heat flux switch irrespective of background climate. A drawback of this method is that it provides no information about the sign and absolute   Figure 6. The ensemble statistics for the surface mass balance for region 5 (west 1000-2000 m) shows the 5/95, 25/75, 33/66 and 50 quantiles in progressively darker shading. Black points represent the ensemble mean and the grey points correspond to the rest of the ensemble, apart from outliers (max 5 per bin allowed), which are removed to improve readability. Each plot represents the range of one parameter with ↵ fs , ↵ fi and ↵i in the top row, DSH, r lh/sh and ✏atm in the middle and Q Lon/of f , ↵ and ⇣max at the bottom. As Q Lon/of f and ↵ have two and four discrete values, the parameter range is not split in 20 intervals. is necessary as the ensemble was not created with parameters at regularly spaced intervals. It shows the behavior of the entire or a sub-ensemble with regard to one parameter. This is most interesting for the parameters where the effect changes sign depending on the atmospheric conditions. In Figure 6 the impact of the various parameters on the PD SMB is shown for region 5. This is the western region of Greenland ranging from 1000 -2000 m. The dominant parameter is the atmospheric emissivity ✏ atm . Over the range of Rephrase. This is worth expanding on. Maybe use phrases like "emissivity and surface mass balance are inversely correlated", "non-linear relationship between emissivity and surface mass balance", etc.
spread is decreasing. This will be discussed in detail later. Ice albedo ↵ i , liquid water content ⇣ max and r lh/sh have almost no impact. The SMB increases in the presence of turbulent latent heat flux due to the heat loss of sublimation in this region. With all albedo modules the ensemble has a wide spread, but the variation is smallest for the time dependent decay (Oerlemans and Knapp, 1998) which also has the highest median mass balance. The parametrization based on temperature, wetness and time (Aoki) has the lowest median SMB. 5 The strong impact of the emissivity is due to the larger annual average of the incoming long-wave radiation relative to the other fluxes. At very low atmospheric emissivity ✏ atm the largest energy source during the PD climate is reduced, and hardly any melt occurs in this region. In the absence of melt, mass balance response is only due to the sublimation and since this process has a modest absolute impact on the SMB the spread of the ensemble is low. Vice versa at large ✏ atm values a lot of grid cells will experience a heat up and potential melt of the snow-pack early leading to a positive feedback effect 10 with lower albedo. In agreement with the GSA (fig. 3) the firn albedo ↵ fi is almost as important as fresh snow albedo ↵ fs .
The ELA is located in this region and snow temperatures are quite warm, resulting in a large impact of snow albedo decay and its parametrization. The constant albedo parametrization has quite a wide spread, as it is very sensitive around the ELA, changing instantly from ↵ fs to ↵ fi when the snow-pack reaches the melting point. While the albedo parameterizations based on Oerlemans and Knapp (1998) and Bougamont et al. (2005) have a temporal decay relative to the instant one in the constant 15 case, with Bougamont having a snow temperature dependent increased decay rate and even more in the wet case. This leads to the low SMB tail of it relative to the other. The last albedo parametrization decays faster than the other two at warmer temperatures, leading to already lower albedo than all other parametrizations before the melting point is reached. Additionally, there is a difference between the models in case of snowfall. Oerl, const and Boug reset to ↵ fs and decay thereafter, while Aoki increases the albedo after snowfall incremental with the amount of snowfall.

20
During the LGM the western region between 1000 -2000 m (not corresponding to the identical geographical area as PD due to topographic changes) shows positive, but low SMB with a much lower spread of the ensemble. There are three distinct changes: The impact of ↵ fi and ↵ is drastically reduced, QL results in a decrease of SMB and the importance of ↵ fs increases relative to ✏ atm . The climate in the west of Greenland in the LGM is characterized by a much lower air temperature, slightly more annual mean radiation and lower precipitation. The lower air temperature reduces the impact of ✏ atm and 25 produces a more positive SEB, which in turn results in lower snow temperatures leading to slower snow albedo decay, no melting and almost no impact of associated parameters. The incremental increase of snow albedo with snowfall with the Aoki parametrization, gives slightly lower albedo in the dry climate. During the LGM sublimation prevails over the entire year, but in the absence of melt it results in a mass loss rather than a mass gain as in PD (via cooling and associated reduced melt).
The turbulent sensible heat flux may either be a heat loss or heat gain for the surface, making it particularly interesting. The I Maybe start with e.g. "In the PD climate, the largest energy source for melt is incoming short-wave radiation. As atmospheric emissivity decreases, less energy for melt is available. (Or similar) little surface ...warming, and potentially early melting of the snowpack, occurs in many grid cells... Be more precise, quantify.

Reference Figures
The strength of...
Rephrase "the figure is limited to" to "Only ... "is shown in Figure 7 because..."  show a distinct decrease in surface mass balance spread of the ensemble with increasing D SH and slight decrease of the mean (excluding region 6). The smallest spread of the ensemble is found in the high altitude-regions 9, 10, 11. The general slightly negative trend for the SMB with D SH is a result of higher air-temperatures than snow surface temperatures. This shows a very strong influence at the lower regions of the west and south, where warm air advection is frequent, resulting in increased melt due to the heat supplied by the sensible heat flux. The regions are also quite moist, in particular the south east (region 4), 5 therefore the latent heat flux during summer is also positive. The north-east (region 2, 3, 6, 7) of Greenland is colder and drier resulting in lower temperature and moisture differences between surface and atmosphere. This reduces the effect that D SH has on the SMB and SEB. The ensemble spread is narrower at higher altitudes (7-11) due to a generally positive mass balance, so an increased energy input due to any parameter will mainly impact the snow temperature, but in the absence of melt the mass balance does not change significantly. Still the change in temperature alters sublimation accounting for the remaining 10 variability. This is pronounced in region 10, at low exchange coefficients melt is still possible, if other parameters result in a strong positive energy input, leading to a skewness towards negative SMB. At a stronger exchange the tight coupling with the heat reservoir of the earth limits melt and the skewness towards negative SMBs vanishes, while the median stays almost constant. The parameter D SH is not the main driver for the SMB in this region, rather it acts as a buffer of the SMB. At strong turbulent sensible heat exchange the surface temperature will be buffered by the air-temperature heat resovoir. This effect is 15 visible for most of the regions, with air-temperatures most of the year below the melting point (2,3, 5, 7-10).
The earlier discussion mainly focused on the impact of D SH on the surface mass balance via the turbulent sensible heat flux Q SH , but it also impacts the turbulent latent heat flux (eq. 12). In a similar manner to Figure 7 the turbulent latent heat flux alone was analyzed, but is not shown here. Q L gets more negative with increasing exchange coefficient. There are two effects present, the increased surface water vapor pressure due to higher surface temperature as result of Q SH and the actual 20 water vapor exchange rate. The annual average of the ensemble over Greenland is a negative turbulent latent heat flux, meaning sublimation occurs more often than condensation. r lh/sh does increase Q L in its absolute value and therefore leads to more sublimation, but the effect on the mean is also lower than on the whole ensemble. At larger exchange rates more mass can be moved and therefore the variation over the ensemble increases as larger surface temperatures are prevailing.
The average SMB increases at low elevation if QL on, and decreases above 2000 m and is almost constant above 3000 25 m, where the annual average of the latent heat flux is almost zero. The first is a result of reduced melt due to a negative latent heat flux (sublimation), and it is therefore most pronounced in the north-east. The latter is a region where melt only occurs for a couple of extreme ensemble members, so the sublimation is the only mass loss and therefore the SMB decreases. The other parameters show analogous behavior to the GSA: An increase in ✏ atm decreases the SMB. The effect is not purely linear, rather increasingly more melt in case of a high emissivity. This is mainly related to the all year around impact of this parameter 30 altering ice exposure and albedo decay. The snow albedos increase the SMB, with the fresh snow albedo being more important.
The SMB ensemble plots do not show any dependency on ice albedo, liquid water content and the mean is also unaffected by Not a useful statement.
is negligible above 2500 instead of above 3000 m and melt tails are limited to below 2000 m. ✏ atm has less influence during the colder period, as well as the snow related parameters.
Additionally, the sensitivity of the 10 m firn temperature was assessed. This was limited to locations where actual firn cores have been taken . The firn temperature at these mostly central locations is sensitive to three parameters: ✏ atm , ↵ fs and D SH . The first two have a rather linear impact with increasing temperature with increasing emissivity and decreasing 5 albedo. The turbulent heat exchange coefficient has a similar effect on the 10 m firn temperature as on the SMB in region 5. At small D SH the 10 m temperature has a larger variability depending on the other two sensitive parameters, while at large values the air-temperature buffers the snow temperature. The sensitive parameters are the main drives of the SEB, which determines the 10 m temperature, in those areas as discussed previously for the SMB. The only difference to regions at higher elevation of the SMB is the insignificance of QL . The turbulent latent heat flux has only a minor importance for the SEB, but in the 10 absence of melt it is a requirement for SMB changes. These findings highlight, that by tuning the model for the 10 m firn temperature only information about the sensitive parameters could be extracted.

Discussion
The model sensitivity in this study is evaluated based on parameter uncertainties. Cloud cover and the associated atmospheric radiation are big uncertainties in the climate state. This is represented by the large parameter uncertainty of the atmospheric 15 emissivity ✏ atm . The SMB is most sensitive to this parameter under PD conditions. The sensitivity of the SMB is not drastically different during the LGM, but lower atmospheric temperatures reduce the impact of Q LW in and ✏ atm , while also leading to fewer areas where melt and runoff occurs. The relative contribution of the turbulent latent heat flux to the SMB increases drastically during the LGM (4 to 15 % of the total mass flux), making it the model's most sensitive parameter in large parts of the ice-sheet. This is due to the absence of melt, similar to the highest elevations during PD climate. Additionally, SMB values 20 have a smaller magnitude, as precipitation is less during the LGM and therefore the relative contribution of the vapor fluxes to the SMB is larger. For an accurate modeling of the SMB over the glacial cycle an inclusion of the turbulent latent heat flux is necessary, which may not be as important for a warmer climate.
The sensitivity reported is the result of two components: the absolute strength of a particular flux and the parameter uncertainty range of the parameter. The incoming long-wave radiation is the biggest single energy source for the surface during PD, 25 ranging from twice the incoming solar radiation around the margin to about 1/3 in the interior. The impact of the atmospheric emissivity is decreasing from the coast inwards for multiple reasons: Firstly, temperatures are higher around the margin leading to increased incoming thermal radiation. Secondly, cloud cover is larger and therefore solar radiation is reduced. Thirdly, at negative SMB the rate of melting depends greatly on the total energy flux and albedo decays faster at a warmer snowpack/ice if more long-wave radiation reaches the surface during winter. It is important to differentiate between the sensitivity of the del Last E.g.: "In this study we asses model sensitivities due to parametric uncertainties." largest The Greenland-wide mass balance is most sensitive to the atmospheric emissivity during PD ( fig. 2), while during the LGM the SMB shows increased sensitivity to the fresh snow albedo, the choice of albedo parametrization, and the turbulent latent heat flux. This is a result of lower atmospheric temperatures and an ice sheet with less melt increasing the importance of the vapor fluxes similar at high elevations at PD. While it would be acceptable for the integrated SMB to not include Q L in a warm climate using this model, it is not during a colder and dryer climate state. Additionally, this means that although cloud cover 5 uncertainty may be similar during the colder period of the LGM the sensitivity of this model towards this uncertainty is lower.
At the local scale even in PD climate the vapor flux (Q L ) is up to 1/3 of the total SMB, despite its small impact on the Greenland wide integrated scale. Additionally, in the absence of melt the sublimation and condensation are the only changes to the SMB with a fixed precipitation forcing. Neglecting these processes will result in fundamental biases over long-term simulations of the Greenland ice sheet, in particular in its interior. This needs to be considered when tuning surface mass 10 balance models for long time scales. A tuning for the Greenland wide SMB will mainly constrain the most sensitive parameters, which constrain certain key regions of high mass turnover, but not for the bulk surface. Furthermore, the temporal change of the sensitivities on the local scale indicates that models of reduced complexity may fail drastically for other time periods (absence of Q L for example).
The sensitivity in this study is not defined relative to the absolute flux, it does not tell us that it is most sensitive to the long- 15 wave radiation, but rather that the uncertainty of the long-wave radiation has a larger impact on the SMB than the uncertainty in the incoming solar radiation on the SMB uncertainty. The impact each energy flux has on the absolute SMB has to be analysed separately. The long-wave radiation dominates, followed by the solar radiation. The increased sensitivity of the QL switch during the LGM is mainly due to its increase on the absolute SMB. In the absence of melt, and with reduced precipitation sublimation accounts for a larger portion of the absolute SMB.

20
It is beyond the scope of this model to resolve all the physical processes. We use a simple parametrization for the incoming long-wave radiation which does not accurately represent reality. The atmospheric emissivity is neither constant in space nor time. Area distributed values may work during the observational period, but changes in the atmospheric circulation may change these patterns over the glacial cycle. We conclude that the overall model uncertainty can effectively be improved by changing the simplified representation of the long-wave radiation flux as a function of atmospheric temperature and emissivity for either 25 a more sophisticated parameterization or to use long-wave radiation as climate model input. If the model is to be tuned for the Greenland-wide SMB it will be biased towards the melt regions around the margin and therefore the atmospheric emissivity.
Where possible parameters should be calibrated via quantities which they are sensitive for. Neither of our parameters are to be assumed constant in space, as albedo for example depends very much on impurities and snow temperature, but unless the uncertainty of the incoming long-wave radiation is reduced, it is justifiable to work with optimized values from PD. The found 30 model sensitivity towards the parameters is to be set into context with the assumed forcing (climate) uncertainty.

I
Rephrase to make clear that you mean models that calculate SMB don't need to include Q_L strongly depends on

Conclusions
The surface mass and energy balance model BESSI has been improved by accounting for turbulent latent heat flux and snow aging. The sensitivity of the model towards these changes and other uncertain parameters was assessed with a variance based sensitivity method based on two ensembles with a total of 16500 simulations. Warm present day climate and the cold last glacial maximum were used to study the change of the model response under different boundary conditions. The sensitivity 5 analysis reveals that the inclusion of the turbulent latent heat flux is a necessity to simulate the local SMB and the integrated SMB over the entire Greenland ice sheet. The relative importance of sublimation and condensation is larger in the dry and cold climate of the LGM, as air temperature and precipitation is lower.
The uncertainty associated with cloud cover and atmospheric emissivity creates a SMB model uncertainty. With the change in circulation during the last glacial a changing energy input from the atmosphere to the surface will result in a SMB response.

10
Simple surface mass balance models like enhanced temperature index models are likely to create a bias as solar insolation changes are accounted for while the impact of the cloud cover on other components of the energy balance are neglected or implicitly included in the PDD factor, this is even more true as even under PD conditions the impact of clouds on the SMB is highly variable (Van Tricht et al., 2016;Hofer et al., 2017). The sensitivity study further reveals that the uncertainty of the SMB as a result of the atmospheric radiation decreases in colder climate. 15 The current model version does not use wind fields, although this impacts the SEB via the turbulent fluxes. The strength of the turbulent exchange does not have a big impact on the SMB and the approach to neglect wind speed variability is therefore justified in the context of more uncertain parameters.
We find that uncertainties in the ice albedo, liquid water content and differences of the turbulent fluxes are of minor importance for our and likely also similar models. In order to reduce model uncertainty most effectively, first the larger energy 20 sources of short wave and long wave radiation need to be constrained via the snow albedo and the atmospheric emissivity.
Additionally, the github branch also contains the analysis and plotting scripts. Upon acception the branch will be archieved via zenodo.
Author contributions. TZ implemented the model changes, conducted the ensemble simulations, the sensitivity studies and the data analysis and wrote the main part of the manuscript. AB contributed to all aspects of this study.

25
Competing interests. The authors declare that they have no conflict of interest. Here you introduce a new idea/model that has not yet been discussed.