2.1. Model description

We used the one-dimensional multilayer snowpack model Crocus (Reference VionnetVionnet and others, 2012) that is currently implemented in SURFEX (‘Surface Externalisée’ in French), a surface modelling platform developed by Météo-France in cooperation with the scientific community (Reference MassonMasson and others, 2013). Crocus computes the energy and mass balance of the snowpack and accounts for surface melting, internal melting/refreezing, compaction, snow metamorphism, snow aging, near-surface densification and sublimation enhancement during snowdrift. As the model domain was limited to the thickness of the firn layer, bottom melting was not represented in the simulations. Crocus discretizes the firn stratigraphy with up to 50 vertical layers ranging in thickness from 0.05 m close to the surface, to 10 m at depth. Each layer is approximately characterized in terms of its density, liquid water content, temperature, thickness and snow grain parameters (age, dendricity, sphericity and grain size). The initial age of each of the layers, which is used in computing the surface albedo, was determined based on the rate of accumulation of ~ 0.23 ± 0.07 m w.e. a^–1 over the surveyed area of Devon Ice Cap (Reference Colgan and SharpColgan and Sharp, 2008). Given that rate, a thickness of 10 m w.e. represents 43.5 years of net accumulation. The parameters dendricity, sphericity and grain size are used in the computation of surface albedo, transmission of solar radiation, and the impact of wind drift on compaction and metamorphism of the surface layers.

In the absence of direct observations of snow grain parameters, we tested the model using several values for dendricity, sphericity and grain size. Although the initial values may affect how rapidly the model converges, they change rapidly and converge quickly on the same computed values. As well, new precipitation creates new layers close to the surface, which progressively reduces the impact of previous surface layers with arbitrary characteristics. As a result, we found that varying the initial values for snow grain parameters did not change the outcomes of the simulations. As the snow grain parameters are updated at each time step through snow metamorphism (Reference VionnetVionnet and others, 2012), we used the model’s default parameter values to initialize the physical properties of snow layers. The number and thickness of layers vary in time during simulations to represent the evolution of snowpack characteristics resulting from all the processes accounted for in the model.

Crocus computes the density of freshly fallen snow as a function of wind speed, U, and air temperature, T _a, using

(1)

where T _fus is the temperature of the freezing point for water, a_ρ = 109 kg m^–3, b_ρ = 6 kg m^–3 K^–1 and c_ρ = 26 kg m^–7/2 s^–1/2. The minimum fresh snow density was set to 50 kg m^–3. If the snowdrift option is activated, the last term on the right-hand side of Eqn (1) is removed.

The densification of snow is expressed as

(2)

where D is the layer thickness, σ the vertical stress, dt the model time step, and η the snow viscosity (which is a function of the snow temperature, density, liquid water content and snow microstructure properties). Additional information on the computation of the vertical stress and snow viscosity can be found in Reference VionnetVionnet and others (2012).

The evolution of the snow temperature profile is given by

(3)

where c _i is the heat capacity of ice, T _sn is the temperature of the snow layer, L _f F is the heat released by refreezing of water and L _f M is the heat flux related to melt. The last term on the right-hand side represents the penetration of net shortwave radiation. The effective thermal conductivity of the snow (W m^–1 K^–1) is computed using the formulation of Reference YenYen (1981):

(4)

where ρ _w is the density of water.

Vertical percolation of meltwater is gravitationally driven (bucket model formalism), and the water flow solution procedure starts from the top layer and proceeds downward (Reference VionnetVionnet and others, 2012). Crocus assumes homogeneous infiltration of water into firn. Vertical percolation of surface meltwater through a given layer is allowed when the computed liquid water content of the layer exceeds its maximum liquid-water-holding capacity, assuming that the layer can no longer freeze liquid water. The model’s default maximum liquid-water-holding capacity (W _{liq-_max}; kg m^–2), is expressed as 5% of the total pore volume (Reference VionnetVionnet and others, 2012):

(5)

