Introduction

The Greenland ice sheet (GrIS) is the largest reservoir of permanent snow and ice in the Northern Hemisphere (∼7 m sea-level equivalent) and is highly sensitive to climate changes (e.g. Reference Bamber, Ekholm and KrabillBamber and others, 2001; Reference BoxBox and others, 2006). Satellite-based observations indicate that in 2007 and 2010 the GrIS had record surface melt extent (Reference MoteMote, 2007; Reference TedescoTedesco, 2007; Reference SteffenSteffen and others, 2008; Reference TedescoTedesco and others, 2011; J.E. Box and others, http://www.arctic.noaa.gov/reportcard/greenland.html) accompanied by unusually high freshwater runoff to the ocean (Reference Mernild and HasholtMernild and Hasholt, 2009). The extent and duration of surface melting are important for a number of reasons. An altered melting regime can produce substantial differences in surface albedo and energy and moisture balances, especially because wet snow absorbs up to three times more incident solar energy than dry snow (Reference SteffenSteffen, 1995). Further, an altered melt regime can influence the ice sheet’s surface mass balance (SMB) and runoff and its dynamic and subglacial sliding processes. Mechanisms that link climate, surface hydrology, internal drainage and ice dynamics are poorly understood, and numerical ice-sheet models do not simulate these changes realistically (Reference Nick, Vieli, Howat and JoughinNick and others, 2009). This was one of the main reasons why the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4; Reference SolomonSolomon and others, 2007) did not include future predictions for GrIS mass loss contribution related to global sea-level rise. Without a deeper understanding of the surface/ dynamic-sliding interactions and the subsequent hydrodynamic response of ice sheets, neither a best estimate nor an upper bound of future sea-level rise is reliable. Therefore, to quantify and predict GrIS mass loss, including runoff conditions, knowledge about its long-term variations in surface melt extent and trends, including the temporal and spatial distributions, melt intensity and duration, needs to be further developed. This surface meltwater runoff represents about half the annual mass loss from the GrIS (Reference Zwally and GiovinettoZwally and Giovinetto, 2001; Reference Rignot and KanagaratnamRignot and Kanagaratnam, 2006), with iceberg calving generating approximately the other half.

In this paper, we examine the duration, extent and magnitude of the GrIS surface melting from 1960 through 2010. Melting is simulated using SnowModel, a spatially distributed meteorological snow and ice model (Reference Liston and ElderListon and Elder, 2006a,Reference Liston and Elderb; Reference Mernild and ListonMernild and Liston, 2010) forced with observed meteorological data for the period when reliable forcing data are available. We compare results with observations during the satellite era (1979–2010) and we simulate GrIS melting conditions and trends from 1960 to 2010. Further, we examine whether the 2010 record melt season lies outside the natural melt variability since 1960.

Methodology and Data

SnowModel (Reference Liston and ElderListon and Elder, 2006a,Reference Liston and Elderb; Reference Mernild, Liston, Hasholt and KnudsenMernild and others, 2006) is a spatially distributed system for modeling meteorological conditions, snow evolution, snow and ice melting and runoff, given surface meteorological forcing. It simulates surface energy and moisture exchanges, including snow and glacier melt, multi-layer heat- and mass-transfer processes in snow (e.g. snowpack temperature and density evolution) and ice hydrological features such as meltwater routing from the surface through snow, firn and glacial ice. SnowModel routines have been described and tested by, for example, Reference Liston and ElderListon and Elder (2006a,Reference Liston and Elderb), Reference Liston, Haehnel, Sturm, Hiemstra, Berezovskaya and TablerListon and others (2007), Reference Liston and HiemstraListon and Hiemstra (2008, in press), Reference Mernild, Liston, Hiemstra and SteffenMernild and others (2008, Reference Mernild, Liston, Hiemstra, Steffen, Hanna and Christensen2009, Reference Mernild, Liston, Hiemstra and Christensen2010, Reference Mernild, Liston, Hiemstra, Christensen, Stendem and Hasholt2011) and Reference Mernild and ListonMernild and Liston (2010) and used successfully in a variety of snow and glacier environments, including the Arctic and Antarctic.

Table 1.

Meteorological input data for the GrIS melt extent simulations. Stations were operated and data were provided by the Danish Meteorological Institute (DMI)), University of Colorado at Boulder (CU), Geological Survey of Denmark and Greenland (GEUS), University of Copenhagen (UC) and University of Utrecht (UU). Parameters are Ta: air temperature; Rh: relative humidity; Ws: wind speed; Wd: wind direction; and P: precipitation. For station locations see Figure 1

