Fluctuations in the Greenland ice sheet (GrIS) surface mass balance (SMB) and freshwater influx to the surrounding oceans closely follow climate fluctuations and are of considerable importance to the global eustatic sea level rise. A state-of-the-art snow-evolution modeling system (SnowModel) was used to simulate variations in the GrIS melt extent, surface water balance components, changes in SMB, and freshwater influx to the ocean. The simulations are based on the Intergovernmental Panel on Climate Change scenario A1B modeled by the HIRHAM4 regional climate model (RCM) using boundary conditions from the ECHAM5 atmosphere–ocean general circulation model (AOGCM) from 1950 through 2080. In situ meteorological station [Greenland Climate Network (GC-Net) and World Meteorological Organization (WMO) Danish Meteorological Institute (DMI)] observations from inside and outside the GrIS were used to validate and correct RCM output data before they were used as input for SnowModel. Satellite observations and independent SMB studies were used to validate the SnowModel output and confirm the model’s robustness. The authors simulated an ∼90% increase in end-of-summer surface melt extent (0.483 × 106 km2) from 1950 to 2080 and a melt index (above 2000-m elevation) increase of 138% (1.96 × 106 km2 × days). The greatest difference in melt extent occurred in the southern part of the GrIS, and the greatest changes in the number of melt days were seen in the eastern part of the GrIS (∼50%–70%) and were lowest in the west (∼20%–30%). The rate of SMB loss, largely tied to changes in ablation processes, leads to an enhanced average loss of 331 km3 from 1950 to 2080 and an average SMB level of −99 km3 for the period 2070–80. GrIS surface freshwater runoff yielded a eustatic rise in sea level from 0.8 ± 0.1 (1950–59) to 1.9 ± 0.1 mm (2070–80) sea level equivalent (SLE) yr−1. The accumulated GrIS freshwater runoff contribution from surface melting equaled 160-mm SLE from 1950 through 2080.
The Greenland Ice Sheet (GrIS) is the Northern Hemisphere’s largest terrestrial permanent ice- and snow-covered area and a reservoir of water, from a hydrological perspective (e.g., Box et al. 2006; Fettweis 2007; Richter-Menge et al. 2007; Mernild et al. 2008d, 2009a,b), containing between 7.0-m and 7.4-m global sea level equivalent (SLE) (Warrick and Oerlemans 1990; Gregory et al. 2004; Lemke et al. 2007). It is essential to predict and assess the impact of future climate on the GrIS, which is believed to be influenced by human activities (Albritton et al. 2001). We must establish the present and future state of the GrIS surface melt extent and surface mass balance (SMB), including freshwater flux, to detect warning signs indicative of its future response (Hanna et al. 2008). Variability in mass balance closely follows climate fluctuations; the mass balance was close to equilibrium during the relatively cold 1970s and 1980s and lost mass rapidly as the climate warmed in the 1990s and 2000s with no indication of deceleration (Rignot et al. 2008). The GrIS is a useful indicator to ongoing climatic variations and changes, and it is suggested that the GrIS responds more quickly to climate perturbations than previously thought, particularly near the margin in southern Greenland (Velicogna and Wahr 2006).
The climate appears to be changing. Observations indicate that the most pronounced temperature increase occurs at higher northern latitudes, which have increased at almost twice the global average rate in the past 100 years (Alley et al. 2007). Since 1957 air temperature for the Arctic has increased more than 2°C (available online at http://www.giss.nasa.gov/). The warming was accompanied by an increase in precipitation of ∼1% decade−1 (ACIA 2005). Simulations by atmosphere–ocean models for areas north of 60°N project an increased mean surface air temperature of 2.5°C by the mid-twenty-first century and 4.5°–5.0°C by the end of the twenty-first century (ACIA 2005; Alley et al. 2007).
A response to an altered climate has already been observed on the GrIS, manifested by thinning along its periphery (primarily in the south) and a slight thickening of ∼2–5 cm yr−1 in the interior (e.g., Krabill et al. 1999, 2000, 2004; Johannessen et al. 2005; Thomas et al. 2006; Zwally et al. 2005). Changes in air temperature result in large changes in the surface melt extent. The GrIS passive microwave satellite-derived surface melt extent increased 17.6 × 103 (1973–2007; Mote 2007) and 40.0 × 103 km2 yr−1 (1992–2005, Tedesco 2007); for 2007 a record GrIS melt extent occurred (e.g., Mernild et al. 2009a,b). Further, for the area above 2000-m elevation the 2007 melt index, defined as the melting area times the number of melting days, was 153% greater than the average for the period 1988–2006, setting a new record (Mote 2007; Tedesco 2007; Mernild et al. 2009b). In contrast to 2007, snowmelt over the whole GrIS in 2008 was not significant at high elevations. Melt extent in 2008 was, however, above the 1979–2007 average, with the 2008 updated melt extent trend of approximately 16 × 103 km2 yr−1 (Tedesco et al. 2008).
Numerous GrIS mass-balance studies using airborne laser altimetry and models (the positive-degree approach and energy balance) suggest a balance ranging between 25 and −60 km3 water equivalent (w.eq.) yr−1 (1961–2003), from −50 to −100 km3 w.eq. yr−1 (1993–2003) and a reduction at even higher rates between 2003 and 2005 to a loss of ∼270 km3 w.eq. yr−1 in 2007 (Lemke et al. 2007; Rignot et al. 2008; Mernild et al. 2009b). Analyses of the Gravity Recovery and Climate Experiment (GRACE) satellite data show mass loss of 75–129 km3 w.eq. yr−1 (2003–05) and losses ranging from 150 to 270 km3 w.eq. yr−1 (2002–07) (Velicogna and Wahr 2005, 2006; Chen et al. 2006; Lutchke et al. 2006; Ramillien et al. 2006). This indicates an accelerating GrIS mass loss in the 1990s up through the beginning of the twenty-first century, equivalent to a net global eustatic sea level rise of ∼0.5 mm SLE yr−1 (Velicogna and Wahr 2006).
Nearly half of the mass lost from the GrIS originates by surface melting and subsequent freshwater runoff into the ocean. The other half is from iceberg calving and geothermal melting (e.g., Lemke et al. 2007; Mernild et al. 2008c). Calculated runoff losses are provided by Janssens and Huybrechts (2000), 281 km3 yr−1 (1953–2003); Mote (2003), 278 km3 yr−1 (1988–99); Box et al. (2006), 396 km3 yr−1 (1995–2004); Hanna et al. (2008), 351 km3 yr−1 (1995–2007); Fettweis (2007), 304 km3 yr−1 (1979–2006); Mernild et al. (2008d), 392 km3 yr−1 (1995–2005); and Mernild et al. (2009a,b), 397 km3 yr−1 (1995–2007). Increases indicate an accelerating GrIS runoff, probably playing a potential role in ocean salinity, sea ice dynamics, the global eustatic sea level rise (e.g., Dowdeswell et al. 1997; ACIA 2005; Box et al. 2006; Alley et al. 2007), and thermohaline circulation (THC) of the Greenland–Iceland–Norwegian Seas (e.g., Broecker et al. 1985; Broecker and Denton 1990; Su et al. 2006). Accelerating GrIS runoff could perturb the THC by reducing the density contrast driving the circulation (Rahmstorf et al. 2005). Any weakening of the circulation in response to increased GrIS runoff induced by global warming (Gregory et al. 2005; Swingedouw et al. 2006) could reduce heat inflow to the Greenland–Iceland–Norwegian Seas and subsequently reduce the warming in the region, including northwest Europe.
This study attempts to improve our quantitative understanding of the past, present, and future (131-yr perspective, 1950–2080) GrIS surface melt extent and its related water balance components. Specifically, we address changes in the SMB and the influx of freshwater to the ocean as a contribution to the global eustatic sea level rise. The goal of this study is to apply a well-tested approach—a state-of-the-art snow-evolution modeling system, SnowModel (e.g., Liston and Elder 2006a,b; Liston et al. 2007; Mernild et al. 2006b, 2008c)—to the GrIS for the period from 1950 through 2080 based on the Intergovernmental Panel on Climate Change (IPCC) A1B climate scenario. The climate scenario is used in a high-resolution regional climate model (RCM)HIRHAM4 developed by the Danish Meteorological Institute (DMI) (Christensen et al. 1996; Stendel et al. 2008). The RCM output data was calibrated and tested using in situ meteorological observations obtained from the GrIS [Greenland Climate Network (GC-Net); 1995–2005] and the coast [World Meteorological Organization (WMO) DMI meteorological stations] before being used as meteorological forcings for SnowModel. SnowModel was tested by coincident passive microwave satellite images and SMB studies. We performed the GrIS model simulations for the 131-yr period (1950–2080) with the following objectives: 1) assess the HIRHAM4 RCM meteorological driving data against in situ meteorological observations; 2) quantify the GrIS end-of-summer maximum surface melt extent and long-term trends; 3) estimate and analyze the GrIS water balance components, including the SMB and freshwater runoff; and 4) quantify the GrIS freshwater runoff and accumulated runoff contribution to global sea level rise.
2. Study area
Greenland is the world’s largest island, and the GrIS the Northern Hemisphere’s largest terrestrial permanent ice- and snow-covered area (1.834 × 106 km2), which covers approximately 85% of the island (Fig. 1). Greenland is roughly 2600 km long, up to approximately 950 km wide, and the ice sheet’s maximum altitude is more than 3200 m MSL. The total ice sheet volume is 2.85 × 106 km3, equivalent to an average global sea level rise between 7.0 and 7.4 m SLE (Warrick and Oerlemans 1990; Gregory et al. 2004; Lemke et al. 2007).
The climate in Greenland is arctic (Born and Böcher 2001). In the northern parts of the GrIS, winter air temperatures can drop below −70°C, whereas on the east Greenland land strip, summer temperatures can briefly rise above 25°C (Mernild et al. 2008b). The observed GrIS mean annual air temperature (MAAT) is −13.3°C, covering a nonsignificant MAAT warming of ∼1.8°C for the period 1995–2005 (based on data from the 10 coastal meteorological stations, Fig. 1; Table 1, stations 16–25). In southern and southeastern Greenland, the observed annual precipitation is ∼2400 mm w.eq. yr−1, while the northern desertlike areas hardly receive any precipitation (<200 mm w.eq. yr−1) (e.g., Born and Böcher 2001; Serreze and Barry 2005). Many of the island’s characteristics cause considerable contrast in its weather conditions—including complex coastal topography, elevation, distance from the coastal area, marginal glaciers, and ice caps—and the GrIS, which makes the climate vary appreciably even over short distances. Temperature inversions are a common feature for Greenland coastal areas (Hansen et al. 2008; Mernild et al. 2008a) and for the GrIS (Putnins 1970).
3. Models and methods
Throughout the Arctic, rough terrain, harsh climatic conditions, and remote locales are commonly cited reasons for lacking knowledge and adequate data. Furthermore, logistical constraints make it difficult to collect extensive observations of snow, sublimation, evaporation, and snow and glacier melt. Collecting runoff measurements have typically been considered impossible. Also, scattered Arctic meteorological stations and limited winter and summer GrIS mass-balance observations have resulted in sparse and unreliable data related to the spatial and temporal distribution of snow precipitation, sublimation, surface snow and ice melt across much of the GrIS, and runoff to the ocean. Such key components are essential to hydrological research efforts, and there is a clear need to explore issues associated with data sparseness and modeling capabilities.
Likewise, there are several kinds of uncertainties related to climate projections using simulations with coupled atmosphere–ocean GCMs (AOGCMs). Apart from uncertainties in future greenhouse gas and aerosol emissions and their conversion to radiative forcings, there are uncertainties in global and, in particular, regional climate responses to these forcings owing to, for example, different parameterizations [discussed in detail by Stocker et al. (2001)]. There is also large natural variability on the regional scale (consider, e.g., the North Atlantic Oscillation), so it is difficult to determine which part of the response of a model is due to anthropogenic forcing and to natural variability (solar, volcanic, but also unforced), respectively.
This implies that there is no single “best” model to use in an assessment of Arctic (or Greenland) climate changes, although some models clearly perform better than others (e.g., Christensen et al. 2007b; Walsh et al. 2008). However, most of the uncertainties mentioned in the previous paragraph can be quantified by using ensembles of model simulations rather than one particular model. Here, we are limited by the availability of one realization of a downscaling experiment only, where the behavior of the driving GCM model is amazingly realistic in describing present-day conditions in the Arctic in general and around Greenland in particular (Walsh et al. 2008; Stendel et al. 2008).
a. Choosing the model configuration
To assess how representative our results may be, given our particular regional–global model setup, we note as an important starting point that the ECHAM5 GCM is one of the best performing state-of-the-art models when it comes to representing the present climate. Walsh et al. (2008) in their study of GCMs participating in the most recent Coupled Model Intercomparison Project (CMIP) found that for Greenland, the Arctic, and the Northern Hemisphere, the ECHAM5 model outperforms most other IPCC class models with respect to all three parameters studied: temperature, precipitation, and mean sea level pressure. In the presence of sea ice, snow coverage, and frozen grounds, the interpretation of a climate change signal from GCM simulations is very sensitive to the realism of simulated present-day conditions (e.g., see Christensen et al. 2007b, 2008). Visual inspection of temperature bias maps in Walsh et al. (2008, Fig. 8) documents that sea ice coverage for present-day conditions in the Arctic is very well depicted by ECHAM5, whereas most other GCMs tend to show severe biases reflecting either too much or too little sea ice in many locations of the Arctic. In Christensen et al. (2007b), the climate change signal for 21 CMIP models (including ECHAM5) are displayed. As pointed out in Christensen et al. (2008), it is interesting to note, however, that the projected climate signals to some degree are caused by quite different mechanisms. Observing in an extract of the performance of the model from Walsh et al. (2008), they found a common feature for most of the models, reflected by the ensemble mean, of a clear tendency to simulate too much sea ice in the Barents Sea in winter, with the ECHAM5 model being a clear exception. At the same time, the greatest warming by the end of the century is simulated exactly over this region in the ensemble mean as well as by the individual models. In the National Center for Atmospheric Research and Geophysical Fluid Dynamics Laboratory models, for example, this is partly reflecting the present-day bias, while in the ECHAM5 model, this apparently cannot be the case. Furthermore, in general, the largest warming occurs in the area with too much ice (strong cold bias) under present-day conditions. Thus to some extent, the results at the regional scale are clearly subject to systematic errors in present-day simulations. Using an ensemble of models masks this deficiency. Therefore, maps of warming must be carefully analyzed and cannot be used at face value in a region with nonlinear feedbacks, such as in the presence and absence of sea ice.
Given these caveats, we are confident that ECHAM5-based simulations are as good as any possible simulation based on a random choice of a GCM. If all GCMs were used for downscaling, the resulting distribution would partly (perhaps even largely) be due to outlier models with a poor representation of Arctic and Greenland climate (especially sea ice) conditions. It is beyond the scope of the present paper to quantify in more detail the uncertainty of our results owing to of the GCM–RCM choice made here.
b. HIRHAM4 RCM
The IPCC A2 and B2 climate scenarios (available online at http://www.ipcc.ch) have been used in the HIRHAM4 RCM (Christensen et al. 1996, 2001; Bjørge et al. 2000; Christensen and Christensen 2007) for several climate change time-slice experiments for present and future conditions, with the third climate configuration of the Met Office Unified Model (HadCM3) AOGCM (Gordon et al. 2000; Pope et al. 2000) or ECHAM4 AOGCM (Roeckner et al. 1996; Christensen et al. 2007a,b) as boundary conditions. These simulation areas cover central and northern Europe. The performance of the model for arctic conditions has been found to be state of the art in many aspects (e.g., Christensen and Kuhry 2000; Dethloff et al. 2002; Kiilsholm et al. 2003). In a recent HIRHAM4 study Stendel et al. (2008) has set up the model to conduct a transient climate change experiment representing the period 1950–2080 for the IPCC scenario A1B, covering Greenland and adjacent sea areas (Fig. 1). This A1B scenario was run on a 25-km grid cell increment with 19 vertical levels, using boundary conditions from ECHAM5 AOGCM (Roeckner et al. 2003). While high-resolution regional climate simulations to date mainly have been run as time-slice experiments, we present results of a transient simulation covering 1950–2080. All forcing data have been taken from the transient ECHAM5–Max Planck Institute Ocean Model (MPI-OM1) run. The A1B experiment, as described in the IPCC Fourth Assessment Report (AR4), runs (runs were done in a worldwide coordination with a clearly defined common model setup) begin in 2000, and the AOGCM uses outputs from a detailed simulation of the twentieth-century experiment as initial conditions in 2000 (e.g., Randall et al. 2007).
Of course, it would be desirable to investigate an ensemble of RCM simulations [different RCMs forced by different GCMs or at realizations of a single GCM, as done in the PRUDENCE project, see Christensen et al. (2007a)]. This was, however, impossible owing to the lack of computer capacity, so we had to restrict ourselves to this particular configuration. The results from Walsh et al. (2008) indicate that the chosen GCM is a sensible choice when only one realization can be offered. Déqué et al. (2007) indicates that the role of choosing a particular ensemble member is insignificant compared with choosing the GCM.
c. SnowModel description
SnowModel (Liston and Elder 2006a) is a spatially distributed snowpack evolution modeling system specifically designed to be applicable over the wide range of snow landscapes, climates, and conditions found around the world. It is made up of four submodels: MicroMet defines the meteorological forcing conditions (Liston and Elder 2006b); EnBal calculates the surface energy exchanges, including melt (Liston 1995; Liston et al. 1999); SnowPack simulates snow depth and water equivalent evolution (Liston and Hall 1995); and SnowTran-3D is a blowing-snow model that accounts for snow redistribution by wind (Liston and Sturm 1998, 2002; Liston et al. 2007). Although other distributed snow models exist (e.g., Tarboton et al. 1995; Marks et al. 1999; Winstral and Marks 2002), the SnowTran-3D component allows for application in arctic, alpine (i.e., above treeline), and prairie environments that make up 68% of seasonally snow-covered areas in the Northern Hemisphere (Liston 2004). SnowModel also simulates snow-related physical processes at spatial scales ranging from 5 m to global and temporal scales ranging from 10 min to a whole season. Simulated processes include 1) accumulation and loss from snow precipitation, and blowing-snow redistribution; 2) loading, unloading, and blowing-snow sublimation; 3) snow density evolution; and 4) snowpack ripening and melt. SnowModel was originally developed for glacier-free landscapes. For glacier surface mass-balance studies, SnowModel was modified to simulate glacier ice melt after winter snow accumulation had ablated (Mernild et al. 2006b, 2007).
MicroMet is a quasi-physically based meteorological distribution model (Liston and Elder 2006b) specifically designed to produce the high-resolution meteorological forcing distributions (air temperature, relative humidity, wind speed, wind direction, precipitation, solar and longwave radiation, and surface pressure) required to run spatially distributed terrestrial models over a wide range of landscapes in a physically realistic manner. MicroMet uses elevation-related interpolations to modify air temperature, humidity, and precipitation following Kunkel (1989), Walcek (1994), Dodson and Marks (1997), and Liston et al. (1999). Temperature and humidity distributions are defined to be compatible with the observed lapse rates. Wind flow in complex topography is simulated following Ryan (1977) and Liston and Sturm (1998). Solar radiation variations are calculated using elevation, slope, and aspect relationships (Pielke 2002). Incoming longwave radiation is calculated while taking into account cloud cover and elevation-related variations following Iziomon et al. (2003). Precipitation is distributed following Thornton et al. (1997). In addition, any data from more than one location, at any given time, are spatially interpolated over the domain using a Gaussian distance-dependent weighting function and interpolated to the model grid using the Barnes objective analysis scheme (Barnes 1964, 1973; Koch et al. 1983). Liston and Elder (2006b) and Liston et al. (2007) have performed a rigorous validation of MicroMet using various observational datasets, data denial, and geographic domains. Furthermore, MicroMet has been used to distribute observed and modeled meteorological variables over a wide variety of landscapes in the United States: Colorado (Greene et al. 1999), Wyoming (Hiemstra et al. 2002, 2006), Idaho (Prasad et al. 2001), and arctic Alaska (Liston et al. 2002, 2007; Liston and Sturm 1998, 2002); Norway: Svalbard and central Norway (Bruland et al. 2004); east Greenland (Mernild et al. 2006b, 2007); the Greenland ice sheet (Mernild et al. 2008c,d, 2009a,b); and near-coastal Antarctica (Liston et al. 1999; Liston and Winther 2005). As an example, the GrIS validations of MicroMet-simulated meteorological data indicate substantial correlation with independently observed GrIS meteorological data from, for example, the Swiss camp [located within 50 km from JAR1 (station 7) at 1140 m MSL (Table 1, Fig. 1)]. MicroMet-generated air temperature, relative humidity, and precipitation values account for 84%, 63%, and 69%, respectively, of the variance in the observed 1995–2005 daily averaged dataset. The wind speed has less strong correlations, but the results remain respectable (>50% variance) for representations of GrIS meteorological processes (Mernild et al. 2008d).
EnBal performs standard surface energy balance calculations (Liston 1995, Liston et al. 1999). This component simulates surface (skin) temperatures and energy and moisture fluxes in response to observed and/or modeled near-surface atmospheric conditions provided by MicroMet. Surface latent and sensible heat fluxes and snowmelt calculations are made using a surface energy balance model of the form:
where Qsi is the solar radiation reaching the earth’s surface, Qli is the incoming longwave radiation, Qle is the emitted longwave radiation, Qh is the turbulent exchange of sensible heat, Qe is the turbulent exchange of latent heat, Qc is the conductive energy transport, Qm is the energy flux available for melt, and α is the surface albedo. Details of each term in Eq. (1) and the model solution are available in Liston (1995) and Liston et al. (1999). In the presence of snow or glacier ice, surface temperatures greater than 0°C indicate that energy is available for melting. This energy is computed by fixing the surface temperature at 0°C and solving Eq. (1) for Qm. Energy transports toward the surface are defined to be positive.
SnowPack is a single-layer snowpack evolution and runoff/retention model that describes snowpack changes in response to precipitation and melt fluxes defined by MicroMet and EnBal (Liston and Hall 1995; Liston and Elder 2006a). Its formulation closely follows Anderson (1976). In SnowPack, the density changes with time in response to snow temperature and weight of overlying snow (Liston and Elder 2006a). A second density-modifying process results from snow melting. The melted snow reduces the snow depth and percolates through the snowpack. If the snow temperature is below freezing, any percolating/liquid water refreezes and is stored in the snow (in the “pores”) as internal refreezing. When saturated snow density is reached, assumed to be 550 kg m−3 (Liston and Hall 1995), actual runoff occurs. This provides a method to account for heat and mass transfer processes, such as snowpack ripening, during spring melt. The density of new snow from additional accumulation is defined following Anderson (1976) and Liston and Hall (1995). Static-surface (nonblowing snow) sublimation calculated in EnBal is used to adjust the snowpack depth; blowing-snow sublimation is calculated in SnowTran-3D (Liston and Elder 2006a).
SnowTran-3D (Liston and Sturm 1998; Liston et al. 2007) is a three-dimensional submodel that simulates snow depth evolution (deposition and erosion) resulting from windblown snow based on a mass-balance equation that describes the temporal variation of snow depth at each grid cell within the simulation domain. SnowTran-3D’s primary components are a wind flow forcing field, a wind shear stress on the surface, snow transport by saltation, snow transport by turbulent suspension, sublimation of saltating and suspended snow, and accumulation and erosion at the snow’s surface (Liston and Sturm 2002). Simulated transport and blowing-snow sublimation processes are influenced by the interactions among available snow, topography, and atmospheric conditions (Liston and Sturm 1998). SnowTran-3D simulates snow depth evolution and then uses the snow density simulated by SnowPack to convert it to the more hydrologically significant snow water equivalent (SWE) depth. Deposition and erosion, which lead to changes in snow depth [Eq. (2)], are the result of changes in horizontal mass transport rates of saltation Qsalt (kg m−1 s−1), changes in horizontal mass transport rates of turbulent suspended snow Qturb (kg m−1 s−1), sublimation of transported snow particles Qυ (kg m−2 s−1), and the water equivalent precipitation rate P (m s−1). Combined, the time rate of change in snow depth ζ (m) is
where t (s) is time; x (m) and y (m) are the horizontal coordinates in the west–east and south–north directions, respectively; and ρs and ρw (kg m−3) are snow and water density, respectively. At each time step Eq. (2) is solved for each individual grid cell within the domain and is coupled to the neighboring cells through the spatial derivatives (d/dx, d/dy). SnowTran-3D simulations have previously been compared against observations in glacier and glacier-free alpine, Arctic, and Antarctic landscapes (Greene et al. 1999; Liston et al. 2007; Prasad et al. 2001; Hiemstra et al. 2002, 2006; Liston and Sturm 2002; Bruland et al. 2004; Mernild et al. 2006b, 2007, 2008c,d, 2009a,b).
5) SnowModel input
To solve the equations, SnowModel requires spatially distributed fields of topography and land cover in addition to temporally distributed point meteorological data (air temperature, relative humidity, wind speed, wind direction, and precipitation). Meteorological data was obtained from the HIRHAM4 RCM model (1950–2080) based on the IPCC scenario A1B and from observations from meteorological stations located within the simulation domain (1995–2005). For this study, observed data are obtained from 25 meteorological stations: 15 stations from the GC-Net, 9 from the WMO station from the near coast operated by the DMI, and 1 by the Danish National Environmental Research Institute and the University of Copenhagen (Fig. 1 and Table 1). Simulations were performed on a one-day time step, although snow and ice melt and blowing snow are threshold processes that may not be accurately represented by this time step. We recognize that the use of daily averaged atmospheric forcing variables, instead of hourly values, will produce a smoothing of the natural system. Therefore, daily simulated melt (ablation) and blowing-snow processes (accumulation) were tested against hourly simulated ablation and accumulation values from a test area, the Mittivakkat Glacier (31 km2) in southeast Greenland (Mernild and Liston 2010), and remain significant (p < 0.01, where p is level of significance), with an average difference of 2%, 3%, and 8% for the glacier winter, summer, and net mass balances, respectively.
Snow precipitation measurements include uncertainties, especially under windy and cold conditions (e.g., Yang et al. 1999; Liston and Sturm 2002, 2004; Serreze and Barry 2005). Solid and liquid precipitation measurements at the DMI meteorological stations (Fig. 1 and Table 1, stations 16–18 and 20–25) were calculated from Helman–Nipher shield observations corrected according to Allerup et al. (1998, 2000). Solid (snow) precipitation was calculated from snow depth sounder observations at the other stations (Fig. 1 and Table 1) after sounder data noise was removed; these data are assumed to be accurate within ±10%–15% (Mernild et al. 2007, 2009b). Snow depth sounder observations were partitioned into liquid (rain) and solid (snow) precipitation at different air temperatures based on methods employed at Svalbard (Førland and Hanssen-Bauer 2003). For air temperatures below −1.5°C, sounder data were considered to represent solid precipitation and for temperatures above 3.5°C precipitation is considered liquid; linear interpolation was used to calculate snow and rain fractions at temperatures between these limits. Snow depth increases at relative humidity <80% and wind speed >10 m s−1 were removed to better distinguish between the proportions of real snow accumulation based on precipitation events and blowing-snow redistribution (Mernild et al. 2007, 2009b). Remaining snow depth increases were adjusted using a temperature-dependent snow density (Brown et al. 2003) and hourly snowpack settling.
Greenland topographic data for the model simulations were provided by Bamber et al. (2001), who applied “correction” elevations derived by satellite imagery to an existing radar-altimetry digital elevation model (DEM). The image-derived correction was determined from a high-resolution (625 m) grid of slopes inferred from the regional slope-to-brightness relationship of 44 Advanced Very High Resolution Radiometer images covering all of Greenland (Scambos and Haran 2002). For the model simulations, this time-invariant DEM was aggregated to a 5-km gridcell increment and clipped to yield a 2830 km × 1740 km simulation domain that encompassed all of Greenland. The GrIS terminus was confirmed or estimated by using aerial photos and maps (1:250 000, Geodetic Institute, Denmark).
SnowModel is a surface model producing first-order effects of climate change; it does not include glacio-hydrodynamic and sliding routines. Using a time-invariant DEM could be inappropriate. Therefore, a 1950–2080 assessment of GrIS volume, area, and maximum and average heights was performed using the simulation code for polythermal ice sheets (SICOPOLIS) (Greve 1997a,b, 2005), a state-of-the-art 3D dynamic/thermodynamic shallow-ice approximation model. On the basis of the IPCC A1B scenario for the period, there was a small change in the GrIS. By 2080 the volume differed 3% (5.01 × 104 km3), melt area changed 8% (4.88 × 104 km2), and heights shifted <1% (maximum height: 17 m and average height: 3 ± (8) m). These discrepancies fall well within the uncertainties of this study. Furthermore, for this study there is only a one-way nesting between HIRHAM4 (the atmosphere) and SnowModel (the surface), not taking into account, for example, the positive albedo feedback associated with snowmelt and the fact that wet snow absorbs as much as three times more incident solar energy than dry snow (Steffen 1995).
Each grid cell within the domains was assigned a U.S. Geological Survey land use/land cover system class according to the North American Land Cover Characteristics Database, Version 2.0 [available from the USGS Earth Resources Observatin and Science (EROS) Center’s Distributed Active Archive Center, Sioux Falls, South Dakota]. The snow-holding depth (the snow depth that must be exceeded before snow can be transported by wind) was assumed to be constant. Albedo was assumed to be 0.8 for snow (Table 2). Realistically, snow albedo changes with time and surface characteristics (Pomeroy and Brun 2001); thus, the model will likely underestimate energy available for surface melting. Therefore, SnowModel simulations with a fixed snow albedo of 0.8 was tested against simulated variable snow albedo from a test area, the Jakobshavn Isbræ drainage area in west Greenland [for information about the variable albedo routines see Mernild et al. (2009c)], indicating a mean annual variable snow albedo of 0.7 and a difference of up to ∼15% in SMB and runoff. When the snow is ablated, GrIS surface ice conditions are used. Ice albedo was invariant and assumed to be 0.4. The GrIS ablation zone is characterized by lower albedo on the margin and an increase in albedo toward the equilibrium line altitude (ELA), where a veneer of ice and snow dominate the surface (Boggild et al. 2006). The emergence and melting of old ice in the ablation zone 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 windblown sediment (Boggild et al. 2006). Particles on or melting into the ice change the area-average albedo, increasing melt. User-defined constants for SnowModel are shown in Table 2 [for parameter definitions see Liston and Sturm (1998, 2002)]. All fjord and ocean areas within the domain were excluded from model simulations. Further, changes in glacier storage based on changes in supraglacial storage (lakes, pond, channels, etc.), englacial storage (ponds and the water table), subglacial storage (cavities and lakes), meltwater routing, evolution of a runoff drainage system, and changes based on iceberg calving, tidal response where ice meets ocean, and geothermal melting are not calculated in the SnowModel simulations, even though they might influence the contribution of runoff.
d. Satellite images
Detection of surface melt at large spatial scales is effectively accomplished by using satellite microwave data. The daily GrIS snowmelt extent is mapped (25-km gridcell increment) using passive microwave satellite observations that discriminate wet from dry snow. The criterion for melt is 1% mean liquid water content by volume in the top meter of snow (Abdalati and Steffen 1997). The center part of the GrIS is the area where the melting threshold of the cross-well ground-penetrating radar microwave algorithm did not show any melt. The end-of-summer maximum observed spatial surface melt distribution at the GrIS was used to validate SnowModel melt simulations.
e. HIRHAM4 RCM validation and uncertainty
Before the HIRHAM4 RCM output data can be trusted for use as input data for further modeling, it needs to be tested and calibrated against observed meteorological data since RCM output biases can be large. Stendel et al. (2007) provide some basic validation of the current simulation. However, because HIRHAM is running in a full climate mode, that is, the driving GCM only knows about the state of the atmosphere–ocean system from the external drivers (sun, aerosols, and greenhouse gases), whether actually realized (1950–2000) or projected (2001–80) we need for our purpose to tie in this single realization with the observed climate system. We have excellent data for verifying SnowModel covering the period 1995–2005. For a bias adjustment or calibration of the HIRHAM results a 10-yr period is relatively short; however, we have assessed the role of this short period by an additional calibration in which the model years were 1980–90 and the observed years 1995–2005. The resulting offset in precipitation is on average 42 mm w.eq. (or 7%) and the temperature difference was 1.5°C (or 10%) for 1980–90 with respect to the calibration period 1995–2005. Relative humidity and wind are both insignificantly changed. Mean monthly offset between the RCM modeled output and the observed meteorological data were further estimated for the period 1995–2005 (see Fig. 2 and Table 1 for station information). These mean monthly (1995–2005) offset values were added to the daily RCM meteorological parameters to correct each variable (air temperature, relative humidity, wind speed, and corrected precipitation) for the 1950–2080 period before being used as meteorological forcings for SnowModel. To assess the performance of the adjusted SnowModel simulated spatial distributed meteorological data the spatial distributed meteorological data were tested against in situ meteorological observations (with data not used for calibration) spanning 1995–2005 (see Table 1 for stations used for calibration and validation). We have ranked the data for each period and compared the ranked numbers. This illustrates the ability of HIRHAM to capture the span of realized parameters for the period of concern and, therefore, also gives a rough estimate about the calibration method. Ideally we should use longer periods and address classical climatological values, but this is beyond the scope of the present work, and some of the results are provided elsewhere (e.g., Stendel et al. 2007). Validations of the simulated GrIS meteorological data (air temperature, relative humidity, and wind speed) indicate substantial correlation with in situ–observed meteorological data from different meteorological stations on the GrIS—JAR1, Humbolt, Saddle, and Summit—at different elevations and with in situ–observed precipitation from outside the GrIS—Station Nord, Danmarkshavn, Ittoqqortoormiit, and Ikerasassuaq—at different latitudes (see Fig. 3 and Table 1 for station information). Modeled air temperature values account for 98%–99% of the variance in the observed 1995–2005 mean monthly dataset. The relative humidity, corrected precipitation, and wind speed have the same or slightly fewer strong correlations, but results remain respectable for relative humidity (between 85% and 96%), wind speed (between 83% and 98%), and for precipitation (between 89% and 98%) for representations of the GrIS meteorological processes (Figs. 3a–3d).
The most obvious model bias is the systematic dry bias of the simulated near-surface humidity, particularly when humidity is high. This is related to a general difficulty of representing coastal climate on model land points in HIRHAM (Stendel et al. 2008). For precipitation, we note that, with the exception of Danmarkshavn, the model captures the span of observed precipitation rather well (perhaps surprisingly so), given the short period of comparison. For temperature, the annual cycle is clearly the dominant feature in the explained variance. However, we also note that the agreement over the full span of temperatures is within what seems to be acceptable for our purpose (see, however, Stendel et al. 2008).
f. SnowModel validation and uncertainty
Few quality observations for spatial in situ snow-evolution, snow and ice surface melt, and glacier net mass-balance are available in Greenland, including from the GrIS. SnowModel accumulation and ablation routines have been tested quantitatively [simulations based on observed meteorological data; for further information see Mernild et al. (2006a,b, 2007, 2008d)] at local scale (from east and west Greenland) and regional scale (from the GrIS) using observations from snow pit depths; glacier winter, summer, and net mass balances; depletion curves; photographic time lapses; satellite images (microwave satellite–derived melt extent); and different parameterizations such as melt index and ELA. A maximum discrepancy between modeled and observed SWE depths of 7%, glacier mass balances of 7%, snow cover extent of 7%, and GrIS melt discrepancy between melt and nonmelt boundaries of 32(±24) km occurs (Mernild et al. 2009a,b). However, in northeastern Greenland, the discrepancy can be up to 160 km where the distances among meteorological stations is far (Fig. 1). In this study, SnowModel-simulated melt extent were compared against concurrent passive microwave satellite–derived melt extent and previous GrIS SMB studies.
SnowModel, like all models, possesses uncertainties owing to processes not represented by the modeling system. For example, routines for simulating the air temperature inversion layer and variable snow and ice albedo are not yet included. In addition, changes in the GrIS area, size, and height according to glacier dynamical processes and subglacial geothermal bottom melting and sliding are not calculated in the model routines. Based on the uncertainties in the modeled results from previous Greenland SnowModel simulations, including the GrIS, it is reasonable to assume that this GrIS SMB study is influenced with a similar maximum uncertainty of 7% for SWE depth, snow cover extent, and SMB components (Mernild et al. 2006a,b, 2007, 2008d).
4. Results and discussion
a. Regional climate model trends 1950–2080
The regional climate model adjusted meteorological data for the 1950–2080 GrIS (air temperature, relatively humidity, wind speed, and precipitation) are illustrated in Figs. 4a–d. The GrIS is divided into four subareas (I–IV, Fig. 4a). The greatest changes in predicted MAAT of 5.6°C occurs in northeast Greenland (area I) (significant, p < 0.01); this is likely due to the projected change in sea ice extent and thickness, particularly off the east coast of Greenland. The lowest warming, 3.6°C, occurs in area III (significant p < 0.01), southwest Greenland (Fig. 4a), where sea surface temperatures are changing only marginally (see Stendel et al. 2008). Overall, MAAT increased a significant 4.8°C from 1950 through 2080 (Fig. 4a). Patterns of annual minimum average (1981) and annual maximum average (2080) temperature distribution (Fig. 5a) show the variable extent of low interior temperatures (<−20°C) and higher temperatures in coastal Greenland (>0°C). Temporally, the average change in summer (June–August) temperature—temperatures affecting the ablation processes—is 3.1°C (significant, p < 0.01). Six of the coldest summers occurred in the first decade (1950–59), while the six warmest summers were in the last decade (2070–80) of the simulations (Table 3). We note here that this behavior underlines the general aspect of the simulation, namely, that the steady warming despite decadal variations is quite robust, which is not likely to be altered if another ECHAM5 ensemble member had been chosen. The average change in summer temperature (3.1°C) is below the average change in MAAT (4.8°C). A winter (December–February) average change of 5.9°C (significant, p < 0.01) is sizeable. Identical seasonal trends were identified in observations by Box (2002) and Sturm et al. (2005) from the 1970s through 1990s.
From 1950 to 2080 relative humidity increased 1.2% on average (significant, p < 0.01) (Figs. 4b and 5b). Average wind speed decreased slightly, <0.1 m s−1 (insignificant, p < 0.25); the largest reduction, −0.2 to −0.5 m s−1, occurs on the GrIS northeastern interior (Figs. 4c and 5c). Modeled precipitation increased 80 mm w.eq. on the GrIS (significant, p < 0.01), with the lowest gain of 57 mm w.eq. in northwest Greenland (area IV) and the greatest increase of 160 mm w.eq. in southeast Greenland (significant, p < 0.01) (area II, Figs. 4d and 5d) due to projected changes in cyclonic systems. The overall trend for predicted climate (1950–2080) is warmer and wetter, where MAAT will increase from −14.8° (1950–59) to −10.1°C (2070–80) and average precipitation from 600 (1950–59) to 770 mm w.eq. yr−1 (2070–80) (Table 4).
b. SnowModel melt extent simulations
The simulated end-of-summer GrIS melt extent is illustrated in Figs. 6a –6c. To examine annual melt extent spatial variation, 1996 and 2007 were selected randomly from the observation period (1979–2007) and the assessment indicates a reasonable degree of similarity between the observed (passive microwave) and modeled melt distributions (Fig. 6a). The differences among spatial annual simulated and observed GrIS melt boundaries average 51(±34) km with a maximum distance of ∼180 km. Modeled end-of-summer maximum melt extent from the observation period are, on average, overestimated by ∼10% (Table 4) when compared with satellite observations. The interannual discrepancy, likely due to a fixed albedo for snow and ice and a mismatch in modeled and observed resolutions, ranges from ∼303 600 km2 (∼17%) in 1998 to ∼7200 km2 (∼1%) in 1991.
The GrIS-simulated surface melt and nonmelt extent are further shown on a decadal basis for the period 1950–59 through 2070–80 (Fig. 6b). The average 1950–59 end-of-summer melt extent was 30% (0.542 × 106 km2) and 56% for 2070–80 (1.025 × 106 km2), indicating an average maximum difference of ∼90% (0.483 × 106 km2). The greatest difference in melt extent occurs in the southern part of the GrIS. To the northwest (area IV) and northeast (area I) of the GrIS, the changes in melt extent are less pronounced (Fig. 6b). Further, for 1950–59 and 2070–80, surface melt occurred at elevations as high as 2550 and 3050 m MSL, respectively. The distribution of the amount of simulated melt days is further shown for the periods 1950–59 and 2070–80 (Fig. 6b), indicating a significant average increase of 28 melt days for the GrIS (R2 = 0.74, p < 0.01). For the period 1950–59 the maximum number of melt days was 126 increasing to 242 for 2070–80. The greatest number of melt days is seen in the southeastern part of the GrIS (area II). The largest change (∼50%–70%) in the number of melt days was visible in the eastern part (areas I and II) of the GrIS and is lower (∼20%–30%) to the west (areas III and IV) for 1950–2080 (Fig. 6b). The reason is likely the projected change in sea ice extent and thickness in adjacent seas.
A time series of the simulated end-of-summer GrIS surface melt extent from 1950 through 2080 is illustrated in Fig. 6c. The percentage of total modeled melt extent is shown for four years: 1961, the year with the lowest melt extent in the simulation period (1950–2080); 1983, lowest melt extent since the satellite observations began in 1979; 2007, greatest melt extent since the satellite observations began; and 2077, the year with the highest melt extent in the simulation period. Simulated melt extent varies from 0.389 × 106 km2 (21% of the total GrIS area) in 1961 to 1.204 × 106 km2 (66%) in 2077, indicating an increasing GrIS melt extent through the period.
c. Water balance components
Throughout the year, surface processes such as snow accumulation and redistribution, evaporation, sublimation (including blowing-snow sublimation), and surface melt affect the GrIS water balance [Eq. (3)]. The yearly water balance equation for the GrIS can be described by
where P is the precipitation input from snow and rain (and possible condensation), E is evaporation (liquid to gas phase (atmosphere) flux of water vapor), SU is sublimation (snow blowing; solid to gas phase with no intermediate liquid stage), R is runoff, and ΔS is change in storage (ΔS is also referred as SMB) from changes in glacier storage and snowpack storage. 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 as the residual value.
The RCM–SnowModel SMB precipitation for 1995–2004 falls within the range of other studies (Table 5). The greatest average difference is 15 km3 yr−1: it is not surprising given vast uncertainties in measuring snow precipitation. Measuring snow precipitation typically includes errors, especially under windy and cold conditions (e.g., Yang et al. 1999; Liston and Sturm 2002, 2004; Serreze and Barry 2005). Snowfall in the Arctic is most often connected with strong winds and typically takes the form of fine snowflakes (Sturm et al. 1995). As a result, wind easily lifts and redistributes the snowflakes according to exposure and local topography, and it is sometimes difficult to distinguish between a period of snowfall and a period of drifting snow.
RCM–SnowModel simulated GrIS runoff estimate (1995–2004, Table 5) was highest compared to the other studies. The maximum difference was 59 km3 yr−1 and the minimum difference was 3 km3 yr−1. SnowModel runoff routines take retention and internal refreezing into account when meltwater penetrates through the snowpack. These routines do have a significant effect on the SMB runoff. The role of meltwater retention in terms of the overall GrIS mass balance indicates that runoff is overestimated between 20% and 29% if no retention/refreezing routines are included (1995–2004) (Mernild et al. 2008c). The overestimation corresponds with previous values of ∼25% estimated by the Janssens and Huybrechts (2000) single-layer snowpack model (used by, e.g., Hanna et al. 2002, 2005, 2008; Table 5). The lack of retention/refreezing routines in SnowModel [used in this paper and Mernild et al. (2008d), Table 5] leads to an overestimation of ocean runoff, a consequent overestimation of global sea level rise, and may explain the larger difference among SnowModel simulated runoff and the other studies.
For SMB (1995–2004), the average RCM–SnowModel simulated values were lowest, 32 km3 yr−1 lower than Mernild et al. (2008d; a study based on observed data only) and 51 km3 yr−1 lower than Box et al. (2006). Compared with the study by Hanna et al. (2008), the RCM–SnowModel SMB was 207 km3 yr−1 lower. The lower SnowModel-simulated GrIS SMB values are due to the incorporation of evaporation and sublimation values of 142 km3 yr−1 in the SMB calculations [see Eq. (1)], where sublimation alone accounts for 64 km3 yr−1 on average: a value in the lower end of the Box and Steffen (2001) observed GrIS sublimation values of 62(±23) to 120(±65) km3 yr−1 (1995–2000). SnowModel simulated evaporation and sublimation accounted for 26% of the total GrIS ablation processes, indicating variations in the range from 134 in 2003 to 153 km3 yr−1 in 1999.
Table 4 presents the decadal GrIS surface melt conditions and the water balance components [Eq. (3)] for the period 1950–2080. The melt index (the area above 2000-m elevation where the greatest changes in melting occur) increased 138% (1.96 × 106 km2 × days), and the end-of-summer maximum melt extent grew 89% (4.83 × 105 km2). The trend in melt extent is illustrated in Fig. 6c. Over time, more of the GrIS surface area melted increasing from 0.542 × 106 (1950–59) to 1.025 × 106 km2 (2070–80), and the melting occurred for a longer duration during the ablation season. Increasing decadal temperatures largely explain the variance between MAAT and melt extent (R2 = 0.79) and melt index (R2 = 0.89), indicating that rising temperatures influence the ablation processes and melt conditions.
Modeled ELA provides a useful metric of accumulation and ablation’s net influence on the SMB (Table 4). On the GrIS from 1950–2080, the decadal ELA is changing in elevation from 1150 to 2060 m MSL, an average elevation increase of ∼70 m MSL decade−1. Values of ELA correlate highly with MAAT (R2 = 0.94, p < 0.01) and with precipitation (R2 = 0.73, p < 0.01) (Table 4). Location of the ELA is largely tied to changes in MAAT and subsequent changes in melt extent and melt index. The spatial location of ELA is influenced by local topography, regional variations in precipitation regimes, dominant cyclonic systems, and latitude.
The SMB trend from 1950 to 2080 (Fig. 7; Table 4) integrates accumulation (snow precipitation) and ablation (evaporation, sublimation, and runoff) over the GrIS. It is the manifestation of increased precipitation and ablation.
Interannual variability in precipitation and ablation caused sizeable SMB fluctuations with correlations of R2 = 0.46, p < 0.01 and R2 = 0.93, p < 0.01, respectively (Table 4). SMB fluctuations were largely tied to changes in ablation processes. Fluctuation patterns illustrated in Fig. 7, which were almost identical to the trends described by Rignot et al. (2008), indicated the highest balance in the 1970s and early 1980s with subsequent rapid losses as temperatures warmed. In Table 4 the interdecadal trend and variability in precipitation (R2 = 0.73, p < 0.01), evaporation/sublimation (R2 = 0.85, p < 0.01), runoff (R2 = 0.94, p < 0.01), and SMB (R2 = 0.86, p < 0.01) rates possessed significantly high correlations throughout the simulation period. Precipitation rose 133 km3, evaporation/sublimation 73 km3, and runoff 391 km3—leading to enhanced average SMB losses of 331 km3 (Fig. 7). Throughout the simulation period the decadal SMB varied from 256 (1960–69) to −99 km3 yr−1 (2070–80), averaging 79(±129) km3 yr−1 (Table 4). SMB values below zero (negative SMB value) occur from the period 2040–49 through 2070–80 (Table 4). A negative SMB developed in response to high ablation values (Table 4), ranging from an average of 706 km3 yr−1 (of which 74% was runoff) for 2040–49 to 870 km3 yr−1 (of which 77% was runoff) in 2070–80. Our SMB is similar to the Fettweis et al. (2008) SMB loss estimates generated from mean SMB values of 24 AOGCMs (using projections of temperature and precipitation anomalies from AOGCMs) performed for the IPCC Fourth Assessment Report for 2010–80. The RCM–SnowModel simulated SMB is, on average, ∼180 km3 yr−1 below the AOGCMs mean values and is similar to the lowest AOGCMs 2080 projected SMB of −100 km3 yr−1.
Sublimation plays an important role in the annual high-latitude hydrological cycle. Previous GrIS studies (e.g., Box and Steffen 2001; Mernild et al. 2008d) have shown that as much as 12%–23% of the annual precipitation may be returned to the atmosphere by sublimation. In arctic North America, studies by Liston and Sturm (1998, 2004), Essery et al. (1999), and Pomeroy and Essery (1999) indicate that 5%–50% of the annual solid precipitation was returned to the atmosphere by sublimation. For the GrIS (1950–2080), modeled annual sublimation averaged 74 km3 yr−1, which is ∼47% of the total average for evaporation and sublimation of 158 km3 yr−1 (Table 4) and is 12% of the total precipitation of 677 km3 yr−1. SnowModel simulated results were in the lower end of the Box and Steffen (2001) observed GrIS sublimation values of 62(±23) to 120(±65) km3 yr−1 (12%–23% of the total precipitation), even though the observed values are from the period 1995–2000.
The average GrIS runoff from the period 1950 through 2080 is 442(±134) km3 yr−1 (Table 4, Fig. 7). During this time runoff accelerated ∼30 km3 decade−1 to a runoff value of 668 km3 yr−1 (2070–80). The average GrIS runoff of 442 km3 yr−1 is comparable to approximately 1000 icebergs (density 917 kg m3) with dimensions 1 km × 1 km and an ice thickness 500 m. The GrIS runoff equals a specific runoff of 7.6(±2.3) l s−1 km−2 yr−1, equivalent to an average rise in global eustatic sea level of 1.2 mm SLE yr−1, changing from 0.8 ± 0.1 (1950–59) to 1.9 ± 0.1 mm SLE yr−1 (2070–80) (Table 4 and Fig. 7). The accumulated GrIS freshwater runoff is 160 mm SLE from 1950 through 2080. In addition to enhanced runoff, GrIS may shed mass by iceberg calving and geothermal melting. Thus, simulated GrIS freshwater runoff might underestimate the mass lost by half (Lemke et al. 2007; Mernild et al. 2008c).
In terms of our general satisfaction with these model results, it is important to be clear about the assumptions and potential deficiencies of this modeling study. In these simulations we have assumed a mean monthly offset value corrected to each meteorological variable, a time-invariant DEM, a fixed albedo for snow and ice, and no routines for the influence of air temperature inversions on snowmelt and glacier mass-balance simulations. We also recognize that the use of daily averaged atmospheric forcing variables will produce a smoothing of the natural system compared with higher temporal resolutions. Our understanding of the GrIS freshwater flux to the ocean is still far from complete. Detailed climate–cryospheric interactions are being examined at finer scales at the GrIS Kangerluassuaq drainage area, west Greenland, to estimate the freshwater influx to the ocean before upscaling routines to the entire GrIS.
5. Summary and conclusions
These GrIS simulations reveal continued warming and dramatically increased ablation amount and extent from 1950 to 2080. Over the period of simulation, surface runoff increased from 285 (1950–59) to 668 km3 yr−1 (2070–80). The GrIS freshwater runoff will be a factor in global sea level rise, equivalent to an average rise of 1.2 mm SLE yr−1, and a cumulative increase of around 160 mm SLE in this particular model setup under an IPCC A1B emission scenario.
Realistic simulations are required to better predict GrIS SMB loss and the impacts of this loss for the North Atlantic Ocean since it plays an important role in determining the global thermohaline circulation, salinity, sea ice dynamics, and the global eustatic sea level rise. There is a high degree of agreement between GrIS simulations and recorded observations as well as simulated GrIS SMB values and previous modeling studies. However, SnowModel does not yet include routines for variable snow albedo and for the influence of air temperature inversions on snowmelt and glacier mass-balance simulations. These improvements are forthcoming and will likely bolster modeling efforts. In this work, we have not considered feedback processes from a changing GrIS to the atmosphere, which are also likely to influence simulated surface air temperatures and thereby impact the resulting melt rates.
Another uncertainty that we have partly ignored here is the spread in model projections of the climate of the future. We acknowledge that more than 20 IPCC-type GCMs have been analyzed with respect to their projection in climate change, that is, by the IPCC (Christensen et al. 2007b), showing a wide range of results with the Arctic exhibiting an even higher lack of confidence than any other region. We attribute a substantial part of this uncertainty to imperfections of various models, particularly with respect to the representation of arctic processes. In our work, we employed only one model, which we have identified as one of (if not) the best GCMs in representing present climate conditions, ECHAM5. This model simulates a future climate not far diverged from the ensemble mean of 21 IPCC class models. Our results presented here are representative of state-of-the-art modeling, but are not comprehensive in assessing the entire range of possibilities.
Very special thanks to the three anonymous reviewers for their insightful critique of this article. This work was supported by grants from the University of Alaska Presidential IPY Postdoctoral Foundation and the University of Alaska Fairbanks (UAF) Office of the Vice Chancellor for Research and conducted during the first author’s IPY Post Doctorate Program at the UAF. Special thanks to the faculty of science and Institute of Low Temperature Science, Hokkido University, Japan, for hosting the first author from April through July 2008 and to Colorado State University, Cooperative Institute for Research in the Atmosphere, for hosting the first author through September and October 2008, and February 2009. Thanks to Martin Stendel and John Cappelen at the Danish Meteorological Institute (DMI) for providing HIRHAM4 RCM data and WMO meteorological data obtained from the near coast, and to the Cooperative Institute for Research in Environmental Science (CIRES), University of Colorado at Boulder, for providing meteorological data from the Greenland Climate Network (GC-Net) automatic weather stations. Also, thanks to Theodore Scambos, CIRES, University of Colorado at Boulder, for providing the Greenland digital elevation model.
* Current affiliation: Climate, Ocean, and Sea Ice Modeling Group, Computational Physics and Methods, Los Alamos National Laboratory, Los Alamos, New Mexico.
Corresponding author address: Dr. Sebastian H. Mernild, Climate, Ocean, and Sea Ice Modeling Group, Computational Physics and Methods (CCS-2), Los Alamos National Laboratory, Mail Stop B296, Los Alamos, NM 87545. Email: firstname.lastname@example.org