where D is the snow layer thickness, LWC is the liquid water content, and ρ _w, ρ _i and ρ are the water, ice and snow layer densities, respectively. Under this parameterization, the model neglects the formation of capillary barriers such as thick impermeable ice layers. Instead, once an ice layer has reached its maximum liquid-water-holding capacity, the model forces meltwater to flow through it. A percentage of total pore volume smaller (larger) than 5% decreases (increases) liquid retention, which enhances (reduces) vertical percolation and refreezing of meltwater at depths. As Crocus is used primarily to simulate snowpacks that are warmer and moister than that which characterizes Devon Ice Cap, the model was tested for maximum liquid-water-holding capacities ranging between 1% and 6% of pore volume, with increments of 1%. Comparison between observed density profiles and profiles simulated using values between 1% and 3% of pore volume revealed a difference of 10% between simulated and observed density profiles. When values between 4% and 6% of pore volume were used, however, differences between observed and simulated density profiles increased to 15–17%. Given these results, and, in an effort to keep the value of the maximum liquid-water-holding capacity close to Crocus’s default value, results presented here are from simulations that used a maximum liquid-water-holding capacity of 3%.

2.2. Initial conditions

Firn density profiles for model runs were initialized using field observations from 30 April 2004 at Sites 1 and 2, and 30 April 2006 at HB 9–1 and HB 13–7 (Fig. 2; Table 1). Simulations ran continuously at a time step of 15 min until 1 May 2012 at all sites. Initial density profiles were derived from 20 and 18 m firn cores drilled in 2004 at Sites 1 and 2, respectively, and from 4.5 m cores drilled in 2006 at HB 9–1 and HB 13–7. Detailed analyses of these cores are presented in Reference Bezeau, Sharp, Burgess and GasconBezeau and others (2013). Density profiles were derived by dividing the mass of each core segment by its volume. Core segments were 0.04–0.25 m in length, and consisted of the relatively homogeneous material between density transitions along the core. Average density measurement error is estimated to be 40 kg m^–3. In the absence of direct observations from HB 9–1 and HB 13–7, initial firn density profiles for depths between 5 and 25 m were prescribed to firn density profiles measured in 2006 at those sites. Initial density profiles below the surface at HB 9–1 and HB 13–7 were set to increase from 525 kg m^–3 at 5 m depth to 560 kg m^–3 at 25 m depth, and from 625 kg m^–3 at 5 m depth to 650 kg m^–3 at 25 m depth, respectively. These profiles are based on profiles measured in the upper percolation area of the Greenland ice sheet (Reference Brown, Bradford, Harper, Pfeffer, Humphrey and Mosley-ThompsonBrown and others, 2012). Initial density profiles were also extended from 20 and 18 m depth at Sites 1 and 2, respectively, to a depth of 25 m by assuming an initial constant density of 550 kg m^–3. The effect of including deeper layers in the evolution of the density profiles was evaluated, but did not change the outcome of the simulations.

Table 1.

General information about the four field sites

In the absence of detailed firn temperature profiles for 2004, 16 day average spring 2012 firn temperature measurements were used to represent the end-of-winter temperature profiles. These were made with 10 k Ω negative temperature coefficient thermistors (accuracy ± 0.1ºC) located every 1 m to a depth of 15 m at Sites 1 and 2, and at 10 and 15 m depth at HB 9–1 and HB 13–7 (Reference Bezeau, Sharp, Burgess and GasconBezeau and others, 2013). At HB 9–1 and HB 13–7, spring 2012 firn temperature measurements were available only for depths of 10 and 15 m. At these two sites, vertical firn temperature profiles were extrapolated between the surface and 10 m below the surface assuming the same vertical temperature gradients as those derived from temperature measurements at Sites 1 and 2. At all sites, we assumed a constant firn temperature of –17.5ºC below 15 m. One 24 hour average firn temperature measurement made at 10 m below the snow surface at Site 1 in 2004 allowed us to determine that between 2004 and 2012 the firn temperature at 10 m depth at Site 1 increased by 3.8ºC, from –21.3 to –17.5ºC.

Empirical laws of snow metamorphism and snow compaction are dependent on firn temperatures and the resolution of the snow temperature profile. In an effort to quantify the effect of different firn temperature initializations on our simulations and assess the model’s sensitivity to initial firn temperature profiles, we tested the model with initial firn temperatures up to 3.8ºC cooler than the values used in the initiation of the simulations presented in this paper (i.e. the magnitude of the measured firn temperature warming at 10 m depth at Site 1 between 2004 and 2012). We found that this led to a difference of <1% in the simulated firn density profiles, while it did not affect the magnitude and direction of the changes in firn temperature profiles over time. Although 2012 firn temperature profiles were warmer than those of 2004, they remained well below the melting point of ice and firn, and a relatively short simulation (8 years) produced non-significant differences in model outcome between different initial firn temperature profiles. We therefore assumed that using the 2012 firn temperature profiles to initialize firn temperature profiles at all sites would not significantly alter the firn density profile evolution during the simulations. As simulations were initialized for 30 April 2004/2006 (prior to the beginning of the melt season), dry snow conditions were assumed for the start of all simulations.

