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⁻³ (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⁻³ (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⁻³ 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 ₀ is valid for ρ < 550 kg m⁻³ and

where k ₁ is valid for 550 kg m⁻³ < ρ < 840 kg m⁻³. k ₀ and k ₁ 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 ₀ and k ₁, are given as:

where T is the near-surface air temperature in Kelvin and R = 8.314 J K⁻¹ mol⁻¹ 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 k_z 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. ρ ₀ and ρ _e = 960 kg m⁻³ 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, ρ ₀. 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 ρ ₀ = 300 kg m⁻³ and c/L = 0.006°C⁻¹, 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. 1. PR = 0,

2. 2. ,

3. 3. ,

4. 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 ‘p_{p 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.

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 km³ a⁻¹ 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 (km³ a⁻¹ 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⁻¹ 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.