Calibration and evaluation of a high-resolution surface mass-balance model for Paakitsoq, West Greenland

Abstract

Modelling the hydrology of the Greenland ice sheet, including the filling and drainage of supraglacial lakes, requires melt inputs generated at high spatial and temporal resolution. Here we apply a high spatial (100 m) and temporal (1 hour) mass-balance model to a 450 km² subset of the Paakitsoq region, West Greenland. The model is calibrated by adjusting the values for parameters of fresh snow density, threshold temperature for solid/liquid precipitation and elevation-dependent precipitation gradient to minimize the error between modelled output and surface height and albedo measurements from three Greenland Climate Network stations for the mass-balance years 2000/01 and 2004/05. Bestfit parameter values are consistent between the two years at 400 kg m^-3, 2^°C and +14% (100 m)^-1, respectively. Model performance is evaluated, first, by comparing modelled snow and ice distribution with that derived from Landsat-7 ETM+ satellite imagery using normalized-difference snow index classification and supervised image thresholding; and second, by comparing modelled albedo with that retrieved from the MODIS sensor M0D10A1 product. Calculation of mass-balance components indicates that 6% of surface meltwater and rainwater refreezes in the snowpack and does not become runoff, such that refreezing accounts for 31% of the net accumulation.

1. Introduction

The Greenland ice sheet (GrIS) is the largest terrestrial permanently snow- and ice-covered area in the Northern Hemisphere and is highly sensitive to climate change (Box and others, 2006; Fettweiss, 2007). It is important to assess the impact of climate change on the GrIS as the temperature rise in high northern latitudes is strongly correlated with global warming, has increased at almost twice the global average rate over the past 100 years and is likely to continue at this rate into the future (Solomon and others, 2007; AMAP, 2009). At Aasiaat, on Greenland’s west coast, 2010 was the warmest year since records began in 1951, with records also set for winter, spring, May and June temperatures (Box and others, 2009;Cappellen, 2010;Tedesco and others, 2011).

Although the GrIS was thought to have been in near balance with the colder climate of the 1970s and 1980s, it has responded rapidly to post-1990 warming (Hall and others, 2008; Hanna and others, 2008; Rignot and others, 2008). The GrIS is currently experiencing a net mass loss as a result of increased wastage along its margins (Alley and others, 2010;Chen and others, 2011;Zwally and others, 2011). This is partly due to the acceleration of outlet glaciers (Rignot and Kanagaratnam, 2006; Howat and others, 2011) but also to increased surface meltwater runoff (Mote, 2007; Hanna and others, 2008;Van den Broeke and others, 2009) and associated increase in the size of the melt area (Tedesco, 2007;Fettweiss and others, 2011a).

Previous modelling studies have shown that a 1°C rise in surface air temperature over Greenland produces 20–50% more melt (Oerlemans, 1991;Braithwaite and Olesen, 1993;Ohmura and others, 1996;Janssens and Huybrechts, 2000; Hanna and others, 2005), so a predicted rise in air temperature of 2–5°C would approximately double total melt (Mernild and others, 2008). Although feedbacks due to increased surface melt are thought to be complex and the sign of the forcing is unknown, increased surface melt is likely to result in reduced albedo and increased incident energy absorption, leading to further melt (Tedesco and others, 2011). Thus, modelling the surface mass balance of the GrIS is a key component in predicting how the ice sheet will respond to projected future climate change. To predict detailed changes along specific margins, or of particular outlet glaciers, high spatial resolution models are required.

High spatial resolution surface mass-balance models are also required for studies concerned with calculating the delivery of surface meltwater to the ice-sheet bed via crevasses and moulins (Catania and others, 2008;Das and others, 2008; Van de Wal and others, 2008). It is thought that the seasonal filling and subsequent drainage of supraglacial lakes on the GrIS (Box and Ski, 2007;McMillan and others, 2007;Sneed and Hamilton, 2007;Selmes and others, 2011) may play a key role in linking the surface melt signal to ice motion by supplying the volume of water needed to propagate crevasses to the base of the ice sheet (Alley and others, 2005;Van der Veen, 2007;Das and others, 2008). Rapid supraglacial lake drainages may be accommodated by the subglacial drainage system through temporary spikes in subglacial water pressure (Schoof, 2010). These increases in subglacial water pressure reduce basal friction and drive transient ice-sheet accelerations (Iken and Bindschadler, 1986;Hooke and others, 1989;Mair and others, 2002).

The direct motivation for our study comes primarily from a desire to model accurately the filling and therefore subsequent drainage potential of supraglacial lakes, and the delivery of meltwater to moulins across the GrIS. To this end, it is necessary to have a good representation of the ice- sheet surface in the form of a high-resolution, accurate digital elevation model (DEM), and to model accurately surface melt and runoff rates, and the subsequent routing of water across the ice-sheet surface at a high (hourly) temporal resolution. These requirements are particularly critical since it is the variability in magnitude and timing of meltwater inputs to the subglacial drainage system, rather than the total water volume, which is thought to have the greatest influence on subglacial water pressures and therefore ice velocities (Schoof, 2010;Bartholomew and others, 2011; Colgan and others, 2011).

For our study, we model melt and runoff for the Paakitsoq region, West Greenland, using a surface energy-balance (SEB) model coupled to a subsurface model of the upper snow/firn/ice layers. The model is of high resolution with a gridcell size of 100m and is forced using a full range of meteorological variables, predominantly from the JAR1 Greenland Climate Network (GC-Net) automatic weather station (AWS) (Steffen and Box, 2001). Coastal precipitation data are obtained from the ASIAQ Greenland Survey.

For the sake of simplicity and to reduce computational expense, our melt/runoff model is driven by extrapolating the locally measured data across the model surface using elevation gradients, rather than, for example, downscaling the output from climate reanalyses or a regional climate model (RCM). This is justified because our model domain is small (450km²) and narrow (7.5 km), and there is a good sampling by AWSs.

We calibrate the model by comparing modelled and measured point data of surface height change and albedo changes at the GC-Net stations JAR1, JAR2 and Swiss Camp (Steffen and Box, 2001). Subsequently, we evaluate model performance in two ways. First, we compare the modelled snowline positions at various stages of the summer to the snowline positions delineated from Landsat-7 Enhanced Thematic Mapper Plus (ETM+) satellite imagery. Second, we compare the modelled daily albedo over the model domain with the daily surface albedo retrievals from the NASA Terra platform Moderate Resolution Imaging Spectroradiometer (MODIS) sensor M0D10A1 product.