2.3. Atmospheric forcing

Crocus requires inputs of air temperature, specific humidity, wind speed, shortwave and longwave radiation, precipitation rate and atmospheric pressure. Three-hourly surface data from the North American Regional Reanalysis (NARR; Reference MesingerMesinger and others, 2006) were used to force the model between 2004 and 2012 (Fig. 2). The NARR model uses the US National Centers for Environmental Prediction (NCEP) Eta Model (32 km/45 layers) together with the Regional Data Assimilation System (RDAS) to produce 3 hourly, daily and monthly data available from 1979 for 29 distinct pressure levels over North America at a horizontal resolution of 0.3˚ (32 km) at the lowest latitude (Reference MesingerMesinger and others, 2006). Local observations of wind and air temperature from automated weather stations (AWS) at Sites 1 and 2 between 2004 and 2010, and shortwave and longwave radiation budgets from a net radiometer at Site 2 between 2007 and 2010 were used to evaluate and adjust NARR fields prior to using them to drive simulations. NARR data adjustment was performed primarily for the melt season as there were data gaps during the winter seasons resulting from intermittent instrument failure and riming. Details of meteorological field data collection are provided in Reference Gascon, Sharp and BushGascon and others (2013a). Examples of the comparison between AWS and NARR 2 m air temperature, and of the meteorological forcings over one year at Site 2 are shown in Figures 3 and 4, respectively. NARR temperature fields were downscaled from Site 1 to HB 9–1, and Site 2 to HB 13–7, using an environmental lapse rate of 4.9˚C km^–1 (Reference GardnerGardner and others, 2009). Under conditions of flat and uniform topography, the wind field and the incoming shortwave radiation time series were assumed to be uniform between sites. Using these meteorological forcing data, the evolving albedo, surface fluxes and surface energy balance were calculated for each time step at each site using Crocus.

Fig. 3.

Comparison between Site 2 AWS (solid) and NARR (dotted) 2 m air temperature between May and September for (a) 2004, (b) 2005, (c) 2006, (d) 2007, (e) 2008, (f) 2009 and (g) 2010. The average correlation coefficient is 0.8. For reference, days 150 and 250 are 30 May and 7 September, respectively.

Fig. 4.

Example of the atmospheric forcing time series for 2010 at Site 2: (a) 2 m air temperature, (b) wind speed, (c) net shortwave (SW) and longwave (LW) radiation, (d) precipitation, (e) pressure and (f) specific humidity. For reference, days 150 and 250 are 30 May and 7 September, respectively.

Annual end-of-winter snowpack depth and density measurements from snow pits (Reference Bezeau, Sharp, Burgess and GasconBezeau and others, 2013) were used to adjust NARR precipitation values at all four sites. Based on linear scaling between NARR precipitation and field observations, we found that NARR underestimated precipitation by 30–50%, with the highest underestimation occurring at Site 1. Model simulations presented in this paper were forced with scaled NARR precipitation forcing.

2.4. Calculation of the mass percent difference

Simulated density profiles, and mass and thickness changes between 2004 and 2012 were compared to field observations. Model density profiles were compared with measured profiles from 4.5 m cores drilled in 2008 and 2011 (Site 1, HB 9–1 and HB 13–7), and from 15 m firn cores drilled in 2012 (all sites). Detailed descriptions of these cores are provided in Reference Bezeau, Sharp, Burgess and GasconBezeau and others (2013). The masses of firn and ice in each core were calculated by multiplying each layer’s thickness by its density, and summing the products down to 4.5 or 15 m depth. Identical calculations were applied to the simulated density profiles to obtain simulated masses. The mass percent differences between simulated and observed profiles were then calculated for each site and each year by taking the difference between simulated and observed masses, and dividing by the observed mass. Both simulated and observed density profiles were linearly extrapolated to uniform thickness layers of 0.05 m to facilitate comparison.