Atmospheric forcing required by SnowModel is provided by MicroMet (Reference Liston and ElderListon and Elder, 2006b), which assimilates and interpolates time series of air temperature, relative humidity, wind speed and direction and precipitation from surface meteorological stations near or within the simulation domain. MicroMet uses known relationships between meteorological variables and the surrounding landscape (e.g. topography and surface characteristics) to distribute these variables in physically plausible and computationally efficient ways (e.g. Reference Liston, Haehnel, Sturm, Hiemstra, Berezovskaya and TablerListon and others, 2007; Liston and Hiemstra, 2011). Data are interpolated horizontally to a regular grid using a Barnes objective analysis scheme (Reference BarnesBarnes, 1964, Reference Barnes1973; Reference Koch, desJardins and KocinKoch and others, 1983) that applies a Gaussian distance-dependent weighting function. Interpolation weights are determined objectively as a function of data spacing and distribution. At each time-step, air temperature, solar radiation, albedo, outgoing longwave radiation, latent heat flux, sensible heat flux, snowmelt, sublimation, snowmelt runoff, snow depth and snow water equivalent are calculated and made accessible to SnowModel. The SnowModel simulated melt extent (for each daily time-step) was determined by summing the number of gridcells where simulated melt occurred.

Melt simulations were performed on daily time-steps based on: (1) a 5 km digital elevation model (Reference Bamber, Ekholm and KrabillBamber and others, 2001) and land cover from the US Geological Survey (USGS) North American Land Cover Characteristics Database, version 2.0 (http://edcdaac.usgs.gov/glcc/na_int.html) and (2) observed atmospheric forcing data from 56 meteorological stations located both in coastal areas and on the GrIS (Fig. 1; Table 1). The simulations covered the period September 1959 through December 2010. User-defined constants used in the simulations are listed by Reference Mernild, Liston, Hiemstra, Steffen, Hanna and ChristensenMernild and others (2009).

Fig. 1.

Greenland simulation domain with topography (500 m contour interval) and the location of the coastal and GrIS meteorological tower stations (red dots). Station specifications are provided in Table 1.

Using the increase in microwave emissivity as liquid water forms in previously dry snow (Reference Stiles and UlabyStiles and Ulaby, 1980), passive microwave satellite sensors were used to determine surface and near-surface melt occurrence for comparison with the SnowModel results. A microwave emission model was used to determine the brightness temperature (T _B) associated with 1% volumetric water content each year for locations across the GrIS using the dynamic threshold algorithm given by Reference Mote and AndersonMote and Anderson (1995). These modeled T _B values are used as threshold values to distinguish melt from non-melt during the summer of a given year. While Reference Mote and AndersonMote and Anderson (1995) used the 37 GHz frequency, this study used the 18 GHz scanning multichannel microwave radiometer (SMMR) and 19.35 GHz Special Sensor Microwave/Imager (SSM/I), following Reference MoteMote (2007).

The 18 GHz horizontally polarized channel from SMMR (January 1979 to June 1987; http://nsidc.org/data/nsidc-0007.html) and 19.35 GHz channel from SSM/I (July 1987 to September 2010; http://nsidc.org/data/nsidc-0001.html) T _B data were extracted from the US National Snow and Ice Data Center polar stereographic daily grids with a nominal 25 km gridcell increment. The emission model simulated the T _B values associated with melt for each gridcell each year. The melt area was determined by summing the total area covered by cells where the observed T _B exceeded the modeled threshold T _B. A comparison of SSM/I melt occurrence with infrared land surface temperatures from the Moderate Resolution Imaging Spectroradiometer (MODIS) in 2006 and 2007 showed that the SSM/I indicated 5–8% greater melt area; the two products agreed for 93% of non-cloud-contaminated cells.

Results and Discussion

Since 1960, the maximum modeled GrIS melt extent has been in 2010 and the minimum in 1972 (Fig. 2a). For 2010, the maximum melt extent was ∼9.5 × 10⁵ km², or 52% of the GrIS surface. This slightly exceeded the previous record of 51% set in 2007 and was significantly above the 51 year average of 34%. The modeled standard deviation in melt extent from 1960 to 2010 was 9% ((∼6.2 ± 1.7) × 10⁵ km²). For 1972, a melt extent of 17% (∼3.1 × 10⁵ km²) was simulated (Fig. 2a). In contrast to the simulated values, the satellite-derived area that underwent melting was more extensive in 2002 than in 2010 (Fig. 3a), indicating surface melt on 80% (∼13.2 × 10⁵ km²) of the ice sheet. The SSM/I satellite-derived maximum 2002 melt extent is different from the results based on cross-polarized gradient ratio (XPGR) algorithms (Reference Abdalati and SteffenAbdalati and Steffen, 1997; Reference SteffenSteffen and others, 2008 and updated), where 2010 was a record year, close to values from 2007 (personal communication from K. Steffen, 2010). A possible reason for the extremely high SSM/I satellite-observed 2002 melt extent might be the greater sensitivity of the dynamic threshold algorithm to melt onset compared with the XPGR algorithm (Reference Anderson, Mote, Abdalati and ColbeckAnderson and others, 1996). Furthermore, while the dynamic threshold algorithm used here showed a slightly greater total area underwent melt in 2002, it also showed that 2007 and 2010 had more frequent melt and longer melt seasons than 2002. For 2010, simulated surface melting occurred at elevations as high as 3000 m a.s.l., whereas for 1972 surface melting occurred only up to 2400 m a.s.l. Simulated melt extent was distributed largely as expected, with no melting at higher elevations in the interior of the GrIS (Fig. 2a), while marginal regions had surface melt 76–100% of the time during the summer (May to September), with melting most pronounced in southwest Greenland (Fig. 2b).

Fig. 2.

(a) Maximum simulated GrIS surface melt extent for 1972 (the year with minimum melt extent within the time series 1960–2010) and 2010 (the year with maximum melt extent). The satellite-derived boundary between melting and melt-free areas is shown for 2010 (bold black curve). (b) Simulated 2010 melt frequency in percentage of total melt days. (c) The difference between 2010 simulated melt duration and the 1960–2010 mean, in days.

Fig. 3.

(a) Time series of simulated annual maximum GrIS melt extent (1960–2010) and satellite-derived maximum melt extent (1979–2010). (b) Time series of seasonal mean (May through September) simulated (1960–2010) and satellite-derived (1979–2010) GrIS melt extent and standard deviation. (c) Mean Greenland summer air temperature anomaly with standard deviation (1960–2010). (d) Unsmoothed and smoothed Atlantic multi-decadal oscillation index (http://www.esrl.noaa.gov/psd/data/timeseries/AMO/). (e) Simulated and satellite-derived seasonal melt extent in 2007 and 2010.

In Figure 2a, the maximum simulated 2010 melt extent is compared with satellite-derived melt extent. The simulated boundary between melting and melt-free areas is highly consistent with the satellite-derived boundary, with an average discrepancy of 40 ± 35 km. The discrepancy was larger (up to ∼160 km) in the northeast and northwest and in the southern part of the interior. Differences between satellite-derived and simulated results were most pronounced where the distance between meteorological stations was largest. However, the general correspondence between simulations and observations supports the use of SnowModel for analyzing melting conditions before the satellite era (1960–79) when reliable atmospheric forcing data are available.

The difference between the 2010 melt duration and the 1960–2010 mean is illustrated in Figure 2c. The melt duration in 2010 was 41–60 days longer than the 1960–2010 average, with the largest differences occurring at the southwestern and western margins of the ice sheet at elevations up to ∼2000 m a.s.l. According to meteorological observations and US National Centers for Environmental Prediction/US National Center for Atmospheric Research (NCEP/NCAR) reanalysis data, this melting was consistent with anomalous warm airflow from the south during spring and summer, and record high winter air temperature leading to warm pre-melt conditions (Box and others, http://www.arctic.noaa.gov/reportcard/greenland.html), along with somewhat lower-than-average winter precipitation. Low winter snowfall led to earlier exposure of glacier ice and of the previous year’s summer snow surface; these surfaces have a lower albedo than fresh snow, promoting greater solar absorption and increased melting (Reference OkeOke, 1987; Reference Douville, Royer and MahfoufDouville and others, 1995). For example, at Nuuk (64.2° N along Greenland’s west coast) the 2009/10 winter, spring and summer were the warmest since record-keeping began in 1873. Similarly, Aasiaat (69.0° N along Greenland’s west coast) had its warmest year since record-keeping began in 1951 (Box and others, http://www.arctic.noaa.gov/reportcard/greenland.html), and Greenland’s winter precipitation in 2009/10 was ∼35% below the 1960–2010 average.

Melt durations in low-elevation areas along the southeast and east margins were up to 30 and 20 days longer than average, respectively. Along the northeastern portion of the ice-sheet margin, where winter temperatures were lower than average, the 2010 melt duration was up to 20 days less than the long-term average, probably due to regional variability in the atmospheric circulation and a late start to the melt season because of a greater cold content of the snowpack. As a whole, the 2010 GrIS melt duration was 14 ± 16 days longer than the 1960–2010 average, with a longer than average melt season on 89% of the GrIS. Because our model simulations were performed at relatively high resolution, the spatial heterogeneity in duration was greater than that of coarser-resolution satellite observations as described, for example, by Reference MoteMote (2007 and updated) and Box and others (http://www.arctic.noaa.gov/reportcard/green_land.html). Overall, however, the simulated melt duration was consistent with the pattern identified by Box and others (http://www.arctic.noaa.gov/reportcard/greenland.html, fig. GL3) for 1979–2010 (Fig. 2c).

For 1960–2010, time series of simulated maximum melt extent showed interannual variability superimposed on a period from 1960 to 1972 that can be approximated by a linear decreasing trend (p < 0.10, where p is level of significance) of 7% (∼1.2 × 10⁵ km²) of the total area of the GrIS (Fig. 3a). After 1972, the maximum melt extent increased significantly (p < 0.01) by 22% (∼3.7 × 10⁵ km²) of the GrIS area, consistent overall with the increase in GrIS maximum melt extent estimated by Reference Fettweis, Tedesco, van den Broeke and EttemaFettweis and others (2011). Likewise, the satellite-derived maximum melt extent (Fig. 3a) increased significantly (p < 0.01) by 27% (∼4.4 ± 10⁵ km²) from 1979 to 2010 and is in general agreement with the SnowModel results. The correlation between satellite-derived maximum melt extent and the simulated maximum melt extent is r ² = 0.56 (based on linear regression; where r ² is the explained variance), indicating a statistical significance at the 95% level (p < 0.05; i.e. there is <5% probability that such a correlation between the two time series can be produced by chance) and an rmse of 55%.

The maximum melt extents shown in Figures 2 and 3a are snapshots of surface melt. For studying seasonal climate impacts on SMB and runoff, it can be more appropriate to study the mean of integrated seasonal melt extent, as shown in Figure 3b. The simulated mean annual melt extent was 15 ± 5% (1960–2010), and the simulated mean melt extent of 17 ± 4% ((∼3.0 ± 0.7) × 10⁵ km²) (1979–2010) is significant (97.5% quantile) compared with the satellite-derived melt extent of 15 ± 5% ((∼2.6 ± 0. 8) × 10⁵ km²) and within one standard deviation of the observed values (Fig. 3b). For the mean modeled melt extent, a significant (p < 0.01) decreasing trend amounting to 6% of the GrIS area (∼1.1 × 10⁵ km²) occurred from 1960 to 1972. After 1972, the GrIS exhibited a reversal of this trend, with a significantly increasing trend in melt extent of 13% of the GrIS area (2.4 × 10⁵ km²) found in the simulations (1972–2010; p < 0.01) and a slightly smaller increase of 12% (∼2.0 × 10⁵ km²) derived from the satellite observations (1979–2010; p < 0.01), meaning that larger areas of dry snow were turned into percolation facies. The simulated melt extent since 1972 indicates that the melt extent in 2010 averaged more than twice that of the early 1970s. Many studies have demonstrated that the GrIS melt extent is greater than in the past (Reference MoteMote, 2007; Reference TedescoTedesco, 2007; Reference HannaHanna and others, 2008; Reference SteffenSteffen and others, 2008 ; Reference Mernild, Liston, Hiemstra, Steffen, Hanna and ChristensenMernild and others, 2009, Reference Mernild, Liston, Hiemstra and Christensen2010; Reference Fettweis, Tedesco, van den Broeke and EttemaFettweis and others, 2011; Reference TedescoTedesco and others, 2011; Box and others, http://www.arctic.noaa.gov/reportcard/greenland.html) and SnowModel results offer additional confirmation of these trends at a higher resolution. The increase in melt extent is associated with an observed rise in Greenland mean summer air temperature of 2.3°C for 1972–2010 (Fig. 3c) and with a rise in coastal Greenland summer temperatures of 1.7°C based on available long-term weather-station data for 1991–2006 (Reference HannaHanna and others, 2008). Furthermore, the overall trends in mean melt extent (Fig. 3b) seem to correlate with the smoothed trends of the Atlantic multi-decadal oscillation (AMO) index (Fig. 3d). From 1960 to the beginning of the 1970s the smoothed AMO index decreased on average and thereafter it increased towards 2010, analogous with the trends in melt extent for the GrIS. Reference Chylek, Folland, Lesins and DubeyChylek and others (2010) showed that the Arctic detrended temperatures were highly correlated with AMO, suggesting the Atlantic Ocean as a possible source of Arctic climate variability.

The impacts from major episodic volcanic eruptions (e.g. Agung, Indonesia (1963), El Chichón, Mexico (1982) and Pinatubo, Philippines (1991)) do not appear to have a significant impact on melt extent variations (Fig. 3b), even though, as stated by Reference Hanna, Huybrechts, Janssens, Cappelen, Steffen and StephensHanna and others (2005), global dust veils generated by volcanic activity might cool the polar regions and suppress ice-sheet melt.

According to simulations, the GrIS had the lowest mean melt extent of 7% (∼1.2 × 10⁵ km²) in 1972 and the highest mean melt extent in 2010; the satellite observations also affirm that 2010 was a record year with respect to mean melt extent (Fig. 3b). The simulated minimum mean melt extent was due to a 1972 mean summer (June–August) temperature ∼0.7 ± 0.3°C below the 1972–2000 baseline, and the maximum melt extent was due to a 2010 record high mean temperature ∼1.9 ± 0.6°C above the baseline (Fig. 3) in combination with the lower than average winter precipitation. For 2010, the simulated mean melt extent was 28% (∼5.2 × 10⁵ km²) of the total ice-sheet area and the satellite-observed melt extent was 23% (∼4.2 × 10⁵ km²).

The interannual variability in simulated GrIS mean melt extent occurred simultaneously with observed variations in Greenland mean summer air temperatures (at a level of r ² = 0.65, p < 0.05, rmse = 14%) (Fig. 3c). In Figure 3e, time series of seasonal simulated and observed melt extent are shown for 2007 and 2010 (the previous and present record melting years). The variability in seasonal simulated and observed melt extent was very similar. However, by the end of the season in mid-August, the simulations for both years relative to the observations are overestimated by an average of 11 ± 5% of the GrIS area. The average seasonal cycle for 1960–2010 indicated that melt extent in the extreme years 2007 and 2010 was about one to two standard deviations greater than average, and that 1972 (the year with minimum melt extent) was one to two standard deviations lower than average (Fig. 4a).

Fig. 4.

(a) Time series of maximum daily simulated GrIS melt extent from 1972, 2007, 2010 and average 1960–2010, including standard deviations. (b) Maximum daily simulated GrIS melt extent from January through December for 1960–2010.

The trend in seasonal melt extent for the simulation period is shown in Figure 4b. From 1960 to 1972 the melting period (with a >10% melt extent) decreased by an average of 3 days a⁻¹. After 1972 the period of >10% melt extent increased by an average of 2 days a−¹, culminating in an extended melting period of ∼70 days: 40 days in spring and 30 days in spring and autumn (Fig. 4b). Not only did the melting period increase for the GrIS, but also the number of days with a maximum melt extent greater than 30%. The number of days with a melt extent greater than 30% increased from ∼25 days in 1972 to ∼80 days in 2010, indicating an increasing melt index for the GrIS, where the melt index is defined as the melting area times the number of melting days (Reference TedescoTedesco, 2007).

Conclusion

Our findings show that GrIS melt extent reveals unambiguously that the melt period decreased from 1960 to 1972 (before satellite observations were available) and increased thereafter at a rate of 2 days a−¹, giving an extended melting period of 70 days, beginning 40 days earlier in spring and ending 30 days later in autumn. This increased melt period included record melt extent in 2010 and abnormally long melt duration, 41–60 days greater than average (1960–2010) at the southwestern and western margin of the ice sheet, supporting the hypothesis that 2010 was an unusually warm and dry year. Simulations of the GrIS melt extent, such as those presented here, will be crucial for understanding SMB, including runoff conditions, since surface meltwater runoff represents about half of the annual mass loss from the GrIS (Reference Zwally and GiovinettoZwally and Giovinetto, 2001; Reference Rignot and KanagaratnamRignot and Kanagaratnam, 2006), with iceberg calving generating most of the other half.