Meltwater flux and runoff modeling in the ablation area of Jakobshavn Isbræ, West Greenland

Abstract

The temporal variability of surface snow and glacier melt flux and runoff are investigated for the ablation area of Jakobshavn Isbræ, West Greenland. High-resolution meteorological observations both on and outside the Greenland ice sheet were used as model input. SnowModel, a physically based spatially distributed meteorological and snow evolution modeling system, is used to simulate the temporal variability of Jakobshavn Isbræ accumulation and ablation processes for 2000/01–2006/07. Winter snow depth observations and MODIS satellite-derived summer melt observations are used for model validation of accumulation and ablation. The modeled interannual runoff variability varied from 1.81 × 10⁹ m³ (2001/02) to 5.21 × 10⁹ m³ (2004/05), yielding a cumulative runoff at the Jakobshavn Glacier terminus of ∼2.25 to ∼4.5 m w.e. The average modeled Jakobshavn runoff of ∼3.4 km³ a⁻¹ was merged with previous estimates of Jakobshavn ice discharge to quantify the freshwater flux to Illulissat Icefjord. 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 Jakobshavn freshwater flux to the fjord: about 7% (∼3.4 km³ a⁻¹) of the average annual freshwater flux of ∼51.0 km³ a⁻¹ originates from the surface runoff.

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. Box and others, 2006; Fettweis, 2007; Mernild 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 (IPCC, 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. Janssens and Huybrechts, 2000; Krabill and others, 2000, 2004; Zwally and others, 2002; Johannessen and others, 2005; Box and others, 2006; Fettweis, 2007; Mernild and others, 2008, 2009a; Hanna 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. Luckman and Murray, 2005; Chylek 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 (Parizek and Alley, 2004; Lemke and others, 2007; Van 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 (Nick and others, 2009). Fortunately, modeling the GrIS surface mass balance (SMB) is relatively well understood and documented in numerical models (e.g. Box and others, 2006; Fettweis, 2007; Mernild 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. Liston and Sturm, 1998; Pomeroy 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. Hanna and others, 2002, 2005, 2008; Mernild 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. Dowdeswell and others, 1997; ACIA, 2005; IPCC, 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 (Liston and Elder, 2006a; Mernild and others, 2006; Liston 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⁻¹; Rignot and others, 2008) into the fjord, draining approximately 6–7% of the GrIS by area (92 080 km²) (e.g. Luckman and Murray, 2005; Holland 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. Weidick 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 Steffen, 1995; Steffen 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 Born 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 Allerup and others, 1998, 2000).

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:

(1)

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.

SnowModel MODELLING SYSTEM: Description

SnowModel (Liston and Elder, 2006a) is a spatially distributed meteorological and snowpack evolution modeling system. It is made up of five submodels. MicroMet (a quasi-physically based meteorological distribution model) defines the meteorological forcing conditions (Liston and Elder, 2006b). EnBal calculates the surface energy exchanges, including melt (Liston, 1995; Liston and others, 1999). SnowPack simulates heat- and mass-transfer processes and snow-depth and water equivalent evolution (Liston and Hall, 1995). SnowTran-3D is a blowing-snow model that accounts for snow redistribution by wind (Liston and Sturm, 1998, 2002; Liston and others, 2007). Snow-Model also includes a snow data assimilation submodel (SnowAssim; Liston and Hiemstra, 2008) that can be used to assimilate available snow measurements to create simulated snow distributions that closely match observed distributions when and where they occur (Liston and others, 2008) (for a more detailed description of SnowModel, see, e.g., Liston and others, 2008; Mernild and others, 2008). SnowModel was originally developed for glacier-free landscapes. For glacier SMB studies in Arctic coastal regions, SnowModel was modified to simulate: (1) the glacier-ice melt after winter snow accumulation had ablated (Mernild and others, 2006, 2007); and (2) the influence of air-temperature inversions on snowmelt and glacier mass-balance simulations where radiosonde data are present (Mernild and Liston, 2010). For this study, routines for temperature inversion were not included due to the lack of available radiosonde data in the area. SnowModel has been used with highly reasonable results over a wide variety of snow and glacier landscapes in the United States, Norway, East Greenland, the GrIS and near-coastal Antarctica (Liston and others, 2008; Mernild and Liston, 2010 and references therein).

Input

To solve this system of equations, SnowModel requires spatially distributed fields of topography and land cover and temporally distributed point meteorological data (air temperature, relative humidity, wind speed, wind direction and precipitation) obtained from meteorological stations located within the simulation domain. For this study, high-resolution data are obtained from five meteorological stations: four stations (JAR1 to JAR3 and Swiss Camp) from the GC-Net, and one standard synoptic DMI WMO-operated station (Aassiaat) (Fig. 1; Table 1).

Table 1. Meteorological input data for the Jakobshavn SnowModel simulations. Meteorological station data on the GrIS (Swiss Camp, JAR1–JAR3) were provided by the GC-Net from CIRES and coastal meteorological station data (Aasiaat) by the DMI

Mean monthly lapse rates (1997–2005) based on air-temperature observations from the transect between JAR1, JAR2 and JAR3 were used as model input (Table 2). The minimum monthly temperature lapse rate of −10.0°C km⁻¹ occurred in October, and the maximum of −4.5°C km⁻¹ occurred in June. The mean annual Jakobshavn lapse rate of −7.1°C km⁻¹ is in line with the average western GrIS lapse rates of −7.8°C km⁻¹ (Steffen and Box, 2001) and the average GrIS values of −7.5°C km⁻¹ (Mernild and others, 2008).

Table 2. Mean monthly air-temperature lapse rates (°C km⁻¹) for the Jakobshavn area and the G-IS. For the Jakobshavn area, the mean monthly lapse rates are based on data from the transect between the meteorological stations from 1997 to 2005: JAR1 (962 m a.s.l.), JAR2 (542 m a.s.l.) and JAR3 (283 m a.s.l.) (see Fig. 1b for location of meteorological stations). For the GrIS, the mean monthly lapse rates are based on temperature data from the Greenland coastal areas and the GrIS (from 1997 to 2005) (for further information, see Mernild and others, 2008)

Across the Arctic, it is well known that precipitation gauges significantly underestimate solid precipitation because of aerodynamic errors at the gauging station, especially under windy and cold conditions (e.g. Yang and others, 1999). Solid and liquid precipitation measurements at the DMI meteorological station (Fig. 1; Table 2) were calculated from Helman–Nipher shield observations corrected according to Allerup and others (1998, 2000) and used as input for SnowModel. Solid (snow) precipitation for JAR1–3 was calculated from snow-depth sounder observations after the sounder data noise was removed and further used in SnowModel; these data are assumed to be accurate within ±10–15% (for further information about precipitation calculations, see Mernild and others, 2007, 2008).

The simulations span the 7 year period from 1 September 2000 to 31 August 2007. Simulations were performed using a 1 day time-step, although snow- and ice melt and blowing snow are threshold processes that may not be represented accurately by this time-step. Therefore, daily simulated melt and blowing-snow processes were compared against hourly simulated values from a test area, Mittivakkat Glacier (31 km²), southeast Greenland (Mernild and Liston, 2010), indicating significant values (p < 0.01; p is significance level and provided p < 0.10 the comparison is considered significant). The relative difference was 2%, 3% and 8% between daily and hourly simulated glacier winter, summer and net mass balance, respectively. We also recognize that daily averaged atmospheric forcing variables, in contrast to hourly data, smoothed the meteorological driving data.

Greenland topographic data for the model simulations were provided by Bamber and others (2001) and the image-derived correction by Scambos and Haran (2002). For the model simulations, this digital elevation model (DEM) was aggregated to a 500 m gridcell increment and clipped to yield a 343.5 km × 200.5 km simulation domain (the Ilulissat region) (Fig. 1). The domain includes the area of interest (127.0 km × 100.0 km) that encompassed the Jakobshavn Isbræ area (Fig. 1d). The Jakobshavn Isbræ ice front was confirmed or estimated for each year according to position illustrated on the satellite image (http://svs.gsfc.nasa.gov/vis/a000000/a003300/a003395/). SnowModel is a surface model producing first-order effects of climate change. It does not include glacio-hydrodynamic and glacio-sliding routines. Therefore, using a time-invariant DEM will yield melt uncertainties in areas where the GrIS surface elevation is modified by the arrival of relatively warm subsurface ocean water (Holland and others, 2008) and by GrIS dynamic processes. Between 1997 and 2001, for example, NASA’s Airborne Topographic Mapper surveys showed a ∼35 m reduction in surface elevation on the floating ice tongue (Holland and others, 2008). Observations from laser altimeter surveys along tracks on Jakobshavn Isbræ show dynamic thinning with rates of up to 10 m a⁻¹ between 1997 and 2003 (Krabill and others, 2004). The dynamic thinning continues, with increasing rates and with the thinning zone migrating inland as shown from the repeat survey between 2002 and 2006 (Joughin and others, 2008). In the simulations presented herein, dynamic thinning is assumed to be a second-order process and is not accounted for.

Each gridcell within the domains was assigned a US Geological Survey (USGS) Land Use/Land Cover System class according to the North American Land Cover Characteristics Database (e.g. Mernild and others, 2008). The snow-holding depth (the snow depth that must be exceeded before snow can be transported by wind) was held constant due to no changes in the bedrock and glacier/snow characteristics. The snow albedo is calculated using Douville and others (1995) and Strack and others (2004), gradually decreasing the albedo from 0.8 to a minimum of 0.5 as the snow ages. The albedo is reset to 0.8 after 0.003 mm SWE has fallen. When the snow is ablated, the GrIS surface ice conditions are used. Albedo was assumed to be 0.4 for ice; however, the GrIS ablation area is characterized by lower albedo on the margin and an increase in albedo toward the ELA, where a veneer of ice and snow dominates the surface (Bøggild and others, 2006). The emergence and melting of old ice in the ablation area creates surface layers of dust (black carbon particles) that were originally deposited with snowfall higher on the ice sheet. This debris cover is often augmented by locally derived wind-blown sediment. Particles on, or melting into, the ice change the area average albedo, increasing melt (these particles and debris cover were not accounted for in the simulations). User-defined constants for SnowModel are shown in Table 3 (for parameter definitions see Liston and Sturm, 1998, 2002).

Table 3. User-defined constants used in the SnowModel simulations (Liston and Sturm, 1998) for parameter definitions. If no interval is mentioned, the values are fixed throughout the simulation period

Calibration, validation and uncertainty

SnowModel was chosen for this study because of its strength (e.g. Mernild and others, 2009 and references therein) and ease of implementation over new simulation domains. This model demands limited input data, an important consideration in areas like Jakobshavn for which data are sparse due to rough terrain, harsh climatic conditions and remote location. To assess the general performance of SnowModelsimulated distributed meteorological data, snow evolution, snow and ice surface melt, glacier net mass balance and other snow and ice processes, simulated values were tested with independent observations. First, SnowModel/MicroMet distributed meteorological data were compared against independent Greenland meteorological station data both on and outside the GrIS, indicating respectable (84–87% variance for air temperature, 49–55% for wind speed, 49–69% for precipitation and 48–63% for relative humidity) representations of meteorological conditions (Mernild and others, 2008; Mernild and Liston, 2010). Second, only a few validation observations for in situ snow evolution, snow and ice surface melt and glacier net mass balance are available in Greenland. Therefore, SnowModel accumulation and ablation routines were tested qualitatively (by visual inspection) and quantitatively (cumulative values and linear regression) using independent in situ observations on snow-pit depths, glacier winter, summer and net mass balances, depletion curves, photographic time lapses, and satellite images also from in and outside the GrIS (e.g. Mernild and others, 2006, 2009a,b; Mernild and Liston, 2010). A comparison performed between simulated and observed values indicates a 7% maximum difference between modeled and observed snow depths, glacier mass balance and snow-cover extent.

To assess the winter and summer model performance for this Jakobshavn study, the simulated end-of-winter SWE depth and the summer melt extent were compared with Swiss Camp SWE depth and Moderate Resolution Imaging Spectroradiometer (MODIS) satellite-derived melt extent observations, respectively. Swiss Camp SWE depth was measured as an average of ten different measurements at the beginning of May and used to validate and adjust the SnowModel-simulated SWE depth. The simulated SWE was underestimated on average by 13% (155 mm w.e.) and up to 42% (430 mm w.e.) for 2003/04 (Table 4). The iterative precipitation adjustment routines (Mernild and others, 2006; SnowAssim: Liston and Hiemstra, 2008) yield a simulated Swiss Camp SWE depth on 10 May within 1% of the observed SWE depth. The MODIS satellite-derived melt extent (the discrimination between dry and wet snow found by reflected visible and near-infrared radiation) was observed for two days, 7 July 2002 and 20 July 2005, based on a spatial resolution of 1.0 km² (Chylek and others, 2007). The criterion for MODIS-derived snowmelt was a snow grain water thickness >40 μm. The upper part of the Jakobshavn drainage area is the area where the melting threshold of the algorithm did not show any melt. Satellite-derived and simulated average and standard deviated surface-melt discrepancy between melt and non-melt boundaries is 7.8 ± 5.1 km. However, the discrepancy can be up to ∼22 km (Fig. 2a). The difference in boundaries indicates a difference in melt area of ∼10% and in surface melt of ∼5% within the Jakobshavn Isbræ drainage area. It appears that our choice of SnowModel provided estimates of the Jakobshavn surface-melt distribution and related water-balance components agree well with observed values. Nevertheless, it is important to keep in mind the limitation for SnowModel results when: (1) tested against sparse observations; and (2) model uncertainties are determined largely by processes not yet represented by standard routines in the modeling system (e.g. routines for simulating changes in glacier area, size and surface elevation according to glacier dynamic and sliding processes). Further, all ocean and fjord areas within the domain were excluded from model simulations. Also, changes in glacier storage based on supraglacial storage, englacial storage, subglacial storage, meltwater routing and evolution of the runoff drainage system are not calculated in SnowModel, even though they certainly have some influence on runoff magnitudes.

Table 4. Observed and modeled SWE depth (mm w.e.) for Swiss Camp and the lower ablation zone of the Jakobshavn Isbræ drainage area at the beginning of May (10 May) and the end of winter (31 May) (for area specifications, see Fig. 1d)

Fig. 2. (a) MODIS satellite-derived melt extent for 7 July 2002 (DOY 188) and 20 July 2005 (DOY 202) including SnowModel-simulated line of melt extent for the Jakobshavn region. (b) Time series of daily modeled surface-snow and glacier-ice melt for the Jakobshavn Isbræ drainage area for 2001/02 (the year with the lowest annual cumulative surface melt) and 2004/05 (highest annual cumulative surface melt).

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 Mernild and others (2006) and used in mountainous areas of Norway (Young 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.

Table 5. Mean June–August temperature and positive degree-day (PDD) maximum modeled snowline elevation and runoff elevation, day of year (DOY) for first day of summer runoff, number of days with runoff, and runoff period, for 2000/01–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 Janssens and Huybrechts (2000) (used by, e.g., Hanna and others, 2002, 2005, 2008) and of 19–27% by Mernild and others (2008). The SnowPack submodel in SnowModel is similar to that used by Janssens 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 Hanna and others (2008) and Mernild 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 (Mernild 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 Rignot and others (2008).

Table 6. Modeled specific runoff for the lower ablation zone of the Jakobshavn Isbræ drainage area from 2000/01 to 2006/07. The runoff values do not include hydro-glaciological processes such as the sudden release of bulk water

For Jakobshavn, the ice discharge (calving) was estimated by Rignot 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 (Steffen 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.