3.1. Model evaluation

3.1.1. Density profile

Table 2 shows the mean error (ME) and root-mean-square error (RMSE) for each site, from 8 years of simulations for the densities of the full firn columns. The average RMSE ranged between 115 kg m^–3 at Site 2, and 184 kg m^–3 at Site 1. The average ME at Site 1 (117 kg m^–3), on the other hand, was 1.4–2.8 times larger than the average ME at each of the three other sites, suggesting reduced model accuracy at Site 1 relative to the other sites. This can also be seen from the comparison between simulated density profiles throughout the simulations and observations (Figs 5 and 6). Crocus performed much less well at Site 1 than at the other sites, especially at the end of the simulations (Fig. 6).

Table 2.

Comparative statistics in terms of mean error (ME) and root-mean-square-error (RMSE) between modeled and observed density (kg m^–3) for the full firn column at the time of comparison

Fig. 5.

Simulated (solid) and observed (dash) density profile (kg m^–3) of the top 5 m of firn at the end of April for 2008, 2011 and 2012 at (a) Site 1, (b) HB 9–1 and (c) HB 13–7. Maximum water content capacity of 3% is used.

Fig. 6.

Simulated (solid) and observed (dash) density profile (kg m^–3) of the top 15 m of firn at the end of April for 2012 at (a) Site 1, (b) HB 9–1, (c) HB 13–7 and (d) Site 2. Maximum water content capacity of 3% is used.

The evolution of density profiles throughout the simulations suggests that Crocus progressively overestimated firn densities over time in the upper part of the profile, with larger overall overestimation occurring for the upper 4.5 m of firn at higher-elevation sites (Site 1 and HB 9–1) (Figs 5 and 6). The observed differences between density profiles were quantified as the difference between the simulated and observed mass of the full firn columns. The mass percent difference between model and observations for the upper 15 m of the firn for the period 2004–12 ranged from –2% at HB 9–1 to 16% at Site 1, with an average of 10% (Table 3). Across all sites, the 2012 simulations overestimated the mass of the upper 4.5 m of the firn (including winter snowpack) by an average of 15% relative to observations, but they underestimated the mass of the lower part of the firn column (4.5–15 m) by 5%. The simulated mass of the upper 4.5 m of the firn was also an overestimate relative to observations at all sites in all years.

Table 3.

Percent difference between modeled and observed mass per unit area (kg m^–2) of the firn column at time of comparison. Positive values represent model overestimation of mass per unit area compared to observations

Simulated and observed winter snowpack masses were first compared to determine whether differences between the simulated and observed winter snowpack alone could account for the differences in simulated and observed masses of firn. We found that the simulated winter snowpack mass was overestimated compared to observations, with an average percent difference ranging between 3% in 2011 and 5% in 2008, and an average percent difference of 4% for all years (Table 4). A larger simulated winter snowpack (average of 4%) implies that more mass was added to the modeled firn column than was the case in reality. Crocus’s average 4% overestimation of the winter snowpack mass relative to observations is, however, smaller than its 10% average overestimation of the mass of the entire column of firn (Tables 3 and 4). This suggests that the differences between simulated and observed mass and density profiles over the course of the simulation are not solely due to the model’s overestimation of winter snowpack. Although NARR temperature and wind data were corrected using field data, data gaps in our field observations during the winter may have introduced errors in the calculation of the density of fresh fallen snow (Eqn (1)). As a result, the wind-induced surface snow densification may not be represented accurately if wind speed, and thus density of surface snow, is erroneous. The underestimation of snow erosion by wind, and the percolation of meltwater to greater depths than Crocus simulates may also contribute to the calculated differences between simulated and observed mass and density profiles. Although mass is conserved in all firn processes, the representation of percolation as a homogeneous (rather than heterogeneous) process may lead to overestimation of mass near the surface as a result of too little percolation being simulated at depth relative to what would be the case with heterogeneous percolation, which can transport meltwater to greater depths by piping.

Table 4.

Percent difference between simulated and observed winter snowpack mass per unit area (kg m^–2) at time of comparison. Positive values represent model overestimation of mass per unit area compared to observations

