Introduction
Global climate models predict a large warming at high latitudes in the coming centuries (Reference Solomon and SolomonSolomon and others, 2007). The response of the Greenland ice sheet (GrIS) to climate change is a major factor in understanding sea-level changes, and estimates of the current mass balance of the ice sheet have many uncertainties (Reference BoxBox and others, 2006; Reference BougamontBougamont and others, 2007; Reference FettweisFettweis, 2007; Reference HannaHanna and others, 2008). An important component in the investigation of an ice sheet is the annual surface mass balance (SMB) since it has a direct relation to atmospheric variability. Several ice-sheet models have estimated the annual SMB of Greenland without inclusion of the potentially important densification of snow, which enables a separation of volume and mass changes (Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999; Reference Tarasov and PeltierTarasov and Peltier, 1999; Reference GreveGreve, 2005). Reference ColbeckColbeck (1976) did pioneering work on meltwater flow in snowpacks and defined a runoff limit. He stated that before any runoff can occur the meltwater must refreeze within the snow cover and raise the temperature to overcome the cold content. In addition, the pore storage of the snowpack must be water-filled to overcome capillary forces. Based on the work of Reference ColbeckColbeck (1976), Reference Pfeffer, Meier and IllangasekarePfeffer and others (1991) argued for and presented the most important parameters for describing the cold content and pore storage and performed a modelling study on retention of meltwater from the GrIS. Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) suggested a parameterization based on the work of Reference Pfeffer, Meier and IllangasekarePfeffer and others (1991) in order to have retention processes included in large-scale ice-sheet models. The latter two studies did not include densification processes which have a direct effect on the volume of pore storage in a snow layer (Reference ColbeckColbeck, 1976). Neglecting densification of the annual snow layer, especially in the lower percolation zone, could have a large impact on mass changes when calculating the total amount of melt and the subsequent runoff (Reference Zwally and LiZwally and Li, 2002; Reference Reeh, Fisher, Koerner and ClausenReeh and others, 2005; Reference ReehReeh, 2008).
The motivation for including densification processes in the SMB of large-scale ice-sheet models is to estimate the meltwater retention more accurately. Without densification the calculated meltwater retention of Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) could be overestimated in the SMB since the cold content and pore storage depend on the annual near-surface air temperature, and the thickness and density of the annual snowpack (Reference Bøggild, Forsberg and ReehBøggild and others, 2005).
Our study investigates the impact of the inclusion of different retention schemes in the SMB module of a large-scale GrIS model. The retention schemes will deal with a coupling between snow densification and retention processes of the annual snowpack. Furthermore, we investigate the inclusion of a temperature standard-deviation parameterization in order to constrain surface melt calculations modelled from a positive degree-day (PDD) method. The aim is to improve SMB estimates for past, present and future scenarios of the GrIS. The presented SMB parameterization is tested for present-day conditions through comparison to observations.
Methods
Firn-densification model
Two stages of dry snow densification are generally considered. The first stage is densification due to structural rearrangement of grains by grain-boundary sliding. This process is dominant up to a critical density of ∼550 kg m−3 (Reference Herron and LangwayHerron and Langway, 1980). The second stage comprises sintering, which involves transfer of material by sublimation, and molecular diffusion within grains to redistribute mass (Reference PatersonPaterson, 1994; Reference HookeHooke, 2005). At the end of this stage, air bubbles are trapped with no contact with the atmosphere. This is the close-off depth and it occurs at a density of ∼840 kg m−3 (Reference Herron and LangwayHerron and Langway, 1980; Reference PatersonPaterson, 1994).
Reference Herron and LangwayHerron and Langway (1980) developed a simple model, hereafter referred to as the HL model, that describes the densification process when no meltwater is present. The model is forced by the mean annual near-surface air temperature and the mean annual accumulation, which can be obtained from automatic weather stations, snow pits, etc. (Reference Reeh, Fisher, Koerner and ClausenReeh and others, 2005). The basic assumptions are that the mean annual accumulation rate, A, and the near-surface air temperature, T, at a given location are constant with time and that the compression of firn is linearly related to the change in stress from the weight of the overlying snow. This can be expressed as:
where ρ is the density, ρ i = 917 kg m−3 is the density of pure ice and C is a constant. In order to find the densification rate, dρ/dt, we assume that dh/dt = Aρ i /ρ, where A is the accumulation rate ice equivalent. Following the notation of Reference Herron and LangwayHerron and Langway (1980), it is also assumed that the temperature and accumulation rate can be separated in the densification rate formulation:
where k 0 is valid for ρ < 550 kg m−3 and
where k 1 is valid for 550 kg m−3 < ρ < 840 kg m−3. k 0 and k 1 are Arrhenius-type rate constants that only depend on temperature, and a and b are constants. Their values are empirically determined to be a = 1.1 ± 0.2 and b = 0.5 ± 0.2 (Reference Herron and LangwayHerron and Langway, 1980). The Arrhenius-type rate constants, k 0 and k 1, are given as:
where T is the near-surface air temperature in Kelvin and R = 8.314 J K−1 mol−1 is the gas constant. The HL model was developed in order to aid the interpretation of ice cores, and in the current form it is not suitable for modelling seasonal variation, due to the assumptions of constant accumulation rate and near-surface air temperature.
The model of Reference Zwally and LiZwally and Li (2002) (ZL model) is based on the HL model, but uses a different approximation for the Arrhenius-type rate constants:
where E = 883.8|T|− 0.885, K 0G = 8.36|T|− 2.061 and β = 8.0 is an empirical constant for adjusting for differences in grain growth during densification. The densification rate becomes:
where kz is valid for all densities and a = 1.0. The ZL model is suitable for modelling the seasonal variation of snow densification induced by a steady-state temperature cycle (Reference Li and ZwallyLi and Zwally, 2002). We use the ZL model in this study, forced by the parameterized annual near-surface air-temperature cycle and the mean annual accumulation rate.
Integrating Equation (7) over 1 year shows the need for a surface density. Reference Reeh, Fisher, Koerner and ClausenReeh and others (2005) derived a relation for surface density which depends on the mean annual near-surface air temperature, T ma, in °C:
The parameterization is based on density observations derived from shallow ice cores and snow-pit measurements on the GrIS.
Degree-day model
To estimate the amount of meltwater that percolates through the snowpack we use a positive degree-day (PDD) approach. This is a statistical approach that relates positive air temperatures to the melting of snow or ice (Reference BraithwaiteBraithwaite, 1985, Reference Braithwaite1995; Reference ReehReeh, 1991). The PDD factors linearly relate the sum of positive near-surface air temperatures to melt, and the percentage of days with melt is assumed to be equal to the probability that the near-surface air temperature exceeds 0°C (Reference ReehReeh, 1991; Reference BraithwaiteBraithwaite, 1995; Reference HockHock, 2003).
The approach in this study largely resembles the method of Reference GreveGreve (2005) and is briefly summarized here. The percentage of PDDs from the normal probability distribution around the monthly mean temperatures during 1 year is produced in the following way:
where t is the time, T (°C) is the actual near-surface air temperature, σ pdd is the standard deviation of T and T a (°C) is the annual mean near-surface air-temperature cycle represented as:
where t a is 1 year. T mj and T ma are the mean July and mean annual near-surface air temperatures; they are taken from the temperature parameterization developed by Reference Fausto, Ahlstrøm, Van As, Bøggild and JohnsenFausto and others (2009), which depends on altitude, z s, latitude, φ, and longitude, λ.
The annual melt, M, in ice equivalent is then calculated using the PDD integral and PDD factors β snow,ice in the following way:
where β
ice, β
snow are assumed to be different for ice- and snowmelt, and for warm, 
, 
, and cold, 
, 
, climate conditions. The PDD factors for cold and warm climate conditions are assumed to differentiate based on the mean July near-surface air temperature, T
mj, in the following way:
where 
, 
and 
. T
w = 10°C and T
c = −1°C are the limiting temperature values, which are used to calculate the PDD factors, limited by warm or cold conditions (Reference GreveGreve, 2005). The σ
pdd in the PDD integral (Equation (9)) is different to the uniform value of 4.5°C used by Reference GreveGreve (2005).
The σ pdd value is based on a distributed parameterization constructed in a similar way to the temperature parameterization of Reference Fausto, Ahlstrøm, Van As, Bøggild and JohnsenFausto and others (2009) using the same dataset. Monthly values are calculated from hourly temperature observations each month in a given year for the Greenland automatic weather stations (AWSs) (Fig. 1; see Reference Fausto, Ahlstrøm, Van As, Bøggild and JohnsenFausto and others, 2009, for details). Subsequently, the mean summer standard deviations are calculated for each station using all available monthly values from June, July and August (σ jja) (Table 1).
Fig. 1. Locations of automatic weather stations on the GrIS.
Table 1. A comparison between the modelled (σ m) standard-deviation distribution and observed data (σ jja) from the stations. σ jja is the mean summer (June, July, August) standard deviation. The difference (Diff.) is calculated between the modelled and observed data. Acc. and Abl. denote stations located in the accumulation and ablation zones, respectively
A least-squares fit was made to the observed σ jja values, assuming a linear dependence on altitude, z s, latitude, φ, and longitude, λ:
where σ m is the mean summer (June, July, August) standard-deviation parameterization (Table 1). The values of the coefficients for Equation (14) are given in Table 2.
Table 2. Summertime coefficients for Equation (14) and their root-mean-square difference (RMSD) relative to the observed standard deviation
Retention models
In order to calculate how much meltwater is retained in the snowpack, a retention model is needed. Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) describe and compare four different meltwater-retention models based on the SMB model of Reference Huybrechts and de WoldeHuybrechts and de Wolde (1999). They calculate the sum of PDDs using the mean daily temperature; then melt, percolation, refreezing and retention values are calculated. The runoff models add rain and meltwater together in the calculations, so they implicitly take into account the warming of rain in the snowpack. Their most comprehensive retention model is taken from Reference Pfeffer, Meier and IllangasekarePfeffer and others (1991) and accounts for refreezing and pore storage of the snowpack (Reference Janssens and HuybrechtsJanssens and Huybrechts, 2000).
Following Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) the retention of meltwater in the snowpack is defined in terms of a potential- retention fraction, PR, meaning that a certain fraction of the snowmelt and rainfall in relation to the annual snowfall needs to be retained by refreezing or capillary processes before runoff can occur. PR is computed in the following way:
where c is the temperature-dependent specific heat capacity and L is the latent heat of fusion. A, M and P are the snowfall, snowmelt and precipitation, respectively. ρ 0 and ρ e = 960 kg m−3 are densities of the surface snow and water-saturated wet snow, and d is the thickness of the thermally active layer that provides the cold content. Here the thickness of the thermally active layer providing the cold content is limited to the thickness of the annual snowfall with a constant density, ρ 0. It is assumed that retention and refreezing processes will only happen within the annual snow layer and remain passive thereafter. Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) use ρ 0 = 300 kg m−3 and c/L = 0.006°C−1, because c is nearly constant for these temperature variations.
Equation (15) forms the basis of a comparison of different retention schemes choosing different combinations of input parameters. The following retention schemes (‘model setups’) will be compared:
-
1. PR = 0,
-
2.
, -
3.
, -
4.
.
,
,
.A, M and P are calculated in metres ice equivalent. T is assumed equal to the absolute value of the annual mean near-surface air temperature, T = |T ma|, and set to 0°C for positive temperatures. The values of d and the density in all model set-ups depend on how the surface density is defined or calculated. The surface density is defined here as the mean density of the annual accumulated snowfall.
Model set-ups 1, 2 and 3 are identical to ‘pp 0.0’, ‘p pmax’ and ‘p capil’ introduced by Reference Janssens and HuybrechtsJanssens and Huybrechts (2000). In model set-up 4, Equation (8) is chosen for the surface density and is then integrated in Equation (7) to achieve a mean annual snow layer density (ρ smean) for the retention model.
Results
Temperature standard deviation
PDDs show a high sensitivity to changes in the value of the standard deviation, σ pdd. Reference Fausto, Ahlstrøm, Van As, Bøggild and JohnsenFausto and others (2009) demonstrate that a uniform increase of σ pdd from 2.53°C to 4.5°C results in a 33% increase in the modelled melt area over Greenland, so it is important to constrain the σ pdd value with observations. The PDD model uses a new parameterization, which is based on observations, for the standard deviation, σ pdd, of the near-surface air temperature (Equation (14)). This parameterization takes into account and reflects more closely the spatial variability of the near-surface air temperature over the GrIS.
The altitudinal component dominates the standard-deviation field, and the effect of the latitudinal and longitudinal components is not very large, as indicated by the three fits listed in Table 2. The fit that depends on altitude, latitude and longitude is chosen as the parameterization for the standard deviation, because it has the smallest rootmean-square difference (RMSD), which indicates the best fit. Figure 2 presents the standard-deviation field computed from the parameterization (Equation (14)). The highest standard deviations are found in the parts with high surface elevation (σ pdd ∼ 5°C), and smallest standard deviations are found at low elevation (σ pdd < 2°C). The standard deviation has a distinct annual cycle, with the largest values during winter and the smallest during summer (not shown). The smaller values during summer can be explained by melt of snow and ice. The surface temperature of a melting snow or ice surface is fixed at 0°C, limiting the variability of the near-surface air temperature. Except for summer, no melt occurs and the surface temperature variations are not limited by the melting-point temperature, allowing a larger variability of the near-surface air temperature. May, June and September account for the highest standard-deviation values of the ablation season (3.0–6.0°C), and the lowest values (<2.0°C) are found in July and August. The temperature variation arises mainly from the diurnal cycle and from synoptic-scale disturbances (Reference Lefebre, Gallée, van Ypersele and HuybrechtsLefebre and others, 2002; Reference Fausto, Ahlstrøm, Van As, Bøggild and JohnsenFausto and others, 2009).
Fig. 2. Parameterized standard deviation in °C for the summer (June, July, August).
Figure 3 shows the observed standard deviation from the AWSs and the difference between the standard-deviation parameterization and the observed values relative to elevation (also listed in Table 1). The standard deviation has a strong dependence on altitude. The stations at lower altitudes show a larger scatter in their calculated values than the stations located at higher elevations. This is due to the higher exposure of low-elevation AWSs to atmospheric variability over land and ocean. This is illustrated at Sermilik1, the southernmost and most exposed of the AWSs on the GrIS. Here, the observed standard deviation is relatively large, causing the largest mismatch in our fitting procedure. Overall, the parameterization approximates the standard deviation within 1°C, and has increasing accuracy with elevation, which benefits the meltwater retention scheme, as this mainly concerns the lower part of the accumulation zone.
Fig. 3. (a) Standard deviation calculated from hourly temperature observations relative to elevation. (b) Difference between the standard-deviation parameterization and the observed values from the AWS relative to elevation.
Potential retention fraction
Table 3 lists the values of the corresponding SMB components for the retention schemes in relation to the total model precipitation that amounts to 575 km3 a−1 ice equivalent. Melt, runoff and retention are calculated in order to quantify the changes when densification is accounted for in the SMB modelling. Model set-up 1 is without retention and all melt is lost by runoff; it has the highest value of all the set-ups. It is considered as the lower bound in relation to the importance of the different components in the other retention schemes. Model set-up 2 calculates the largest retention value and the lowest runoff, which is due to the constant and uniform value of PR. The distribution of runoff is limited to the ablation zone, and the distribution of the retention is limited to the accumulation zone; 40% of the calculated meltwater is retained. Model set-up 3 takes into account the cold content and the pore storage of the snowpack. The distribution of runoff now comes from both the accumulation and ablation zones, and 29% of the calculated meltwater is retained (Table 3). Model set-up 4 includes densification of the annual snow layer. This has a significant effect on the thickness of the thermally active layer, d, and on the pore storage through the calculation of the mean annual snow layer density, ρ smean, in the retention model. The distributed amount of runoff is here slightly higher from the accumulation zone than in model set-up 3, and 23% of the total meltwater production is retained within the snowpack (Table 3). The influence of densification on retention in the annual snowpack is clearly due to less efficient pore storage. Moreover, the thermally active layer for the cold content gets thinner.
Table 3. Calculated present-day SMB (km3 a−1 ice equivalent) for the GrIS using four model set-ups. Abl. and Acc. denote the ablation and accumulation zones, respectively
Table 4 lists the observed and modelled SMB values for the location of seven AWSs. To test the performance of the new retention scheme, the SMB is calculated using the approach of Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) and model set-up 4 for observed standard deviation values (σ pdd = σ jja), which is then compared to observations. Comparing model set-ups 3 and 4 with observation indicates a better performance for model set-up 4 based on the RMSD (Table 4). Model set-up 4 performs better, especially at Swiss Camp, JAR1 and Petermann ELA, because the stations are near the equilibrium-line altitude (ELA), which is close to the lower accumulation zone where the largest amount of retention is found. At the other station locations, the results are the same; this is because both models calculate the same value for PR (PR = 0). Model set-up 4 is also tested using a different standard deviation to that of Reference GreveGreve (2005) (σ pdd = 4.5°C). Here the RMSD indicate a better performance for the distributed standard-deviation parameterization (σ m).
Table 4. Comparison between the modelled and observed SMB (m a−1 ice equivalent) at the location of the AWSs in the ablation zone for model set-ups 3 and 4. The modelled SMB uses σ pdd values derived from observations (σ pdd = σ jja, Table 1) and Equation (14) (σ pdd = σ m, Table 1), and the value of Reference GreveGreve (2005) (σ pdd = 4.5). RMSD is the root-mean-square difference relative to the observed values
The total retention calculated from model set-up 4 is plotted in Figure 4a. The highest retention values are calculated in the lower accumulation zone near the ELA. The lower accumulation zone is clearly recognizable as a band of high retention values all around the perimeter of Greenland. The southwestern and northeastern parts of the ice sheet have a wide band with high retention. This could be due to a relatively small change in elevation in these areas, a small annual accumulation rate or a combination of the two. The difference in calculated total retention from model setups 3 and 4 is plotted in Figure 4b. This plot shows the direct effect of including densification processes in the retention schemes. Positive differences mean that the retention in model set-up 3 is greater than in model set-up 4. Here the highest differences are seen along the southeastern margin, which could be due to the relatively high accumulation rate in that area. Large differences are also seen along the southern and southwestern margins, where the highest melt rates are measured. Common to all the positive differences is, however, that they are located in the lower accumulation zone near the ELA.
Fig. 4. Total retention in metres ice equivalent for (a) model set-up 4 and (b) the difference between model set-ups 3 and 4. Brown denotes land.
The PR of liquid water on the ice sheet calculated from model set-up 4 is plotted in Figure 5a. Differences between model set-ups 3 and 4 are evident around the runoff limit in the lower percolation zone and can be seen as sharp changes in PR (Fig. 5b). The effects of densification result in a higher amount of runoff from both the ablation and accumulation zones, because densification has a negative effect on the ability of the annual snowpack to retain water. The potential retention and difference as a function of elevation is also plotted in Figure 6 as a scatter plot. Figure 6a clearly illustrates the large regional differences in Greenland, due to the relatively broad span of equal elevation values where the model calculates both PR = 1 and PR values close to 0. From the difference plot (Fig. 6b), it is evident that the largest differences in PR are found close to the 1500 m elevation line, which is the altitude for the lower accumulation zone for most places in Greenland. Similarly to the difference plot of the total retention (Fig. 4b), positive values here also mean a reduction in the snow layers’ ability to retain water. The densification processes have the largest influence on retention around the 1500 m elevation line.
Fig. 5. Potential retention fraction for (a) model set-up 4 and (b) the difference between model set-ups 3 and 4. Brown denotes land.
Fig. 6. Potential retention fraction as a function of elevation for (a) model set-up 4 and (b) the difference between model set-ups 3 and 4.
Retention processes reduce runoff from snow and should be included in SMB modelling (Reference Pfeffer, Meier and IllangasekarePfeffer and others, 1991; Reference Bøggild, Forsberg and ReehBøggild and others, 2005). Previous studies have dealt with meltwater retention calculated on the basis of the available energy and pore storage, which depend on the snow density of the surface layer. The inclusion of snow-density calculations in the surface layer has a significant effect on surface-meltwater retention, and the influence on the total runoff is a larger contribution to sea level. The densification has a negative effect on the ability of the annual snowpack to retain water, because it reduces both the thickness of the thermally active layer, assumed to contain the cold content, and the volume of pore storage for meltwater within the annual snow layer (Reference Pfeffer, Meier and IllangasekarePfeffer and others, 1991; Reference Reeh, Fisher, Koerner and ClausenReeh and others, 2005). Only dry snow densification is considered, since the retention model calculates the potential retention based on the cold content and pore storage before the onset of melt. The simple model based on Reference Herron and LangwayHerron and Langway (1980), that describes and calculates the densification, has been tested and used for many model studies of the GrIS (Reference Li and ZwallyLi and Zwally, 2002; Reference Zwally and LiZwally and Li, 2002; Reference Li, Zwally, Cornejo and YiLi and others, 2003, Reference Li, Zwally and Comiso2007).
Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) use their SMB model, using model set-up 3, to calculate that ∼80% of meltwater on the ice-sheet domain produces runoff, with the remainder refreezing in the snowpack and on the ice surface as superimposed ice. Here model set-up 3, which uses the retention scheme of Reference Janssens and HuybrechtsJanssens and Huybrechts (2000), calculates that 71% of the total melt from the PDD model produces runoff (Table 3). The reason for this discrepancy is that we use slightly different PDD factors and different temperature and standard-deviation parameterizations.
Discussion and Conclusions
Coupling snow-densification and meltwater-retention processes in the SMB for a large-scale ice-sheet model of Greenland achieves a separation of volume and mass changes of the surface layer, in order to determine the surface melt contribution to runoff. The separation of volume and mass changes is also important for the interpretation of repeat surface altimetry studies because the short-term (less than a few years) variability of snow densification in the surface layers will change the surface elevation (Reference Zwally and LiZwally and Li, 2002; Reference Li, Zwally, Cornejo and YiLi and others, 2003, Reference Li, Zwally and Comiso2007).
Superimposed ice formation and meltwater refreezing in the snowpack are key processes which should be incorporated into mass-balance models applied to the GrIS. Neglecting these processes results in a larger runoff estimate. In this study, the SMB of GrIS was calculated using the PDD approach, taking into account the processes of meltwater retention in the annual snowpack by pore storage and the meltwater refreezing (cold content) and snow densification. It was shown that including densification in the SMB results in a less efficient retention of meltwater, especially at the 1500 m elevation line in the lower accumulation zone.
The SMB has only been measured directly at a few selected sites. Since available records are short and characterized by interannual variability, a proper validation of the SMB calculation is difficult. However, in an attempt to address this issue, we test the performance of the new SMB by comparing model results with observations. The SMB is calculated for the approach of Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) (model setup 3) and model set-up 4. Moreover, model set-up 4 is also tested using a different standard deviation to that used by Reference GreveGreve (2005). The SMB calculations are listed in Table 4. The results in Table 4 indicate that the SMB calculation of this study perform better (RMSD = 0.61) than those of Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) (RMSD = 0.68). The improvements in the retention values are mostly located in the lower accumulation zone, and the highest values amount to ∼0.6 m ice equivalent along the southeastern margin (Fig. 4b).
The PDD model uses a new parameterization for the standard deviation (σ pdd) of the near-surface air temperature which takes into account the spatial variability over the GrIS. This parameterization is deduced from observations, so it closely reflects the variability of the near-surface air temperature. In this study, the new standard-deviation parameterization may not in full correspond to a climatological mean, due to the limited dataset available. The data, obtained over different periods and sometimes for only a few months, could introduce a bias in the modelled standard-deviation distribution. It is considered that the standard-deviation parameterization is a fair approximation for the period with good AWS coverage (1996–2006) (Reference Fausto, Ahlstrøm, Van As, Bøggild and JohnsenFausto and others, 2009). The parameterization is difficult to validate further in regions where no in situ data is available. The new standard-deviation parameterization is also tested, and the performance in the SMB calculation indicates an improvement compared to the constant value used by Reference GreveGreve (2005). Table 4 compares the observed SMB with that modelled using different values of σ pdd. The comparison indicates a better model performance using the new parameterization (RMSD = 0.61) for the standard deviation compared to the study by Reference GreveGreve (2005) (RMSD = 1.32).
Model parameterization for the standard deviation (σ pdd) has, to our knowledge, never been applied to present-day SMB studies of ice-sheet models of Greenland (Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999; Reference Tarasov and PeltierTarasov and Peltier, 1999; Reference GreveGreve, 2005). It is commonly assumed that σ pdd is approximated with one single value in the interval 4.5–5.5°C (Reference ReehReeh, 1991; Reference Ritz, Fabre and LetréguillyRitz and others, 1996; Reference Tarasov and PeltierTarasov and Peltier, 1999). Reference Fausto, Ahlstrøm, Van As, Bøggild and JohnsenFausto and others (2009) discussed and demonstrated that the PDDs show a high sensitivity to changes of σ pdd and that it is important to constrain σ pdd with observations. The standard-deviation parameterization for present-day summer climate in Greenland quantifies the distributed temperature variation, and the PDDs will reflect the variability of the climate more closely.
Using the mean summer distribution of the standard deviation for year-round ablation calculation may not be the best choice for σ pdd. Using an annual distributed cycle for σ pdd determined for the whole ice sheet would be a better choice, but there are insufficient data available. Another problem arises when modelling over long timescales, as we do not know the standard deviation of the mean temperatures during, for example, the last ice age. Therefore, it is emphasized that the σ pdd distribution (Equation (14)) may not be valid for past and future climates.
The simple assumptions in the PR formulation, which are needed to limit the computation time, have some limitations: Choosing the thermally active layer, d, to be only the annual accumulation layer is necessary to keep computation time low, but will have an influence on the development of runoff because it confines the energy exchange and filling of the pore space to that layer. The present scaling of the thermally active layer may be a source of error in the runoff calculation because field evidence above the equilibrium line shows that surface meltwater penetration is not confined to the annual accumulation layer (Reference Pfeffer, Meier and IllangasekarePfeffer and others, 1991; Reference Reeh, Fisher, Koerner and ClausenReeh and others, 2005). Also, due to the one-step method to calculate the annual accumulation layer, the annual near-surface air temperature is used for the cold content. This will impact the total retention during the melt season because the cold content will be variable (Reference Bøggild, Forsberg and ReehBøggild and others, 2005). Heat loss during the melt season, either by conduction into cold underlying snow or to a cold atmosphere during transient cold spells in the summer, exemplifies seasonal changes that cannot be resolved using the time-step of an ice-sheet model. The runoff calculated through PR in the surface-boundary conditions does not take into account how the water is channelled away from the ice sheet. The actual routing of mobile meltwater is determined, to a certain degree, by surface topography and subglacial topography when it escapes to the base of the ice sheet, either via moulins or crevasses. Mobility of meltwater, however, does not guarantee that the meltwater escapes from the ice sheet (Reference Wankiewicz, Colbeck and RayWankiewicz, 1979; Reference PatersonPaterson, 1994). Determination of the potential errors introduced by these simple parameterizations of the melt, retention and runoff calculations for the surface conditions requires extensive observational data input and more detailed models that resolve the entire melt season. However, comparing the model results presented in this study with observations yields an overall performance of the SMB calculation that is satisfactory for large-scale ice-sheet models of Greenland.
Acknowledgements
Thanks to K. Steffen’s group at CIRES and to M. van den Broeke at the Institute for Marine and Atmospheric Research Utrecht for providing temperature data from the GC-Network and the K-Transect. Many thanks to J. Li, T. Pfeffer and the Scientific Editor, W. Wang, for constructive criticism which improved the manuscript significantly. This paper is published with the permission of the Geological Survey of Denmark and Greenland.