We focus our study on mass-balance years 2000/01 and 2004/05, as during these years there are: (1) nearly continuous meteorological data, so only limited temporal interpolation is necessary;(2) relatively good albedo and surface melt data for model calibration;and (3) good availability of cloud-free Landsat and MODIS imagery during the summer months for model evaluation.

The model output of meltwater runoff per hour for each model gridcell will form the input to a supraglacial hydrology model which routes water: (1) in a saturated layer at the base of the snowpack;or (2) across exposed ice, to moulins (some of which will become active after the filling, then drainage, of supraglacial lakes). Ultimately, these models combined will provide input to a semidistributed, physically based subglacial hydrology model. These aspects of model development form the basis of ongoing work.

2. Study Site

The Paakitsoq region is defined here as the ~2300km² area of West Greenland, northeast of Jakobshavn Isbræ (Fig. 1). We choose this region because of the availability of meteorological data from the three GC-Net stations JAR1, JAR2 and Swiss Camp (Steffen and Box, 2001) and coastal precipitation data from ASIAQ Greenland Survey station 437 (190ma.s.l., ~4km west of the ice margin) (Fig. 1);(2) a relatively high-resolution bed DEM necessary for subsequent subglacial modelling; and (3) proglacial stream discharge data measured at ASIAQ station 437 for subsequent calibration of the full model. We focus on a ~7.5 km wide strip which extends diagonally from southwest to northeast across the Paakitsoq region and includes the GC-Net and ASIAQ stations (Fig. 1). The total area of this strip is 450 km².

Fig. 1. Map of the study site. The Paakitsoq region is delineated by the red box. The Landsat-7 ETM+ image behind is dated 7 July 2001. The strip on which we focus our study is outlined in black. The red triangle marks the ASIAQ 437 precipitation and gauging station. Green boxes are example areas of: (a) snow only; (b) ice only; and (c) snow and ice chosen during the thresholding procedure in Section 4.2. These three boxes correspond to the three example histograms of reflectance values shown below the map, where the y-axes are normalized to the maximum number of cells of a given brightness within each sample area, and the x-axes show brightness numbers.

The Paakitsoq region has predominantly land-terminating ice along its western margin, with a relatively small tidewater glacier that calves into a side arm of Jakobshavn Fjord. Ice elevation varies from near sea level at the western margin to ~1200 m a.s.l. at the easterly edge of the area. Ice flow in the region is predominantly westward, although towards the south it is influenced by Jakobshavn Isbræ and flow occurs in a more southwesterly direction (Mottram and others, 2009;Joughin and others, 2010). The recent retreat of Jakobshavn Isbræ may also be influencing some of the recently observed thinning across the area (Joughin and others, 2004; Krabill and others, 2004).

Previous studies in the region include the mapping of ice surface topography based on panchromatic aerial photography obtained on 10 July 1985 (Thomsen, 1986;Thomsen and others, 1988), and mapping of surface and bed topography using airborne laser altimetry and radar (Mottram and others, 2009). Sohn and others (1998) used a variety of satellite images to monitor the terminus position, and the ice- sheet margins to the north and south, of Jakobshavn Isbræ. To evaluate the potential for hydroelectric power production in the Paakisoq region, Ahlstrom and others (2008) assessed the future availability of meltwater based on a scenario run of a RCM. McMillan and others (2007) investigated surface lake development and drainage over a large part of the Paakitsoq region using a combination of Landsat-7 and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images from 2001, whereas Box and Ski (2007) focused more specifically on calculating temporal changes in depths and volumes of a few lakes through the use of a supraglacial lake-depth retrieval function, based on the correspondence between MODIS reflectance and water depth measured during raft surveys. Several studies in the region have linked glacier acceleration with increases in surface meltwater inputs to moulins, suggesting that hydrau- lically induced basal sliding might explain short-term and seasonal variability in ice velocities (e.g. Zwally and others, 2002;Das and others, 2008;Joughin and others, 2008;Price and others, 2008). Catania and others (2008) established that some, but not all, moulins are associated with surface lakes and may form after lake drainage via hydrofracture. Some moulins were found to be persistent englacial features throughout the winter, probably remaining active over several years (Catania and Neumann, 2010).

3. Surface Mass-Balance Model

The model is based on that of Rye and others (2010) and consists of three coupled components: (1) an SEB model that calculates the energy exchange between the glacier surface and the atmosphere;(2) a subsurface model, simulating changes in temperature, density and water content in the snow, firn and upper ice layers, and hence refreezing and net runoff; and (3) an accumulation model.