In Crocus, water initially froze nearer the surface than in observations, creating a positive model density bias near the surface, and a negative density bias at depth. The assumption of homogeneous percolation implies that the capacity for refreezing meltwater in firn is a function of the total pore volume of the firn layer and of their thermal state. In reality, percolation is heterogeneous, which allows localized deep infiltration of surface melt and the formation of discontinuous ice layers, and implies that some parts of the pore volume are bypassed by water flow (Fig. 1). In the accumulation area, deposition of new snow over each winter added new pore volume surrounded by a matrix of snow crystals. This likely acted as a sink for percolating meltwater at the start of the following melt season. As Crocus models percolation as a homogeneous process, this creates conditions that favor near-surface freezing. Thus, modeling percolation as homogeneous, and as such, not considering deep heterogeneous percolation, appears to have resulted in too much mass being retained in the upper part of the firn profile relative to observations, and in insufficient mass accumulating at depths (Fig. 6).

Improved model performance at lower-elevation sites relative to higher-elevation sites (Figs 5 and 6) may suggest that homogeneous percolation is a more realistic assumption at sites where air temperature is higher, more surface melt is generated, and internal refreezing and ice accretion are more homogeneous. This pattern has been documented using repeat ground-penetrating radar surveys performed over the same transect between Sites 1 and 2 over the period 2007–12. These showed a strong signature of heterogeneous percolation such as pipes and ice lenses at high elevations (closer to Site 1), but more continuous horizontal ice layers at HB 9–1 and below (Reference Gascon, Sharp, Burgess, Bezeau and BushGascon and others, 2013b).

