Introduction

The Greenland ice sheet (GrIS) is the largest terrestrial permanent ice- and snow-covered area in the Northern Hemisphere. The GrIS is a reservoir of water that is highly sensitive to changes in climate (e.g. Reference BoxBox and others, 2006; Reference FettweisFettweis, 2007; Reference Mernild, Liston, Hiemstra and SteffenMernild and others, 2008; http://www.arctic.noaa.gov/reportcard). It is essential to assess the impact of climate change on the GrIS, since the temperature rise at higher northern latitudes is strongly correlated with global warming and confirmed to have increased at almost twice the global average rate over the past 100 years (Reference SolomonIPCC, 2007). Since 1957, Arctic air temperature increases have averaged >2°C (http://giss.nasa.gov/). A response to altered climate has already been observed on the GrIS, manifested by an accelerating surface melt extent, mass loss and freshwater runoff and thinning along the periphery (e.g. Reference Janssens and HuybrechtsJanssens and Huybrechts, 2000; Reference KrabillKrabill and others, 2000, Reference Krabill2004; Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002; Reference Johannessen, Khvorostovsky, Miles and BobylevJohannessen and others, 2005; Reference BoxBox and others, 2006; Reference FettweisFettweis, 2007; Reference Mernild, Liston, Hiemstra and SteffenMernild and others, 2008, Reference Mernild, Liston, Hiemstra, Steffen, Hanna and Christensen2009a; Reference Hanna, Cappelen, Fettweis, Huybrechts, Luckman and RibergaardHanna and others, 2009). Also, at the local scale (e.g. in the Jakobshavn region, West Greenland), the same trend in melt extent and ice thinning are observed (e.g. Reference Luckman and MurrayLuckman and Murray, 2005; Reference Chylek, McCabe, Dubey and DozierChylek and others, 2007).

Efforts to model the GrIS mass balance, its dynamic processes, changes and contribution to the global eustatic sea-level rise still suffer from important uncertainties and limitations (Reference Parizek and AlleyParizek and Alley, 2004; Reference Lemke and SolomonLemke and others, 2007; Reference Van den Broeke, Smeets, Ettema, van der Veen, van de Wal and OerlemansVan den Broeke and others, 2008). The mechanisms that link climate and ice dynamics are poorly understood, and current numerical ice-sheet models do not simulate these changes realistically (Reference Nick, Vieli, Howat and JoughinNick and others, 2009). Fortunately, modeling the GrIS surface mass balance (SMB) is relatively well understood and documented in numerical models (e.g. Reference BoxBox and others, 2006; Reference FettweisFettweis, 2007; Reference Mernild, Liston, Hiemstra and SteffenMernild and others, 2008; http://www.arctic.noaa.gov/reportcard). To estimate the impact of seasonally changing processes on the GrIS surface hydrological cycle, different seasonal processes need to be understood and accounted for. Throughout the GrIS, much of the winter precipitation falls as a solid under windy conditions. As winter progresses, the solid precipitation accumulates on the ground and is frequently redistributed during blowing-snow events. A further consequence of this blowing snow is that significant portions (10–50%) of snow cover can be returned to the atmosphere by sublimation of wind-borne snow particles (e.g. Reference Liston and SturmListon and Sturm, 1998; Reference Pomeroy and EsseryPomeroy and Essery, 1999).

As spring and summer progress, the variation, duration and intensity of snow and glacier melt increases in response to the impact of weather and climate (e.g. insolation, temperature inversions, wind speed) and surface characteristics (e.g. albedo, roughness). The moisture in this system also changes phase (solid, liquid, vapor) throughout the year as part of various physical processes and in response to the available surface and snowpack energy fluxes. It is also important to take into account the role of snowpack meltwater retention. The overall GrIS runoff would be overestimated by approximately 20–30% if no model retention/refreezing routines were included in the model simulations (e.g. Reference Hanna, Huybrechts and MoteHanna and others, 2002, Reference Hanna, Huybrechts, Janssens, Cappelen, Steffen and Stephens2005, Reference Hanna2008; Reference Mernild, Liston, Hiemstra and SteffenMernild and others, 2008). All these seasonally changing processes impact directly the seasonal evolution (mass fluxes) of the GrIS surface hydrological cycle, including the influx of fresh water to the ocean and its subsequent role in controlling global eustatic sea level (e.g. Reference DowdeswellDowdeswell and others, 1997; ACIA, 2005; Reference SolomonIPCC, 2007).

This study attempts to improve our quantitative understanding of the Jakobshavn Isbræ drainage area melt distributions and its surface water balance components, particularly changes in freshwater runoff and net mass balance. The goal of this study was to apply a well-tested state-of-the-art modeling system, SnowModel (Reference Liston and ElderListon and Elder, 2006a; Reference Mernild, Liston, Hasholt and KnudsenMernild and others, 2006; Reference Liston, Haehnel, Sturm, Hiemstra, Berezovskaya and TablerListon and others, 2007), to the Jakobshavn region. SnowModel routines were compared with independent in situ field snow water equivalent (SWE) depth and satellite-derived melt extent observations. We performed model simulations for a 7 year period (2000/01 to 2006/07) with the following objectives: (1) to simulate winter processes related to snow accumulation, snow redistribution by wind and snow sublimation for the Jakobshavn Isbræ drainage area; (2) to simulate summer snowmelt and glacier ice melt for the drainage area; (3) to compare modeled outputs with available independent observational datasets; (4) to generate time series and area distributed runoff fluxes from the seasonal snowpack and the exposed glacier surface to be used as meltwater inputs to hydrological models; (5) to compare the trends in simulated runoff with previous estimates of Jakobshavn ice discharge and merge simulated runoff with ice discharge to quantify the freshwater flux to Ilulissat Icefjord; and (6) to model the net mass balance and the annual variability of the equilibrium-line altitude (ELA) location.

Study Area: Physical Settings, Meteorological Stations and Climate

Jakobshavn Isbræ is located on the west coast of Greenland (69°N, 49°W) ∼55 km east of the town Ilulissat (Jakobshavn) (Fig. 1). The Isbræ is Greenland’s largest outlet glacier and a prolific exporter of ice (∼50 Gt a⁻¹; Reference Rignot, Box, Burgess and HannaRignot and others, 2008) into the fjord, draining approximately 6–7% of the GrIS by area (92 080 km²) (e.g. Reference Luckman and MurrayLuckman and Murray, 2005; Reference Holland, Thomas, de Young, Ribergaard and LyberthHolland and others, 2008). Since the first observation in 1850, there has been an almost continuous recession of the Jakobshavn Isbræ of ∼40 km through the east–west orientated Ilulissat Icefjord. The ice front has receded at a steady rate of about 0.3 km a⁻¹ (1850–1964) (e.g. Reference Weidick and BennikeWeidick and Bennike, 2007), coming to rest ∼15–18 km downstream of the present location in 2001, and then receding again, far more rapidly at ∼3 km a⁻¹, to the present location (http://svs.gsfc.nasa.gov/vis/a000000/a003300/a003395/).

Fig. 1. (a) The Ilulissat region in West Greenland, with the simulation area and the area of interest, including the Jakobshavn Isbræ drainage area (5340 km²). (b) Simulation area including meteorological tower stations. (c) Area of interest, with topography (gray shades, 100 m contour interval) and longitudinal profile. (d) Land-cover characteristics in the area of interest including the four meteorological stations used for air-temperature lapse rates (Swiss Camp (1 140 m a.s.l.), JAR1 (962 m a.s.l.), JAR2 (542 m a.s.l.) and JAR3 (283 m a.s.l.)), the watershed divide and the ELA.

The simulated Ilulissat region (68 872 km²; Fig. 1) includes the area of interest (12 750 km²; Fig. 1a), covering the western part of the GrIS and the lower part (the ablation area) of the Jakobshavn Isbræ drainage area (5340 km²). The area of interest is characterized by elevations ranging up to 1840 m a.s.l. The land cover is dominated by ice in the upper parts of the terrain and ocean/fjord and bare bedrock/vegetation in the lower parts (Fig. 1c and d).

Five meteorological stations are located within the simulation domain, four within the area of interest (Fig. 1d). Station Aasiaat (Egedesminde) (68°42′ N, 52°45′ W; 88 m a.s.l.; a standard synoptic Danish Meteorological Institute (DMI) World Meteorological Organization (WMO) meteorological station) is located within the town of Aasiaat representative of the Disko Bay and fjord conditions. Stations JAR3 (69°23′ N, 50°18′ W; 283 m a.s.l.), JAR2 (69°25′ N, 50°03′ W; 542 m a.s.l.), JAR1 (69°29′ N, 49°41′ W; 962 m a.s.l.) and Swiss Camp (69°34′ N, 49°19′ W; 1140 m a.s.l.) are all part of the Greenland Climate Network (GC-Net) located on the ice sheet and representative of the GrIS conditions (for further information about the GC-Net stations, see Reference SteffenSteffen, 1995; Reference Steffen and BoxSteffen and Box, 2001). Swiss Camp is located close to the W-GrIS ELA (i.e. the elevation where the net mass balance is zero).

The Ilulissat region is considered to be Low Arctic according to Reference Born and BöcherBorn and Böcher (2001). The mean annual air temperature (2000–07) was −7.2°C. Mean annual relative humidity was 83% and mean annual wind speed was 6.1 ms⁻¹. The corrected mean total annual precipitation (TAP) was 810 mm w.e. a⁻¹ (corrected after Reference Allerup, Madsen, Vejen and KajanderAllerup and others, 1998, Reference Allerup, Madsen and Vejen2000).

Water Balance Components

Throughout the year, different surface processes (snow accumulation, snow redistribution, blowing-snow sublimation, surface evaporation and melting) on snow and glacier ice affect the surface glacier mass balance and the high-latitude water balance, including runoff. The yearly water balance equation for a glacier can be described by:

where P is the precipitation input from snow and rain (and possible condensation), E is evaporation, SU is sublimation (including blowing-snow sublimation), R is runoff and ΔS is change in storage (ΔS is also referred to as the SMB) from changes in glacier storage and snowpack storage. Glacier storage also includes changes in supraglacial storage (lakes, ponds, channels, etc.), englacial storage (ponds and the water table) and subglacial storage (cavities and lakes). Glacier storage components are not accounted for in this study. Here η is the water balance discrepancy (error). The error term should be 0 (or small) if the major components (P, E, SU, R and ΔS) have been determined accurately. Here a change in storage is calculated by the residual value.

Results and Discussion

Table 4 shows beginning-of-May observed and modeled SWE depth variations from Swiss Camp. Defined by our precipitation adjustment scheme for Swiss Camp, the areaof-interest average SWE depth on 10 May varied from 440 mm w.e. (2002/03, the year with the lowest SWE depth) to 2020 mm w.e. (2004/05, the year with greatest SWE depth), averaging 1340 ± 480 mm w.e. For the assumed end of winter (31 May, recognized as the end of the accumulation period) the SWE depth showed similar values, of 440 mm w.e. for 2002/03 and 2050 mm w.e. for 2004/05 (arrows in Fig. 3a), averaging 1370 ± 490 mmw.e. In Figure 3a, the assumed end of winter is marked by an arrow illustrating that this does not necessarily correspond to the maximum simulated average SWE depth (simulated end-of-winter), with a maximum difference of 7 days and a difference of average SWE depth <9 mm w.e. The average modeled SWE depth variation for the Jakobshavn drainage area is illustrated in Figure 3a for 2002/03 and 2004/05, indicating an increasing SWE depth throughout the accumulation period (September to May), a decreasing average SWE depth throughout the ablation period (June to August) and an end-of-year net accumulation in SWE depth (Fig. 3a). Here ablation includes phase-change processes like evaporation, sublimation and melting. Within SnowModel, Snow-Tran-3D simulates spatial snow deposition patterns in response to erosion and deposition and EnBal calculates the energy flux available for snowmelt. For the end of year, SWE depth varied from 430 mm w.e. (2002/03) to 1040 mm w.e. (2004/05). A surplus of SWE depth is located above the snowline (defined as the boundary between bare ice and snow cover). Our analysis of the spatial snowline distribution, in response to accumulation and ablation, shows an annual average maximum elevation between 928 ± 26 m a.s.l. (2001/02) and 1397 ± 11 m a.s.l. (2000/01) (Fig. 3d; Table 5). Our analysis of the spatial end-of-winter SWE depth distribution, in response to erosion and deposition (Fig. 3b and c), shows an increasing accumulation from the GrIS margin to the higher inland elevations. This yielded an average SWE precipitation orographic effect for the Jakobshavn ablation area of 71 ± 16 mm w.e. (100 m)⁻¹, while the annual orographic increase was 83 ± 14 mm w.e. (100 m)⁻¹. These values are in line with gradients found in previous East Greenland studies on Ammassalik Island (65° N) by Reference Mernild, Liston, Hasholt and KnudsenMernild and others (2006) and used in mountainous areas of Norway (Reference Young, Bolton, Killingtveit and YangYoung and others, 2006).

Fig. 3. (a) Variation in average modeled SWE depth for the Jakobshavn Isbræ drainage area for the year with the lowest (2002/03) and the highest (2004/05) average end-of-winter SWE depth. The numbers by the arrows indicate the average SWE depth on 31 May (end-of-winter period: accumulation period) and the other numbers the average SWE depth on 31 August (end-of-summer: ablation period). (b) Spatial simulated SWE distribution for the end of winter (31 May 2003). (c) Spatial simulated SWE distribution for the end of winter (31 May 2005). (d) Annual modeled maximum elevated snowline (the boundary between bare ice and snow cover on the glacier surface) from 2000/01 to 2006/07.

Daily modeled time series of surface-snow and glacier-ice melt for 2001/02 (the year with the lowest cumulative melt) and 2004/05 (the year with the highest cumulative melt) are illustrated in Figure 2b. For 2001/02 and 2004/05 the cumulative melt reached 3.83 × 10⁹ m³ and 8.64 × 10⁹ m³, respectively, with maximum daily snowmelt values of 2.26 × 10⁸ m³ and glacier ice melt values of 1.39 × 10⁸ m³. In the early melt period (May and June), surface melt was mainly controlled by snowmelt, whereas later in the season (mid-/late July and August) when the snow cover was largely gone, surface melt was dominated by glacier ice melt. When surface melting is defined by SnowModel, meltwater is assumed to run instantaneously when the surface consists of glacier ice. When snow cover is present, the SnowPack runoff routines take retention and internal refreezing into account when water melts at the surface and penetrates the snowpack. These routines have a significant effect on the runoff lag time and the amount of runoff. Not including retention/refreezing routines in SnowModel would lead to faster outflow of runoff and overestimation of runoff to the ocean and consequent overestimation of the global sea-level rise. If no retention/refreezing routines were included in SnowModel: (1) the initial seasonal runoff would occur 23–85 days earlier; and (2) the Jakobshavn runoff would be overestimated by 65–110%, averaging 80%. This 80% value exceeds previous values for the entire GrIS of approximately 25% estimated by the single-layer snowpack model of Reference Janssens and HuybrechtsJanssens and Huybrechts (2000) (used by, e.g., Reference Hanna, Huybrechts and MoteHanna and others, 2002, Reference Hanna, Huybrechts, Janssens, Cappelen, Steffen and Stephens2005, Reference Hanna2008) and of 19–27% by Reference Mernild, Liston, Hiemstra and SteffenMernild and others (2008). The SnowPack submodel in SnowModel is similar to that used by Reference Janssens and HuybrechtsJanssens and Huybrechts (2000). It does not calculate vertical temperature changes through the snow-pack. For the GrIS in total, the SnowModel retention and refreezing routines indicate that high-runoff years are synchronous with low-precipitation/accumulation years and vice versa. This trend was reported for the GrIS by Reference HannaHanna and others (2008) and Reference Mernild, Liston, Hiemstra and SteffenMernild and others (2008) because higher volumes of meltwater were retained in the thicker snowpack, reducing runoff. It is most pronounced above the ELA, where meltwater does not infiltrate far into the snowpack because of the snowpack’s cold content even during summer. For the Jakobshavn ablation-area study the trend is opposite (R ² = 0.24, where R ² is the explained variance, p < 0.25).

Figure 4a illustrates the time series for the daily surface runoff production from both snow cover and glacier ice throughout winter and summer from 2001/02 to 2006/07, and Figure 4b presents the spatial distribution of cumulative runoff for 2001/02 (lowest annual cumulative runoff) and 2004/05 (highest cumulative runoff). The daily runoff values averaged 0.32 × 10⁸ m³, with a maximum daily value as high as 2.83 × 10⁸ m³, or 8.8 times the average runoff (Fig. 4a). During winter (September/October to May/June), no runoff events were simulated. For 2000/01 to 2006/07 exponential regression indicates an R ² value (the explained variance) of 0.70 (p < 0.01) between modeled daily runoff and mean daily air temperature; high temperatures correspond to high simulated runoff values. The first day for annual runoff varied, by almost 1 month, from the beginning of May (day of year (DOY) 129) to the beginning of June (DOY 159), averaging near the end of May (DOY 142; Table 5). Then a continuously modeled runoff period occurred until September/October, indicating the average number of runoff days to be 97 ± 17. In some areas of Jakobshavn Isbræ (e.g. at the glacier terminus), as much as ∼4.5 mw.e. of runoff was simulated for 2004/05, while only ∼2.25 m w.e. occurred in 2001/02 (Figs 4b, and 5a and b). The amount of simulated runoff decreased with increasing altitude, on average by 218 mm w.e. (100 m)⁻¹, from the ice margin all the way to the runoff discrepancy between the runoff and non-runoff boundaries. The annual runoff boundary was located ∼65 km from Jakobshavn Isbræ and further inland, at 1420–1870 m a.s.l., averaging 1600 ± 170 m a.s.l. (Fig. 5a; Table 5). SnowModel outputs were further used to calculate the annual cumulative runoff. The Jakobshavn runoff ranged from 1.81 × 10⁹ m³ a⁻¹ (2001/02) to 5.21 × 10⁹ m³ a⁻¹ (2004/05), indicating an average annual increase of 0.13 × 10⁹ m³ a⁻¹ (insignificant: R ² = 0.13, p < 0.25) for the simulation period (Figs 4 and 5d). The mean annual Jakobshavn runoff was 3.4 ± 1.1 × 10⁹ m³ a⁻¹ (= 3.4 km³ w.e.a⁻¹). This is equivalent to a specific runoff of 26.1 ± 5.5 L s⁻¹ km⁻² (Table 6), of the same order of magnitude as the value from the Kangerlussuaq drainage basin (67° N, 50° W) located 250 km to the south (for Kangerlussuaq the simulated runoff was validated against runoff observations), and almost four times the average specific runoff of 6.7 ± 1.0 L s⁻¹ km⁻² for the entire GrIS over the period 1995–2005 (Reference Mernild, Liston, Hiemstra and SteffenMernild and others, 2008). This indicates that the Jakobshavn runoff exceeded the spatial specific average runoff for the GrIS.

Fig. 4. (a) Time series of daily modeled runoff for the Jakobshavn Isbræ drainage area from 2000/01 to 2006/07. (b) Spatial simulated runoff distribution for 2001/02 (year with lowest annual cumulative runoff) and 2004/05 (highest annual cumulative runoff).

Fig. 5. (a, b) Cumulative modeled longitudinal runoff profile for 2001/02 (year with lowest annual runoff) (a) and 2004/05 (highest annual runoff) (b) calculated for every second week starting 1 June to 31 August (longitudinal profile in Fig. 1c). (c) Exponential relation between daily runoff and mean daily air temperature. (d) Time series for simulated annual runoff from 2000/01 to 2006/07 and Jakobshavn ice discharge from 2000 and 2004–07 based on data from Reference Rignot, Box, Burgess and HannaRignot and others (2008).

For Jakobshavn, the ice discharge (calving) was estimated by Reference Rignot, Box, Burgess and HannaRignot and others (2008) to be 52.6 ± 7.4 Gt a⁻¹ (= 47.3 ± 6.7 km³ w.e. a⁻¹) (for the years 2000 and 2004–07). For both runoff and ice discharge, the average trends are similar, indicating increasing (insignificant) influx of fresh water to Ilulissat Icefjord for the period 2000/01–2006/07 (Fig. 5d). Based on the few common data points (2004–07; n = 4; Fig. 5d), there is no reason statistically to conclude either that (1) the increasing runoff has any influence on the rapid increasing ice discharge or (2) an increasing flux of surface runoff in a warmer future climate will accelerate the volume of ice discharge. The mechanisms that link changing climate to changing surface conditions, glacier hydrology and ice-sheet dynamics are still poorly understood.

We combined ice discharge with the average Snow Model-simulated surface runoff to deduce the freshwater flux from the Jakobshavn drainage area (losses from geothermal heating/melting were omitted), and found a freshwater flux averaging about 50.7 km³ a⁻¹ to Ilulissat Icefjord. About 3.4 km³ a⁻¹ (∼7%) originated from the surface runoff and 47.3 km³ a⁻¹ (∼93%) from ice discharge, totaling ∼51.0 km³ a⁻¹.For the Jakobshavn drainage area, runoff forms a minor part of the overall freshwater flux to the fjord, whereas further south (e.g. at the Kangerlussuaq drainage area) runoff accounts for 100% of the flux, due to the inland position of the GrIS margin.

To assist with calculating the net mass-balance conditions and the location of the ELA, the water balance (Equation (1)) was divided into two different periods: (1) an accumulation period (September to May; winter period) where accumulation processes (precipitation and snow redistribution, influenced by blowing-snow sublimation) are dominant; and (2) an ablation period (June to August; summer period) where ablation processes (evaporation, sublimation and melting) are dominant. Figure 6 illustrates the simulated net mass-balance variations for the Jakobshavn drainage area. The modeled ELA fluctuated from 990 to 1210 m a.s.l., so that Swiss Camp was positioned within the boundaries of the ELA zone. JAR1 was positioned in the ablation area with a varying net mass balance of −0.8 to −0.2 m w.e. a⁻¹. JAR2 and JAR3 net mass balance varied between −2.4 and −1.1 m w.e. a⁻¹ and between −3.5 and −1.7 m w.e. a⁻¹, respectively. The location of the ELA is sensitive to changes in climate. The modeled ELA is in accordance with observations; however, Stober and Hep-perle (2007) indicated that the ELA had moved to an elevation ∼250 m higher than Swiss Camp by 2006. In general, in the 1990s Swiss Camp was located at the ELA, but the ELA moved to higher elevations at the beginning of the 21st century due to increased melt in the area (Reference Steffen, Zwally, Rial, Behar and HuffSteffen and others, 2006).

Fig. 6. Modeled net mass balance in relation to elevation for the period 2000/01 to 2006/07 for the Jakobshavn Isbræ drainage area.

Summary and Conclusions

Quantifying freshwater runoff where observed datasets are available is relatively easily achieved, albeit possibly subject to significant uncertainties. This study presents simulations of the Jakobshavn surface-melt extent and related water-balance components, focusing on surface runoff for the period 2000/01–2006/07. A robust physically based state-of-the-art snow and glacier ice evolution modeling system (SnowModel) was used. Our SnowModel simulations have been validated against independent in situ observations (accumulation and ablation observations) made on the western GrIS. There is a high degree of agreement between these Jakobshavn simulations and the recorded observations. The simulated Jakobshavn series yielded useful insights into the present conditions of the ice-sheet net mass balance and the interannual runoff variability. The 2000/01–2006/07 mean annual surface runoff was ∼3.4 km³ a⁻¹. In light of missing glacio-hydrodynamic model routines, values from previous studies of the Jakobshavn ice discharge, 52.6 ± 7.4 Gt a⁻¹ (= 47.3 ± 6.7 km³ w.e. a⁻¹), were adapted to provide estimates of the overall freshwater flux to Ilulissat Icefjord (∼51.0 km³ a⁻¹). For both runoff and ice discharge the average trends are similar, indicating increasing (insignificant) influx of fresh water to Ilulissat Icefjord for the period 2000/01–2006/07. This study suggests that surface runoff forms a minor part of the overall freshwater flux to the fjord: ∼7% (∼3.4 km³ a⁻¹) of the average annual freshwater flux of ∼51.0 km³ a⁻¹ originates from surface runoff.