How well is firn densification represented by a physically based multilayer model? Model evaluation for Devon Ice Cap, Nunavut, Canada

Abstract

Evaluation of how accurately snowpack models can capture hydrological processes in firn is critical to determining how well they will simulate future glacier mass-balance changes. Here we compare simulations using the Crocus snowpack model with the evolving firn stratigraphy recorded in 14 cores drilled at four elevations in the accumulation zone of Devon Ice Cap, Canada, during the 2004–12 period of rapid summer warming. Simulations were forced with a combination of surface observations and reanalysis data. Simulations resulted in positive model bias in near-surface density, and negative bias in density at depth compared to observations. Results point to the importance of incorporating heterogeneous percolation in firn in order to improve the representation of meltwater flow, better reproduce observed firn density and temperature profile evolution, and improve simulations of glacier mass balance during periods of climate warming.

1. Introduction

Mass loss from polar ice sheets and ice caps contributes significantly to eustatic sea-level rise, and has been increasing since the beginning of the 21st century (e.g. Hanna and others, 2008; Gardner and others, 2011, 2013; Rignot and others, 2011; Sharp and others, 2011; Svendsen and others, 2013). Meltwater storage in firn (compacted snow that has survived at least one melt season) through refreezing and retention as liquid currently limits the amount of surface melt from polar ice sheets and ice caps that reaches the oceans (Harper and others, 2012). Meltwater can infiltrate firn by either homogeneous or heterogeneous percolation. Homogeneous percolation consists of uniform vertical flow of meltwater, while heterogeneous percolation involves vertical flow of meltwater concentrated in specific areas, typically referred to as piping (e.g. Pfeffer and others, 1991) (Fig. 1). Meltwater storage in firn by refreezing is especially important at higher elevations in the accumulation area, above the superimposed ice zone, as not all the previous winter snowpack melts during the summer and surface meltwater can percolate into and refreeze within cold firn below the surface (e.g. Humphrey and others, 2012; Bezeau and others, 2013).

Fig. 1. Examples of heterogeneous percolation and refreezing: (a) for the top 3 m below the surface at Site 1 in 2011, and (b) at 8–9 m depth at Site 1 in 2011. Photographs are infrared, using a modified Fuji FinePix S9000 digital camera and an infrared filter of 850 nm to magnify contrast between snow textures (Bezeau and others, 2013).

A projected increase in Arctic air temperatures of up to 2^ºC over the next 20 years (Meehl and others, 2014) would initially increase the amount of percolating meltwater that refreezes in the shallow subsurface layers of ice caps. Such changes could raise the near-surface temperature of the ice caps through latent heat release, and significantly increase rates of firn densification (Zdanowicz and others, 2012; Bezeau and others, 2013). Warming of the firn could promote further surface melt and thus internal refreezing. An increase in mean summer air temperature of ~ 4^ºC was measured on Devon Ice Cap, Nunavut, Canada, between 2004 and 2010 (Gascon and others, 2013a). This warming was associated with the rapid densification of the first 4 m of firn below the surface of the accumulation area of the ice cap, and with migration of the upper boundaries of the different snow facies zones to higher elevations after 2005 (Bezeau and others, 2013; Gascon and others, 2013b). This rapid densification occurred by heterogeneous refreezing of meltwater, and is progressively reducing the water storage potential and refreezing capacity of the firn layer of the ice cap. This will eventually lead to an increase in the fraction of meltwater that runs off, unless further internal warming of the firn re-increases firn pore space and its water storage potential.