Crocus density profiles for all sites showed less vertical variability below 4.5 m depth than observed profiles (Fig. 6). The presence of alternating layers of firn and ice below 5 m is evident in observed density profiles and firn stratigraphy (Fig. 6), but is not seen in the simulated density profiles. Crocus allows vertical percolation of surface meltwater through a given layer only when it becomes saturated with respect to the total pore volume (Eqn (5)) and can no longer freeze liquid water. Under conditions of subzero winter firn temperatures, this implies that water-saturated layers would freeze during and/or at the end of the summer, automatically forming ice layers. Over the winter, new snow will bury ice layers formed at (or close to) the surface. During the subsequent summer, new ice layers will be created at (or close to) the surface on the snow that accumulated during the previous winter. If not all of the overwinter snowpack is turned into an ice layer by the end of the melt season, a layer of residual firn will be preserved between two ice layers, as can be observed in some firn cores (Fig. 7). Equation (5), however, implies that when an ice layer forms, its liquid-water-holding capacity is very small, and water is forced to percolate through the ice layer rather than pond or flow laterally above it. This is illustrated in our simulations by the progressive densification of the top layers of the firn column, and the gradual disappearance of residual firn over time (see animated evolution of the firn density profile at HB 9–1 between 2004 and 2010 at http://www.igsoc.org/hyperlink/13j209.avi). The video clearly shows the densification of the top 1–3 m of the firn column during the summer seasons, and the disappearance of residual firn (less vertical variability in density profiles) over time. This, combined with the representation of deep heterogeneous percolation as a homogeneous process, results in the inability to simulate the survival of residual firn bodies between ice layers.

Fig. 7.

Firn stratigraphy showing alternating firn and ice layers between 5 and 10 m depth at HB 13–7 at the end of April 2012: (a) 5–6 m, (b) 6–7 m, (c) 7–8 m, (d) 8–9 m and (e) 9–10 m. Photographs are infrared using a modified Fuji FinePix S9000 digital camera and an infrared filter of 850 nm to magnify contrast between snow textures (Reference Bezeau, Sharp, Burgess and GasconBezeau and others, 2013).

In order to evaluate the modeled firn density profiles at all simulated depths, we compared the observed initial density profiles to their simulated final density profiles to investigate whether firn density changes are observed over the full firn column (Fig. 8). As previously mentioned, the average rate of winter accumulation on Devon Ice Cap is 0.23 m a^–1. Over the 8 simulated years, this sum up to ~2 m w.e. Given an average firn density of 600 kg m^–3 over the first few meters, this accumulation represents the upper ~3.5 m of firn in the final 2012 simulated density profiles. This implies that the first ~3.5 m of firn of the 2012 density profiles at each site have been entirely simulated by Crocus. Below that depth, density profiles initially prescribed are modified at each time step through firn compaction, snow metamorphism and internal refreezing. As seen in Figure 8, the largest differences between the observed initial and simulated final density profiles occur in the upper layers, from the top 4.5 m of firn at HB 9–1, to the top 9 m at Site 2. However, the depth of the zone of greatest changes exceeds the depth of entirely modeled firn at all sites, suggesting that the simulated evolution of firn layers extends into the initially prescribed part of the firn. At HB 9–1, very little change was observed between the 2006 observed initial and 2012 simulated final density profiles. We attribute this to the limited availability of initial observations for that site, which resulted in a larger section of the profile being prescribed on the basis of firn density profiles from the upper percolation area of the Greenland ice sheet (Reference Brown, Bradford, Harper, Pfeffer, Humphrey and Mosley-ThompsonBrown and others, 2012). At HB 13–7, a large section of the initial density profile was also prescribed assuming conditions similar to those on the Greenland ice sheet. Unlike for HB 9–1, however, Crocus was able to simulate firn densification through snow compaction below 5 m depth at HB 13–7 (Fig. 8c). At Site 2, the initial, observed and 2004 and final simulated 2012 density profiles are very similar (Fig. 8d). Site 2 is located at the equilibrium line (Reference Gascon, Sharp, Burgess, Bezeau and BushGascon and others, 2013b), where the net annual mass gain is zero. As such, we expect little vertical displacement of the initial firn stratigraphy. The final modeled 2012 stratigraphy actually appears to have shifted upward, suggesting model overestimation of melt over the 8 year simulation. The final modeled firn densities between 9 and 15 m depth for 2012 are, however, larger than the observed densities from 2004, illustrating some firn compaction over time.

Fig. 8.

Comparison between observed initial (dash) and simulated final (solid) density profiles (kg m^–3) at the end of April for (a) 2004 and 2012 at Site 1, (b) 2006 and 2012 at HB 9–1, (c) 2006 and 2012 at HB 13–7 and (d) 2004 and 2012 at Site 2.

3.1.2. Temperature profile

Densification is a function of rates of snow accumulation, metamorphism, compaction, and pore volume reduction by refreezing. The latter is a function of meltwater percolation and snow/firn temperature. After the parameterization of available pore volume and the vertical flow of meltwater, firn temperature is the major variable that determines whether meltwater will refreeze in a given layer or not. As a result, evaluation of the modeled temperature field may provide additional explanations of the differences between the observed and modeled firn density profiles described in the previous subsection. As our dataset includes full firn temperature profiles at all sites for one spring season only, 2012, and one temperature measurement at a depth of 10 m at Site 1 in 2004 (see Section 2.2), only a limited evaluation of the accuracy of the model’s representation of the processes that control temperature evolution was possible. The simulated firn temperature profiles between 2004 and 2012 are first presented, followed by a brief comparison of observed and simulated profiles.

Figure 9 shows the changes in the temperature profiles that occurred between the initial conditions prescribed for spring 2004, and the modeled final state in spring 2012. We attribute the cooling of the top 1–2 m in 2012 relative to 2004 to year-to-year variations in the air temperature at the time of comparison (end of April). Although the 2012 melt season was warmer than the 2004 melt season, the midApril/mid-May average air temperature on Devon Ice Cap was lower in 2012 (–21.5˚C at Site 2) than in 2004 (–13.1˚C at Site 2). Below 2 m, the simulated warming of firn between 2004 and 2012 reflects the combined effect of warming due to the increase in mean annual air temperature and the release of latent heat during refreezing of meltwater. Larger warming at lower-elevation sites, with the exception of HB 9–1, illustrates a decrease in the modelled magnitude of firn warming with increasing elevation, in agreement with results reported by Reference Bezeau, Sharp, Burgess and GasconBezeau and others (2013). We speculate that the cooling at HB 9–1 could be due to a combination of underestimation of internal refreezing, as shown by the 2% underestimation of the mass of the firn column in 2012 compared to observations (Table 3), and modelling errors.

Fig. 9.

Difference of the temperature profiles between observed initial conditions in spring 2004, and modeled final state at the end of April 2012 at Site 1, and Site 2, and between observed initial conditions in spring 2006, and modeled final state in spring 2012 at HB 9–1 and HB 13–7.

The modeled firn temperature at 10 m depth at Site 1 increased by 0.7˚C between 2004 and 2012. This is an underestimation compared to the observed 3.8˚C increase over the same time period. As refreezing of meltwater releases energy and warms the surrounding firn, an underestimation of modeled firn temperature compared to observations at 10 m depth suggests an underestimation of meltwater refreezing at that depth. We attribute this to the model representation of percolation as a homogeneous process, which, as shown in the previous section, limits the amount of meltwater percolation and refreezing at depth. The lack of observed firn temperatures for both 2004 and 2012 at the other three sites prevents us from evaluating Crocus’s ability to represent the evolution of firn temperatures at all sites. Additional field observations of firn temperatures at different depths over multiple years are thus required to rigorously evaluate Crocus’s simulations of the thermal evolution of the snowpack and firn.

3.2. Implications for mass-balance modeling

A vertical distribution of refreezing in the firn that is different from what is observed may result when Crocus is run for climate-warming scenarios. In our simulations, large amounts of meltwater refroze homogeneously close to the surface whereas, in reality, part of the surface meltwater appears to have percolated heterogeneously to greater depths. Observations show that heterogeneous percolation can lead to vertical transport of meltwater via saturated pipes to depths of at least 10 m in Greenland (Reference Humphrey, Harper and PfefferHumphrey and others, 2012), and similar piping signatures have been observed on Devon Ice Cap. Under a rapid climate-warming scenario, melt can potentially spread quickly to large areas where it was previously very limited or non-existent. As a result, the upper boundaries of the wet snow and percolation zones may migrate rapidly to higher elevations. This has been observed on Devon Ice Cap since 2005 (Reference Bezeau, Sharp, Burgess and GasconBezeau and others, 2013; Reference Gascon, Sharp, Burgess, Bezeau and BushGascon and others, 2013b). Similar migration of the snow facies on the Greenland ice sheet may be observed in the near future since recent extreme melt events (e.g. Reference NghiemNghiem and others, 2012) suggest that the dry snow zone there may disappear within the next decade (Reference McGrath, Colgan, Bayou, Muto and SteffenMcGrath and others, 2013). Without consideration of heterogeneous percolation, however, models could overestimate the extent of near-surface refreezing, leading to an overestimation of the rate of migration of the different snow facies to higher elevations. As the migration of snow facies to higher elevations is associated with a reduction of pore space in firn that will ultimately lead to an increase in surface runoff (Reference Gascon, Sharp, Burgess, Bezeau and BushGascon and others, 2013b), Crocus’s current refreezing parameterizations (and those of any other snowpack model that considers only homogeneous percolation) may lead to an overestimation of surface runoff because they result in too much ice accumulating close to the surface.

In a warming climate, both environmental conditions and firn stratigraphy will change. Rain is likely to be a more common form of precipitation over glaciated regions in summer as a result of higher atmospheric temperatures. Unlike snow, rainwater is already warmer than the freezing point of water, does not require additional energy to cause melt, and may even enhance surface melt through heat transfer to the snow via the albedo feedback. Such a change in precipitation type would lead to an increase of water input to the surface of the snow, which would increase percolation and refreezing of water in firn. Until recently, thick quasi-impermeable ice layers had not been observed in any ice sheet, making it impossible to know whether or how meltwater might percolate through them. Reference Gascon, Sharp, Burgess, Bezeau and BushGascon and others (2013b) argued that the presence of thick ice layers on Devon Ice Cap reduces pore space in firn and can enhance surface runoff. Crocus, however, allows meltwater to percolate through ice layers. This model behavior partially counteracts the overestimation of mass closer to the surface that results from the representation of percolation as a homogeneous process, and may create compensating errors in the model formulation. In addition, Crocus uses a bucket model to simulate meltwater percolation. In comparison to a wetting front model, where the wetting front moves with a certain velocity, Crocus does not have the ability to account for the variable time lag between surface infiltration and percolation at depth. Recent developments, such as the experimental analysis of preferential water flow in a dry snowpack (Reference Katsushima, Yamaguchi, Kumakura and SatoKatushima and others, 2013), the implementation of this process in physically based snowpack models (Reference Wever, Fierz, Mitterer, Hirashima and LehningWever and others, 2014) and consideration of the effects of impermeable ice layers on meltwater percolation in firn (e.g. Reference Østby, Schuler, Hagen, Hock and ReijmerØstby and others, 2013) will contribute to the ongoing improvement of the representation of heterogeneous percolation in snowpack models, and its effects on firn stratigraphy.