A DEM of the ice-sheet surface is required to spatially distribute meteorological data and compute slope angles and aspects as well as any topographic shading. Here we use the ASTER global DEM (GDEM) which has a nominal grid size of 30 m (http://asterweb.jpl.nasa.gov/gdem.asp). We checked the original GDEM for obvious artefacts in the area and none were found. The GDEM quality files for the Paakitsoq region show ASTER stacking numbers here lie between 8 and 12, yielding an accuracy of ±18.2 m at <500 m elevation, and ±13.8 m at >500 m (MacFerrin, 2011). The original data were smoothed with a 6 x 6 cell median filter to remove small-scale noise, then resampled to 100 m resolution using bilinear interpolation.

The model is forced using meteorological observations from the JAR1 GC-Net station where data are available. These include hourly measurements of incoming global shortwave radiation, air temperature, relative humidity and wind speed at 2 m above the ice surface. Due to instrument failure, shortwave radiation data are missing from 1 January to 23 May 2001, so data from Swiss Camp for this period were used instead. Daily precipitation totals (mmw.e.), collected at the ASIAQ station, were divided by 24 to produce hourly totals for model input on days with recorded precipitation. As incoming longwave radiation data were not available for the Paakitsoq region, incoming longwave radiation data were calculated using parameterizations based on the work of Konzelmann and others (1994). Further details of the SEB and subsurface models and the accumulation routine are described in Sections 3.1, 3.2 and 3.3 respectively.

3.1. Surface energy-balance model

The SEB model is based on that originally developed by Arnold and others (1996) and updated by Brock and others(2000), Arnold (2005), Arnold and others (2006) and Rye and others (2010). Here we focus primarily on describing adaptations that have been made to the Rye and others (2010) model for application to the Paakitsoq region.

The model determines the total SEB from five main components:

(1)

where Q _M is the energy available for melt (if surface temperature is at melting point), SW↓ is incoming shortwave radiation, SW↑ is outgoing shortwave radiation (i.e. (1 – surface albedo) * SW↓), LW↓ is incoming longwave radiation, LW" is outgoing longwave radiation, SHF is the sensible heat flux, LHF is the latent heat flux, and GHF is the ground heat flux in the snow or ice (Cuffey and Paterson, 2010).

Shortwave radiation and albedo

Due to lack of detailed cloud-cover records, we follow Oerlemans (1991b), and subsequently Arnold and others (1996), in assuming that diffuse radiation from the sky is one-fifth of the measured global shortwave radiation in all cases. Following the method of Arnold and others (1996), these hourly fluxes of direct and diffuse radiation are then modified at each model gridcell for the local terrain conditions in the case of direct radiation (namely slope angle, aspect, and shading due to surrounding gridcells, although the latter effect is very small over our domain and negligibly affects our results), and for the sky-view factor in the case of diffuse radiation.

Following Greuell and Konzelmann (1994), the surface albedo, α, is calculated as a linear function of the density of the uppermost subsurface grid element, ρ _top:

(2)

where α _snow, ρ _snow, α _ice and ρ _ice are fresh snow albedo, fresh snow density, ice albedo and ice density respectively. However, whereas Bassford (2002) and Rye and others (2010) denote the top 10 cm of snow and/or ice as the uppermost subsurface grid element, we follow Ettema and others (2010) and base our calculations on the top 5 cm of the subsurface grid so that small snowfall events have a more significant short-term influence on the calculations of surface albedo. In the parameterization, changes in surface density represent changes in grain size that affect the surface albedo (Greuell and Konzelmann, 1994). The value for fresh snow albedo is set as 0.82 which is the average snow albedo value calculated from the GC-Net data at JAR1, JAR2 and Swiss Camp from 2000 to 2009 and is also consistent with the snow albedo value used by Ettema and others (2010). If no snow cover exists, the albedo is set equal to the ice albedo value. The ice albedo, which is assumed constant in time and space, is set at 0.48 which is the average albedo measured at the three GC-Net stations from 2000 to 2009 during periods which we assume are snow-free. This value is slightly higher than the value (0.45) used by the MAR model of Fettweis and others (2011b), but is lower than the value (0.50) used by Ettema and others (2010). The value for fresh snow density is set during model calibration (Section 4.1). Snow cells within the model are converted to ice once the density exceeds that of superimposed ice (800 kg m^-3; Wright, 2005). A value of 910kgm^-3 is assumed for ice density (Cuffey and Paterson, 2010).

Longwave radiation

The net longwave radiation is the calculated sum of incoming flux from the sky (LW↓) and the outgoing radiation emitted by the glacier surface (LW↑). Since LW# is strongly dependent on cloud cover, but because we have no cloud- cover data for the Paakitsoq region, it was first necessary to calculate hourly values of cloud amount. We generated an empirical equation to calculate hourly mean cloud amounts for the Paakitsoq region through the combined use of (1) an existing parameterization for calculating LW↓ based on cloud cover (Konzelmann and others, 1994), and meteorological data (including measured LW↓) collected in 2009 from a temporary station on Russell Glacier near Kangerlussuaq, ~300km south of the Paakitsoq region (67°7.240⁰ N, 49°2.040’ W).

For skies that are completely overcast due to low clouds, LW↓ is primarily determined by the temperature of the cloud base (Oke, 1987). For completely overcast skies, therefore, the parameterization of LW↓ should not contain clear-sky emittance. Konzelmann and others (1994) satisfy this condition by describing emittance as the weighted mean of clear-sky and completely overcast emittances:

(3)

where n is cloud amount (n =1 for a completely cloud- covered sky, n = 0 for clear-sky conditions), E_cs is the clear-sky emissivity, E_oc is the emittance of a completely overcast sky, σ is the Stefan-Boltzmann constant (5.67 x 10^-8 Wm^-2 K^-4), T is temperature (K) 2 m above the ice surface and p is a constant.

Konzelmann and others (1994) relate clear-sky emissivity to measured vapour pressure, e, and T at Swiss Camp, through a modified version of Brutsaert’s (1975) equation applied to their data, yielding

(4)

Konzelmann and others (1994) obtained empirically the emissivity of a completely overcast sky, E_oc = 0.952, and the coefficient, p = 4, from hourly means of LW↓, T, e and instantaneous observations of n at Swiss Camp. Thus, the equation for incoming longwave radiation reads

(5)

Although the high power (p = 4) resulted from site-specific cloud climatology, where low cloud amounts were caused by high clouds in the sky, and high cloud amounts by low clouds in the sky, this expression has been used successfully at other sites in Greenland (e.g. Greuell and Konzelmann, 1994;Zuo and Oerlemans, 1996).

Using hourly means of LW↓, T and e measured at the Russell Glacier meteorological station from 26 June to 31 August 2009, hourly mean values for n over the glacier were calculated using Eqn (5). Since we have no measurements of n for the Paakitsoq region but we might expect there to be a good relationship between n and relative humidity, RH, for which we do have measurements, we use the first half of the available data (i.e. 26 June-29 July 2009) from the Russell Glacier dataset and investigate the relationship between n and RH for those data. We obtain

(6)

As an independent test of this relationship, mean hourly values of n over Russell Glacier for the second half of the available data (i.e. 30 July-31 August 2009) were calculated using Eqn (6), and were then used in Eqn (5) to calculate the mean hourly LW↓ fluxes. The root-mean-square error (RMSE) between the calculated and measured average daily LW# fluxes is 19.8 Wm^-2. In comparison, if we use a constant value of n = 0.6 (Konzelmann and others, 1994) the RMSE between the calculated and measured average daily LW↓ fluxes is 51.2 Wm^-2. We appreciate that the datasets used for both parameterizing and testing the relationship between n and RH are relatively short, but the analysis does suggest that in the absence of an alternative approach, RH provides a useful proxy for cloud amount on the GrIS and that the Konzelmann and others (1994) parameterization is a robust one. We therefore use this approach to determine hourly LW↓ at Paakitsoq.

The outgoing longwave radiation emitted from the glacier surface is computed as a function of surface temperature, T_s (i.e. the temperature of the top gridcell), and the StefanBoltzmann constant, σ:

(7)

Turbulent heat fluxes

The sensible and latent turbulent heat fluxes are calculated using the bulk aerodynamic method (Munro, 1990), which requires inputs of T, RH, air pressure and wind speed (Rye and others, 2010). We correct for stability following Dyer (1974). T at each model gridcell is calculated using the atmospheric lapse rate of 6°Ckm^-1 in summer (1 May- 30 September) and 8°Ckm^-1 in winter (1 October-30 April). These values are seasonal averages of Steffen and Box’s calculations of mean monthly temperature lapse rates between JAR1 and Summit GC-Net stations. We verified these lapse rates by calculating the average lapse rates for both summer and winter between JAR2 and Swiss Camp for 2005, which, to one decimal place, are the same as those calculated by Steffen and Box (2001). Air pressure was assumed to decline at a rate of 10kPakm^-1. RH and wind speed are assumed spatially invariant over the glacier surface. The surface roughness lengths for snow (0.00022 m) and ice (0.00066 m) for the Paakitsoq region are taken as the mean values from Arnold and Rees (2003). The model is found to be relatively insensitive to surface roughness values. Furthermore, through the calibration process, small errors introduced by prescribing fixed values for some parameters are compensated for by small adjustments to other parameter values during their optimization (cf. Rye and others, 2010). For each gridcell, the SEB model produces melt (mmw.e.) on an hourly time-step, provided that T_s = 0°C and Q_M>0. This meltwater is then used as input to the subsurface model.

3.2. Subsurface model

We use the subsurface model developed by Rye and others (2010), which built on previous work by Wright (2005), Bassford (2002) and Greuell and Konzelmann (1994). The model is run using a time-step of 15 min in order to maintain numerical stability (Rye and others, 2010).

For each model gridcell, the subsurface model calculates temperature, density and water content on a one-dimensional grid extending at least 25 m vertically from the surface into the ice sheet. We only consider the upper 25 m of snow/firn/ice as this is the depth at which annual temperature oscillations can no longer be detected (Greuell and Konzelmann, 1994). The vertical gridcells range in size from 5 cm near the surface (where temperature gradients are largest) to 200 cm at the grid base (Bassford, 2002). Meltwater generated by the SEB model is able to percolate downward through the grid. Energy (and therefore meltwater) is also added to the subsurface grid through the penetration of SW#, which is attenuated according to Beer’s law (Greuell and Konzelmann, 1994;Bassford, 2002;Bougamont and others, 2005;Wright, 2005;Rye and others, 2010). Refreezing occurs in cells where the temperature is <0°C and the density is less than the density of ice. The cell below receives any residual meltwater if either of these conditions is not met, or if there is excess meltwater after refreezing. Meltwater percolates until it reaches the impermeable snow/ice interface. At this layer, superimposed ice may be formed and is calculated using the approach of Wakahama and others (1976). However, if the rate at which meltwater reaches the snow/ice interface exceeds the rate of superimposed ice formation, then excess water will form runoff from the glacier (Rye and others, 2010). The hourly runoff for each model gridcell is stored as output and will be used as input for a surface meltwater routing model in forthcoming work.

Following Rye and others (2010), densification of snow and firn in the model is driven entirely by surface melting and refreezing. This is a reasonable assumption as the study area is in the wet zone of the GrIS where these two mechanisms dominate the densification process, unlike at higher elevations in the dry zone where settling and packing dominate (Cuffey and Paterson, 2010).

3.3. Surface accumulation model

Hourly precipitation is distributed over the ice-sheet surface using an elevation-dependent precipitation gradient (Rye and others, 2010). Estimates of the elevation-dependent precipitation gradient for the western area of the GrIS from sea level to 2000 m vary widely between 5% and 20% (100 m)^-1 (Ohmura and Reeh, 1991;Bales and others, 2009;Burgess and others, 2010). This suggests that the precipitation gradient also exhibits substantial year-to-year change and/or varies with distance inland (Bales and others, Owing to this uncertainly we treat the precipitation gradient as a tunable model parameter. Once set, however, it is constant in space and time throughout the model run. The fractions of precipitation falling as rain and snow in each gridcell are calculated as a function of air temperature using a threshold temperature (e.g. Oerlemans, 1991b;Bouga-mont and others, 2005;Rye and others, 2010). The value for the threshold temperature for solid/liquid precipitation is treated as a tunable model parameter.

3.4 Initial conditions

The separate model simulations for the mass-balance years 2000/01 and 2004/05 run from 1 September to 31 August. Following Rye and others (2010), each gridcell is initialized as bare ice and all subsurface cells are set to the mean annual air temperature. Before commencing the main model run for each year, the model is spun up for 5 years using that year’s climate data. While our sensitivity tests indicated that years of spin-up were necessary in order for the mass balance and the subsurface temperature profiles to attain a steady state, Rye and others (2010) found 20 years of spin-up were necessary when modelling mass balance for glaciers in Svalbard, and Bougamont and others (2005) found that just months of spin-up was sufficient for a similar model applied to the whole of the GrIS. This difference in required spin-up times between the three studies may be due to the varying initial conditions, which will not influence the steady state, only the time taken for the iterative process to reach a stable solution (Rye and others, 2010).

4. Model Calibration and Evaluation

Calibration involves optimizing model parameters so that modelled results are consistent with the available measurements, in our case continuous recordings of snow depth and albedo made at the GC-Net weather stations JAR1, JAR2 and Swiss Camp (Steffen and Box, 2001) for both 2000/01 and 2004/05. We evaluate the model in two ways. First we test the calibrated model output of snowline position against snowline positions derived from satellite imagery. Second, we compare modelled albedo with satellite-derived albedo over the model domain. Both of these steps are carried out for various stages of the summer for both 2000/01 and 2004/05 In order to model the temporally and spatially varying magnitude of meltwater runoff as accurately as possible, the two mass-balance years are parameterized separately.

We do not have suitable observations to constrain values for three key parameters, so we perform multiple model runs and compare the results with measurements in order to identify suitable values. The parameters are (1) fresh snow density, (2) elevation-dependent precipitation gradient, and (3) the threshold temperature for solid/liquid precipitation. Combinations of parameter values for model runs were chosen systematically from within defined ranges at specific intervals suggested in the literature (Oerlemans, 1991b; Bougamont and others, 2005;Bales and others, 2009; Burgess and others, 2010;Cuffey and Paterson, 2010;Rye and others, 2010) (Table 1). For fresh snow density, we parameterize for the range 100–450 kg m^-3. As snow densification due to settling, compression and the action of the wind is not accounted for by the model, but is instead driven by melting and refreezing alone, a value towards the upper end of this range might be expected during model calibration (Bassford, 2002; Wright, 2005; Rye and others, 2010). The value will represent average snowpack density before the onset of melting rather than the density of fresh new snow (Wright, 2005). Not all combinations of parameter values within these three ranges are tested. Instead, runs with combinations of high, medium and low values of the three model parameters were performed to identify a first tier of model sensitivity: a total of 3³ =27 runs. By analysis of the RMSE between sets of measured and modelled data, regions of parameter space showing the greatest sensitivity were then identified and parameters were tuned further at finer increments as specified in Table 1.

Table 1. Optimal parameter values and the ranges from which they were chosen

4.1. Model calibration using snow depth and albedo measurements

For each of the mass-balance years 2000/01 and 2004/05 we ran the surface mass-balance model many times with different combinations of parameter values for the ~7.5 km wide strip of the Paakitsoq region (Fig. 1; Table 1). For each year, daily variations in modelled surface height and surface albedo for the cells representing the GC-Net stations JAR1, JAR2 and Swiss Camp were compared, where data were available, with surface height measured at these stations by ultrasonic depth gauges (UDGs) and with albedo measured by up- and down-facing pyranometers. Modelled surface height at a point was calculated by multiplying the modelled melt (mm w.e.) for a surface cell by the ratio of water density to modelled density for that surface cell.

Surface height data are available at JAR1 and JAR2 for 2000/01 and 2004/05 but unfortunately are unavailable at Swiss Camp for either year. Albedo measurements are available at JAR1, JAR2 and Swiss Camp for 2000/01, and JAR1 and JAR2 for 2004/05.

The aim during model calibration was to identify the combination of parameter values that minimized the RMSE difference between the measured and modelled data. For each possible parameter set for each year, we calculate two RMSEs, one based on all the available measured surface height data, and one based on all the available albedo data, from the various GC-Net stations, for that year.

4.2. Model evaluation through comparison of modelled and measured snowline position

Snowline delineation from satellite data

We use Landsat-7 ETM+ imagery to delineate the ‘snowline’ position using a combination of the normalized-difference snow index (NDSI) and image thresholding (cf. Gupta and others, 2005). We use the word ‘snowline’ although we realise that in reality it is not a clear-cut boundary, but rather a transitional zone involving patches of snow, ice and slush.

The NDSI technique has been used previously to distinguish between highly reflective snow/ice and cloud (Hall and others, 1995;Konig and others, 2001). The method exploits the different spectral signatures of surfaces at different wavelengths. However, it is much harder to distinguish between snow and ice, and therefore the snowline position, using this method alone. Thus, once cloud was detected using the NDSI, thresholding of the imagery was used to delineate the snowline (see below). Combined, these methods provide an effective and quick way to delineate cloud from snow/ice, and snow from ice. It should be noted that all kinds of snow including slush, firn and fresh snow are included within the ‘snow’ category. This is consistent with the definition of the ‘snowline’ as being the boundary between the wet-snow and superimposed-ice zones, or the boundary between firn and glacier ice at the end of the melt season (Greuell and Knap, 2000;Cuffey and Paterson, 2010; Cogley and others, 2011).

Landsat-7 ETM+ data were acquired for days of our chosen mass-balance years with minimal cloud cover from the United States Geological Survey (http://glovis.usgs.gov). Before the NDSI calculation was made, data were preprocessed. First, an atmospheric correction was applied to raw, 8-bit, level 1b data using the dark pixel subtraction method (Lathrop and others, 1991). Second, each band was converted from digital numbers (DNs) to radiance and then from radiance to reflectance using the equations of Chander and others (2009).

Bands 2 (B2;corresponding to green, i.e. 0.519–601 mm) and 5 (B5;corresponding to the mid-infrared, i.e. 1.547–1.748mm) were then used to calculate the NDSI:

(8)

The NDSI was then used to distinguish between areas of snow/ice and cloud in the images. The technique works best for ‘optically thick’ clouds (e.g. cumulus clouds). The method exploits the high spectral reflectivity of clouds in the visible and shortwave infrared (SWIR) region compared to the relatively low reflectivity of ice and snow at SWIR wavelengths (Greuell and Knap, 2000; Griffin and others, 2005). In the image produced through the NDSI, the thickest cloud is the darkest in colour, and the thinner the cloud, the lighter the grey colour is. Snow is near-white in colour. If cloud is detected then we exercise caution in these areas when delineating the snowline using the subsequent thresholding technique (see below). However, it is not useful to immediately mask off all areas of the Landsat image overlain by cloud, as some cloudy areas are optically thin and give some information on the surface conditions below.

Landsat band 4 (corresponding to the near-infrared band, 1.772–1.898 mm) was selected to generate thresholded snowline images because it is less susceptible to saturation, thus making snow/ice distinction clearer. The brightness threshold is based on the analysis of Landsat band 4 histograms for each individual satellite image. For each image, we identified (a) five areas thought to contain only snow, (b) five areas thought to contain only ice and (c) five areas containing an obvious zone between ice and snow. Only cloud-free areas (identified using the NDSI technique) were selected. The size of these areas ranged from ~50 to ~100km². For illustrative purposes, the Landsat image in Figure 1 is overlain by three boxes containing (a) an area of snow only, (b) an area of ice only and (c) an area containing both ice and snow.

Individual histograms for each of the 15 areas of each image were produced. Again, for illustrative purposes, three histograms are shown in Figure 1, corresponding to the three boxed areas (a-c). For each image, the maximum brightness values in each of the five areas of ice were averaged. Similarly, the minimum brightness values in each of the five areas of snow were averaged. For the five areas containing a snow/ice boundary, the brightness values of the troughs in the bimodal histograms were averaged. The average upper limit for ice and the lower limit for snow were used to help better constrain the bimodal trough value. For example, the maximum brightness value for ice in Figure 1a is 147, and the minimum brightness value for snow in Figure 1b is 140. Therefore the ice/snow threshold value is expected to lie somewhere between 140 and 147. This is confirmed by the minimum brightness value in the trough between the ice and snow peaks in Figure 1c, determined as being 144. We would therefore take 144 as the threshold value for the ice/snow transition, and apply it to the Landsat band 4 brightness values for snowline delineation. Although the threshold value is clear in this example, for histograms with low, broad troughs, or for histograms of Landsat areas containing lakes or large areas of surface water, it is useful to have both the average maximum brightness value for ice and the average minimum brightness value for snow in order to help constrain the limits for the ice/snow threshold value.

Snowline delineation from model output

A key output from the calibrated model runs for each year was gridded cell albedo values for each hour. Cells with an albedo greater than that of ice (set at 0.48) were classified as ‘snow’. ‘Snow’ therefore includes a wide range of albedo values representing snow at different stages of metamorphism and melt, including fresh snow, firn and slush. This is consistent with how we categorize snow in the satellite imagery using the thresholding technique (above).

Comparison of measured and modelled snowline position

For each year (2000/01 and 2004/05) and for each day of the summer for which we have a thresholded satellite image, both the image and the model output consist of a grid containing cells classified as either ice or snow. The cloud cover produced from the NDSI indicated that there was no cloud within the strip on any day for which we were comparing modelled and measured snowline position. To evaluate model performance with respect to the snowline position, we calculate the percentage of cells that are mismatched (i.e. where modelled snow/ice cover does not match observed). For each year, the calibrated mass-balance model is run using the optimal set of parameter values selected in Section 4.1, and the average value of mismatched cells across the two or three available satellite images is calculated.

4.3. Model evaluation through comparison of modelled and satellite-derived albedo data

We compare modelled daily albedo with the daily surface albedo retrievals from the NASA Terra platform MODIS sensor MOD10A1 product (version 005), available from the US National Snow and Ice Data Center (NSIDC) (Hall and others, 2006). Surface albedo at 500 m resolution is calculated using the first seven visible and near-infrared MODIS bands (Klein and Stroeve, 2002;Klein and Barnett, 2003). Original MODIS data downloaded from NSIDC consisted of one tile (1200 km x 1200 km) gridded in a sinusoidal map projection and containing the entire area under study. Albedo values in the product are reported as per cent together with other flag values such as cloud obscuration and self-shadowing. The MODIS reprojection tool (MRT, https://lpdaac.usgs.gov/tools/modis_reprojection_tool) was used to re-project the data in Universal Transverse Mercator (UTM) coordinates and extract daily MODIS albedo data from late spring to late summer for the Paakitsoq region for both 2000/01 and 2004/05. Data tiles with >5% missing data (due to cloud coverage) are discarded. From the remaining data we select five relatively evenly spaced dates, from early in the melt season to late in the melt season, in each year. In order for the MODIS-derived albedo data (500 m) to be directly compared to the modelled albedo data (100 m), we resample the model output to 500m using bilinear interpolation. As ice albedo is set at the constant value in the model (0.48), we ignore all gridcells that have a modelled albedo equal to this. Thus, we focus primarily on snow- covered cells. Subsequently, for each date, the MODIS- derived albedo value for each gridcell is plotted against the modelled albedo value for each corresponding gridcell, and R² values and RMSEs for the relationships are calculated.

5. Results and Discussion

5.1. Model calibration

The lowest RMSEs between the modelled and measured surface height data for 2000/01 and 2004/05 were calculated to be 0.227 m and 0.208 m respectively (Fig. 2). The lowest RMSEs between the modelled and measured albedo data for 2000/01 and 2004/05 were calculated to be 0.084 and 0.118 respectively (Fig. 3). Optimal values for the key parameters of fresh snow density, elevation-dependent precipitation gradient, and the threshold temperature for solid/liquid precipitation that minimized the RMSEs were the same for both years and were 400 kg m^-3, an increase of 14% (100 m)^-1, and 2°C respectively (Table 1).

Fig. 2. Modelled and measured surface height data at GC-Net stations JAR2 (a, c) and JAR1 (b, d) for 2000/01 (a, b) and 2004/05 (c, d).

Fig. 3. Modelled and measured albedo at GC-Net stations JAR2 (a), JAR1 (b) and Swiss Camp (c) for 2000/01, and JAR2 (d) and JAR1 (e) for 2004/05.

Figure 4a and b show the relationships between the modelled and measured surface height and albedo data respectively, for all available measurements across both years. The R² values for the relationships are 0.99 and 0.70 respectively, indicating a very good overall match between modelled and measured data.

Fig. 4. Scatter plots of measured and modelled (a) surface height data and (b) albedo data. Also shown are the best-fit regression lines through each dataset.

A fresh snow density of 400 kg m^-3 is in agreement with the value used by Wright (2005) and Rye and others (2010). A linear elevation-dependent precipitation gradient of an increase of 14%(100m)^-1 is well within the range of possible values suggested by Bales and others (2009) and Burgess and others (2010). A value of 2°C for the threshold temperature for solid/liquid precipitation is consistent with that used by Bougamont and others (2005) and has been widely used in other studies (e.g. Oerlemans and Hoogen-doorn, 1989;Arnold and others, 1996).

The main discrepancies between the modelled and measured surface height data appear to be due to small snowfall events which are measured by the UDG but not captured by the model, as seen during the accumulation period (September-May) in Figure 2d, for example. Furthermore, the model slightly overestimates ablation during the summer, again likely due to the model not capturing small snowfall events during the ablation period. Similarly, the measured albedo values are generally more variable than the modelled values. This is particularly obvious, for example, in March/April 2005 when the effects of snowfall events appear in the measured albedo data at both JAR2 and JAR1 but not in the modelled albedo data (Fig 3d and e). Additionally, the drop from a high snow albedo to a lower ice albedo in late June 2001 at JAR1 may be less rapid in the measured data than in the modelled data due to small snowfall events temporally increasing the surface albedo and delaying the snowpack removal and subsequent exposure of underlying ice.

These differences between modelled and measured surface height and albedo may be mainly due to: (1) real differences in precipitation between the coastal station and the GC-Net sites due to local weather patterns;(2) a fairly high value for fresh snow density (400 kgm^-3), meaning a relatively small increase in surface height is modelled compared to that measured by the UDG;(3) the effect of the model’s relatively simple empirical albedo parameterization (e.g. compared to the more sophisticated, physically based model of Gardner and Sharp, 2010) which averages the density of the top 5 cm of snow/ice, masking the potential increase in albedo from a new snowfall; (4) a single value for the threshold temperature for solid/liquid precipitation not being applicable for all precipitation events; (5) the effect of wind-blown snow not being accounted for by the model;and (6) errors associated with the AWS measurements used to force the SEB model (e.g. Van den Broeke and others, 2004).

5.2 Model evaluation through comparison of modelled and measured snowline position

Following model calibration, evaluation was first undertaken by comparing the modelled snow distribution with that delineated from Landsat-7 ETM+ satellite imagery using the techniques of NDSI classification and supervised image thresholding, as outlined in Section 4.2. When the model was run with the optimal parameter set discussed in Section 5.1, on average the model successfully calculated the distribution of snow and ice for 90.4% and 87.9% of the cells in 2000/01 and 2004/05 respectively (Table 2). As an example, modelled and measured ice/snow distribution for (a) 7 July 2001 and (b) 8 August 2001 are displayed in Figure 5. As shown in Table 2, the majority of the total percentage of mismatched cells (96%) is due to the modelled presence of snow where the Landsat image indicates ice, i.e. the overestimation of snow cover.

Fig. 5. Correspondence between modelled and measured snow and ice distribution for (a) 7 July 2001 and (b) 8 August 2001, produced using the procedure described in Section 4.2.

Table 2. Percentages of gridcells in each of the four categories, and the total percentages of mismatched cells, for two dates in 2001 and three dates in 2005, following comparison of modelled snow and ice distributed with that delineated from Landsat imagery

There are two possible explanations for the small discrepancies between the modelled snowline position and the snow/ice pattern observed in the Landsat imagery. First, the error may lie primarily with the model, rather than with the thresholded imagery. As we have calibrated the model against data at three specific points (JAR1, JAR2 and Swiss Camp), it is not surprising that we cannot reproduce exactly the overall patterns of snow/ice distribution for the entire 450 km² region at different times in the summer. The snow identified through image thresholding is fairly patchy and dispersed, indicating that real snow accumulation patterns are more heterogeneous than modelled patterns, likely due to subtle complexities of the ice-sheet surface topography (i.e. patterns of curvature) which affect patterns of snowfall and redistribution by wind. These patterns are not picked up by the model, which distributes snow according to elevation only, and cannot account for redistribution by wind. For example, examining the locations of the measured snow- patches below the modelled snowline on both 7 July and 8 August 2001 (i.e. the yellow patches in Fig. 5) with reference to the DEM topography shows that the majority are located on the lee side of ridges in the DEM (given that the prevailing wind carrying snow is from the west-southwest (Fettweis, 2007)). Conversely, some of the areas of mismatch above the modelled snowline, where modelled snow is present but not measured (i.e. the light blue areas in Fig. 5), are on the windward side of ridges in the DEM.

The second explanation for the small discrepancies between the modelled snowline position and the snow/ice pattern observed in the Landsat imagery is that the error lies primarily with the thresholded imagery, rather than with the model. For example, we observe that many of the areas of mismatch above the modelled snowline (i.e. the light blue areas in Fig. 5) are linear in appearance and are aligned with many of the blue, wet-looking areas of snow visible in the Landsat imagery. The is particularly obvious when Figure 5a is compared to the Landsat image in Figure 1, which are both dated 7 July 2001. Thus, it is likely that our image- thresholding technique has incorrectly classified these wet, dark areas as ice rather than slushy snow. This problem likely arises because the wavelengths of Landsat band 4 are influenced by changes in grain size and water, meaning that it can be difficult to distinguish between wet snow and ice (Gardner and Sharp, 2010).

5.3. Model evaluation through comparison of modelled and satellite-derived albedo data

The second stage of model evaluation involved the comparison of gridded modelled albedo with the MODIS M0D10A1 albedo product, as outlined in Section 4.3. Ignoring all cells with a modelled albedo of 0.48 (i.e. ice) (an average RMSE of 0.081 is calculated between all cells with a modelled albedo of 0.48 and their corresponding MODIS-derived albedo value), the calculated R² values and RMSEs for the modelled vs MODIS-derived snow albedo data on five dates during both 2000/01 and 2004/05 are shown in Table 3. For example, Figure 6a-c show the relationships between modelled and MODIS-derived snow albedo for three dates in 2001 representing early (5 June) mid- (5 July) and late summer (11 August) respectively. Figure 7 shows the spatially varying difference between modelled and MODIS-derived snow albedo values over the model domain for the same three dates.

Fig. 6. Scatter plots of modelled and MODIS-derived albedo data for (a) 5 June 2001, (b) 5 July 2001 and (c) 11 August 2001. All gridcells with a modelled ice albedo of 0.48 have been removed. The 1: 1 lines and the best-fit regression lines are also plotted for each dataset.

Fig. 7. Spatially varying differences between MODIS-derived and modelled snow albedo values on (a) 5 June 2QQ1, (b) 5 July 2QQ1 and (c) 11 August 2QQ1. All gridcells with a modelled ice albedo of Q.48 have been removed. Positive (negative) values correspond to a higher (lower) MODIS-derived albedo than modelled albedo.

Table 3. Calculated R² values and RMSEs of the relationships between the modelled and MODIS-derived snow albedo data for 5 days in both 2001/01 and 2004/05

In early summer, the modelled albedo spans the complete range of albedo values, from new snow (0.82) to old snow or ice with a thin covering of snow (~0.50), and this is generally mirrored in the MODIS data (Fig. 6a). By midsummer, the modelled and measured albedo values span approximately the same range, but the MODIS albedo values display an increased heterogeneity for a given modelled albedo (Fig. 6b). By late summer, the majority of the snow has melted, leaving only modelled and MODIS albedo values at the lower end of the range (Fig. 6c).

In general, the model overestimates albedo values at the higher end of the range compared to MODIS values, and underestimates albedo values at the lower end of the spectrum compared to MODIS values (Fig. 6). The cells corresponding to these albedos occur in the upper (negative values) and lower (positive values) part of the model domain respectively (Fig. 7). However, the overestimation of modelled albedo values compared to MODIS values at the higher end of the range is less in early summer (Figs 6a and 7a) than in midsummer (Figs 6b and 7b), and the underestimation of modelled albedo values compared to MODIS values at the lower end of the range is greater in early summer (Figs 6a and 7a) than in midsummer (Figs 6b and 7b).

There are two possible reasons for this observation. First, it may be due to the model’s albedo parameterization scheme which calculates surface albedo as a linear function of the density of the top 5 cm of the subsurface grid. Thus, during early summer when melt rates are too limited to melt snow layers away and the subsurface snowpack is below freezing point, continuous rounds of melt (during the day) and refreezing (at night) will drive the density up and the albedo down. For higher melt rates and warmer subsurface snow- pack temperatures in midsummer, refreezing rates are lower, so densities are not increased as much, and albedos are not lowered as much as they are in early summer. Second, this observation may be due to the MODIS product identifying a range of levels of snow water saturation, resulting in an overall decrease in MODIS albedo values compared to values in the model, which does not account for the effect of surface water ponding. This second point is also the likely cause of the increased heterogeneity of MODIS albedos for a given modelled albedo as the summer progresses.

The discrepancies between modelled and MODIS albedo may also be due to unsubstantiated errors associated with the MODIS product (Stroeve and others, 2006) and/or errors associated with the AWS measurements used to force the SEB model (e.g. Van den Broeke and others, 2004).

5.4. Mass and energy budget

Figure 8 shows the calculated average seasonal cycle of the SEB components at the three GC-Net station sites, averaged over the two mass-balance years 2000/01 and 2004/05. As would be expected, the energy involved in melting, Q_M, i.e. the sum of all the incoming and outgoing fluxes, decreases with altitude from JAR2 to Swiss Camp. The net shortwave flux (SW_net), modulated at the surface primarily by albedo, is the dominant factor governing surface melt variability in the ablation area (Van den Broeke and others, 2008), and, on average, decreases towards higher elevations. The shape of the seasonal SW_net cycle at all three sites is not symmetrical. This is due to a gradual decrease in albedo associated with melt at the beginning of the melt season compared to a more rapid increase in albedo at the end of the melt season due to new snowfall. This effect is more noticeable higher up on the ice sheet at Swiss Camp (Fig. 8c) than lower down at JAR2 (Fig. 8a). The net longwave flux (LW_net) is the main energy loss and is strongly dependent on the temperature deficit at the surface, which increases with altitude. The SHF contributes significantly to the energy budget, whereas the LHF is small. Due to refreezing (which here includes both internal accumulation and superimposed ice formation), the GHF is positive. Averaged over the entire model domain, refreezing is calculated to contribute 0.13 m w.e. to the mass budget and net runoff averages –2.22 m w.e. Hence, on average, 6% of all meltwater and rainwater at the surface refreezes in the snowpack and does not become runoff. Averaged over the model domain, accumulation is 0.42 mw.e. Refreezing therefore accounts for 31% of the average net accumulation per year.

Fig. 8. Average seasonal cycle of the surface energy-balance components at (a) JAR2, (b) JAR1 and (c) Swiss Camp, averaged over the two mass-balance years 2000/01 and 2004/05. By definition, Q _M is negative (energy sink), but it is shown as a positive flux here for illustrative purposes.

6. Conclusions

We have applied a physically based surface mass-balance model to the Paakitsoq region of west central Greenland. The melt/runoff component is driven with a full range of meteorological variables collected on the ice sheet, predominantly from the GC-Net JAR1 station. The accumulation component uses precipitation data collected off the ice sheet near the coast from the ASIAQ station. The mass- balance model was calibrated by selecting sets of parameter values which minimized the RMSEs between modelled and measured surface height and albedo data from JAR1, JAR2 and Swiss Camp. We first evaluated the model performance by comparing the modelled snow distribution with that delineated from Landsat-7 ETM+ satellite imagery using the techniques of NDSI classification and supervised image thresholding for various dates during the summers of 2001 and 2005. Second, we evaluated the model through comparison of modelled daily albedo over the model domain with the daily surface albedo retrievals from the NASA Terra platform MODIS sensor MOD10A1 product. Our main findings are:

1. 1. Model calibration showed that the same set of optimal key parameters (fresh snow density (400 kg m^-3), elevation-dependent precipitation gradient (increase of 14% (100 m)^-1) and threshold temperature for solid/liquid precipitation (2°C)) were appropriate for both 2000/01 and 2004/05. This gives us confidence in the physical basis of the model, suggesting it is transferable between years. Calculated RMSEs between the modelled and measured surface height data for 2000/01 and 2004/05 were low at 0.227 m and 0.208 m respectively. The RMSEs between the modelled and measured albedo data for 2000/01 and 2004/05 were just 0.084 and 0.118 respectively. The calculated R² values between all of the available measured and modelled surface height data and all of the measured and modelled albedo data were 0.99 and 0.70 respectively.

2. 2. Model evaluation by comparison of modelled snowline position with that delineated from Landsat imagery shows that the average percentage of mismatched gridcells for 2000/01 and 2004/05 is relatively low at 9.6% and 12.1% respectively, with the majority of this mismatch (96%) due to the overestimation of snow cover. We also note that satellite-derived snow accumulation patterns are more heterogeneous than modelled patterns, likely due to subtle differences in real snow accumulation patterns due to small irregularities in the ice-sheet surface topography and the effect of snow redistribution by the wind. Some areas of error also exist in our thresholding scheme where areas of snow are wrongly classified as ice rather than wet snow.

3. 3. Although there is generally a good correspondence between gridded modelled albedo and MODIS-derived snow albedo, the model tends to overestimate albedo values at the higher end of the range compared to MODIS values, and underestimate albedo values at the lower end of the spectrum compared to MODIS values. This may be due to various reasons including: (1) the model’s albedo parameterization scheme which calculates surface albedo as a linear function of the density; the MODIS product identifying a range of levels of snow water saturation, not accounted for by the model; unsubstantiated errors associated with the MODIS product;and (4) errors associated with the AWS measurements used to force the SEB model.

4. 4. The average seasonal cycles of the SEB components are calculated and the net shortwave flux, modulated at the surface primarily by albedo, is the dominant factor governing surface melt variability in the ablation area. The net longwave flux is the main energy loss. We evaluate the spatial variability in annual net mass balance, runoff, accumulation and refreezing across the model domain and calculate that 6% of all meltwater and rainwater at the surface refreezes in the snowpack and does not become runoff; refreezing accounts for 31% of the average net accumulation.