Evaluation of how quickly this transition towards increased surface runoff will occur requires accurate modeling of firn processes. Such processes are currently represented by parameterizations implemented in snowpack schemes of land surface models (e.g. Vionnet and others, 2012) that are coupled with climate models (e.g. Reijmer and others, 2012; Van Angelen and others, 2013). Obleitner and Lehning (2004) used a physically based snow model to successfully reproduce the formation of superimposed ice in a glacier in Svalbard, but their model evaluation did not extend to the higher portions of the accumulation area where heterogeneous percolation is most important. In the absence of direct observations from all snow facies of the accumulation areas of glaciers and ice caps (percolation, wet-snow and superimposed ice zones), refreezing parameterizations have been validated by comparing results derived using different parameterizations (Reijmer and others, 2012). Although the latter study showed that the six tested refreezing parameterizations agreed with each other (using the parameterization implemented within the Regional Atmospheric Climate Model (RACMO2) as a benchmark), none of them accounted explicitly for heterogeneous percolation of meltwater. Without evaluation against observations, it is impossible to determine whether these parameterizations accurately represent the magnitude and distribution of refreezing. Since refreezing parameterizations are included in models that are being used to predict future glacier mass balance and sea-level change, it is critical to determine whether they are capable of simulating changes in firn processes that are as large and rapid as those that have been observed in the accumulation area of Devon Ice Cap since 2005.

Here we use the physically based multilayer snowpack model Crocus to simulate changes in firn density and temperature profiles over all firn facies of the accumulation area of Devon Ice Cap during the period of rapid summer warming between 2004 and 2012. Model results were compared with the evolving firn stratigraphy recorded in 14 firn cores drilled at four sites in the accumulation zone of the ice cap over the study period (Fig. 2).

Fig. 2. Map of the Devon Ice Cap field sites: Site 1 (1800 m a.s.l.), HB 9–1 (1610 m a.s.l.), HB 13–7 (1490 m a.s.l.) and Site 2 (1400 m a.s.l.). In 2004, Site 1 and HB 9–1 were located in the percolation zone of Devon Ice Cap, while HB 13–7 and Site 2 were located in the wet snow zone. By 2012, Site 1 was still located in the percolation zone, while HB 9–1 was in the wet snow zone, and HB 13–7 and Site 2 in the superimposed ice zone. Squares denote locations of weather stations, net radiometers, and firn-core sites, and circles denote locations of firn-core sites. The black triangle denotes the location of the NARR model results used in this paper.

2. Methods

2.1. Model description

We used the one-dimensional multilayer snowpack model Crocus (Vionnet 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 (Masson 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 (Colgan 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 (Vionnet 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 Vionnet 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 Yen (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 (Vionnet 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 (Vionnet 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 Bezeau 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 (Brown 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 (Bezeau 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; Mesinger 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 (Mesinger 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 Gascon 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 (Gardner 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 (Bezeau 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 Bezeau 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. Results and Discussion

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 (Gascon 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 (Bezeau 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 (Brown 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 (Gascon 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 Bezeau 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 (Humphrey 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 (Bezeau and others, 2013; Gascon 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. Nghiem and others, 2012) suggest that the dry snow zone there may disappear within the next decade (McGrath 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 (Gascon 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. Gascon 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 (Katushima and others, 2013), the implementation of this process in physically based snowpack models (Wever and others, 2014) and consideration of the effects of impermeable ice layers on meltwater percolation in firn (e.g. Ø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.

4. Conclusion

We used a combination of data from firn cores and from NARR to force simulations with the Crocus snowpack model for Devon Ice Cap for the period 2004–12. We compared model simulations with the evolving firn stratigraphy recorded in 14 cores drilled at four elevations in the accumulation zone of Devon Ice Cap over the same period. Simulations with Crocus result in a mass bias of between –2% and +16% relative to observations in the upper 15 m of the firn. Inspection of the differences between simulated and observed density profiles showed that Crocus overestimated the mass of the top 4.5 m of firn and ice at all sites, but underestimated the mass of firn below 4.5 m depth. Annual winter snow accumulation continually added new pore volume close to the surface. As Crocus represents percolation as a homogeneous process, this created conditions that favoured near-surface freezing and insufficient refreezing at depth. However, Crocus’s parameterization of the permeability of ice layers forces meltwater to percolate through them, preventing the build-up of thick impermeable ice layers. Given that percolation in the upper part of the accumulation area of ice sheets tends to be heterogeneous rather than homogeneous, and given the potential importance of thick ice layers in restricting vertical meltwater drainage, it is important for models to simulate both heterogeneous and homogeneous percolation in order to accurately reproduce both firn density and temperature profile evolution and mass balance through periods of climate warming.