The modeling of the surface mass balance (SMB) of the Greenland Ice Sheet (GIS) requires high-resolution models in order to capture the observed large gradients in the steep marginal areas. Until now, global climate models have not been considered suitable to model ice sheet SMB owing to model biases and insufficient resolution. This study analyzes the GIS SMB simulated for the period 1850–2005 by the Community Earth System Model (CESM), which includes a new ice sheet component with multiple elevation classes for SMB calculations. The model is evaluated against observational data and output from the regional model Regional Atmospheric Climate Model version 2 (RACMO2). Because of a lack of major climate biases, a sophisticated calculation of snow processes (including surface albedo evolution) and an adequate downscaling technique, CESM is able to realistically simulate GIS surface climate and SMB. CESM SMB agrees reasonably well with in situ data from 475 locations (r = 0.80) and output from RACMO2 (r = 0.79). The simulated mean SMB for 1960–2005 is 359 ± 120 Gt yr−1 in the range of estimates from regional climate models. The simulated seasonal mass variability is comparable with mass observations from the Gravity Recovery and Climate Experiment (GRACE), with synchronous annual maximum (May) and minimum (August–September) and similar amplitudes of the seasonal cycle. CESM is able to simulate the bands of precipitation maxima along the southeast and northwest margins, but absolute precipitation rates are underestimated along the southeastern margin and overestimated in the high interior. The model correctly simulates the major ablation areas. Total refreezing represents 35% of the available liquid water (the sum of rain and melt).
The Greenland Ice Sheet stores a volume of water that, if released to the oceans, would raise global sea level on average by 7.3 m (Bamber et al. 2001). Recent observations show that the Greenland and Antarctic Ice Sheets are losing mass at an increasing rate (Zwally et al. 2005; Shepherd and Wingham 2007; Rignot et al. 2008a; Velicogna 2009; van den Broeke et al. 2009). A substantial contribution to this mass loss comes from the acceleration and thinning of large outlet glaciers and ice streams in Greenland and West Antarctica (Joughin et al. 2003; Scambos et al. 2004; Rignot and Kanagaratnam 2006). Also, the volume and extent of surface melting in Greenland have increased during recent decades (Steffen et al. 2004; Hanna et al. 2005, 2008; van den Broeke et al. 2009). These observations show that ice sheets respond to atmosphere and ocean warming rather quickly, on decadal or annual time scales.
It is generally accepted that global climate models are not yet ready to properly simulate ice sheet surface mass balance (SMB), owing to climate biases and insufficient horizontal resolution (van den Broeke et al. 2008a). Regional climate models (RCMs) are the state of the art in the modeling of ice sheet SMB (Box et al. 2006; Ettema et al. 2009; Fettweis 2007; Rae et al. 2012; Lenaerts and van den Broeke 2012; Lenaerts et al. 2012), but they suffer from their dependency on a GCM as external forcing. In most cases, RCMs do not update the topography and glacier mask of the GCMs, and there is often a problem of inconsistency between SMB and topography in long-term runs in comparison to earth system models that include interactive ice sheets. Also, when RCM output is used to force ice sheet models, the SMB is calculated at a fixed topography that can be quite different from the evolving topography of the ice sheet in the host model. While other components of the earth system, such as ocean dynamics, land vegetation, and the carbon cycle, are already standard components of global climate models; ice sheets are not. The only studies with an atmosphere–ocean GCM fully coupled with an ice sheet model are those performed with the Hadley Centre Coupled Model, version 3 (HadCM3) (Ridley et al. 2005), ECHAM3/Large-Scale Geostrophic ocean model (LSG) (Mikolajewicz et al. 2007; Vizcaíno et al. 2008), and ECHAM5/Max Planck Institute Ocean Model (MPI-OM) (Vizcaíno et al. 2010). The first two models make use of anomalies for the atmospheric forcing (for the fields of temperature and/or precipitation) and temperature index methods to calculate melt. Only the latter does not require the use of anomalies and uses an energy balance scheme with a simple parameterization for albedo.
The Community Earth System Model (CESM) (Hurrell et al. 2013) now includes a new ice sheet component. In CESM the surface mass balance of glacier surfaces is calculated with an energy balance scheme that includes a sophisticated modeling of surface albedo, as well as percolation and refreezing of meltwater within the snowpack. For this study, the SMB calculation is done in the CESM land surface component, the Community Land Model (CLM), at a horizontal resolution of ~1° at several fixed elevations. Height corrections to temperature and specific humidity are applied to account for differences in atmospheric forcing that are linked to elevation.
In this study we evaluate the performance of the CESM over the Greenland Ice Sheet for the period 1850–2005. We do so by comparing CESM output to that of a state-of-the-art regional climate model, the Regional Atmospheric Climate Model version 2 (RACMO2), and observational data from in situ stations and the Gravity Recovery and Climate Experiment (GRACE). RACMO2 is especially adapted for use over ice sheets. The RACMO2 simulated SMB shows very good agreement with independent in situ SMB observations (N = 265, correlation coefficient r = 0.95) (Ettema et al. 2009). The simulated SMB time series combined with ice flux data (Rignot et al. 2008b) are highly correlated with satellite data from GRACE (Velicogna 2009), with r = 0.99 (van den Broeke et al. 2009).
The outline of this paper is as follows. Section 2 gives a description of the model and the setup of the simulations. Section 3 analyzes the simulated near-surface climate of Greenland and compares it with RACMO2. Section 4 evaluates the surface mass balance of the reference period 1960–2005 by comparison with in situ observations, GRACE data, and output from RACMO2. Section 5 presents the climate and surface mass balance evolution for the whole period of the historical simulation (years 1850–2005). The last section presents a summary and conclusions.
The fully coupled model CESM1.0 is used for these simulations. The Community Atmosphere Model version 4 (CAM4) (Neale et al. 2013) is run with the finite volume dynamical core at a horizontal resolution of 0.9° × 1.25°, which is also the resolution of the land model, CLM (Oleson et al. 2010). CAM4 discretization is done via finite volumes. The horizontal resolution of the Parallel Ocean Program (POP) ocean model (Smith et al. 2010) and the Los Alamos Sea Ice Model (CICE) (Hunke 2010; Jahn et al. 2012) is approximately 1°, and the horizontal resolution of the ice sheet component, the Glimmer Community Ice Sheet Model (Glimmer-CISM) (Rutt et al. 2009), is 5 km. Since the purpose of this paper is to evaluate only the surface mass balance, the ice sheet topography has been kept fixed at present-day measurements (Bamber et al. 2001). Ice sheet flow has not been calculated here but is evaluated in a companion paper by Lipscomb et al. (2013).
The surface mass balance is calculated in the land component (CLM) with an energy balance scheme that includes a sophisticated modeling of snow albedo, the Snow and Ice Aerosol Radiation model (SNICAR) (Flanner and Zender 2006). Each CLM grid cell can contains one or more land surface types or land unit. A new “glacier_mec” land unit type (where “mec” denotes multiple elevations classes) has been added to five previous land units (vegetation, bare soil, urban, lake, and glacier) to compute the surface mass balance of glaciers in several fixed elevation ranges. In the simulation analyzed in this paper, 10 elevation classes have been used, with boundaries at 0, 200, 400, 700, 1000, 1300, 1600, 2000, 2500, 3000, and 10 000 m. The SMB at CISM grid resolution is calculated from the CLM SMB at these elevation classes by horizontal bilinear and vertical linear interpolations. Further details are given by Lipscomb et al. (2013).
Over Greenland, we replaced the CLM default glacier coverage and surface elevation with data from Bamber et al. (2001). Figure 1 illustrates the glacier cover fraction for each grid cell. The total area of the Greenland Ice Sheet in the model is 1.719 × 10 km2. The surface mass balance calculation has been built upon the CLM snowpack mass balance scheme. This provides the advantages of model consistency and a more efficient coupling between atmosphere and land, but also means that a true firn model (with snow transitioning to ice over a depth of tens of meters) is not included. The calculation of the SMB of glaciated surfaces, including the Greenland Ice Sheet, is done over two surface types, snow and bare ice, according to
where SNOW is snowfall, RAIN is rainfall, RU is runoff, and SU is sublimation.
Runoff is calculated as the portion of the available liquid water (ALW) that is not refrozen:
where ME is melt and RF is refreezing. ME is calculated from the sum of surface fluxes (radiative, turbulent, and subsurface heat):
where SWd is downwelling shortwave radiation, is albedo, is the emissivity of the snow/ice surface, σ is the Stefan–Boltzman constant, Ts is surface temperature, LWd is downward longwave radiation, SHF is sensible heat flux, LHF is latent heat flux, and GF is the subsurface heat flux toward the surface.
Retention of liquid water by capillary forces is considered by means of a parameter (irreducible water saturation) (Anderson 1976) in the equations describing the flow of water between snow layers (Oleson et al. 2010). The snowpack can reach a prescribed maximum thickness Hmax over glacier surfaces. Here, Hmax has been set to 1-m liquid water equivalent (LWE). If Hmax is reached, all additional snowfall is considered ice gain. If the snow depth H becomes zero, then additional melt is considered as ice loss. Melted ice does not refreeze. Percolation and refreezing of meltwater within the snowpack can take place only if 0 < H < Hmax. Rainfall falling at a point where H = Hmax is not allowed to percolate but runs off.
The maximum number of snow layers is five, and their thickness evolves during the simulations. The calculation of snow albedo is done using the SNICAR model (Flanner and Zender 2006) and takes into account snow grain size, solar angle, and aerosol content. Bare ice albedo is prescribed with two different values for visible and near-infrared frequencies. These values have been set to 0.60 and 0.40 [see Lipscomb et al. (2013) for details] for the present study.
The atmospheric forcing used in the SMB calculation is provided by CAM4 every half hour. From the atmospheric forcing, elevation corrections are applied to near-surface temperature (with an assumed lapse rate of −6 K km−1) and specific humidity (holding relative humidity constant). Snowfall and rainfall rates, as well as incoming radiation, are not corrected for elevation. Further details are given in Oleson et al. (2010) and Lipscomb et al. (2013).
b. Spinup and simulations setup
Two simulations of Greenland climate and surface mass balance (SMB) are analyzed: a preindustrial simulation with a length of 100 years and a historical transient simulation covering the period 1850–2005. Both have been performed with the fully coupled model CESM1–CISM including interactive calculations of the carbon and nitrogen cycles in the land component and of glacier surface mass balance in different elevation classes.
To obtain initial conditions for the preindustrial simulation, it was necessary to spin up the carbon cycle and the ocean, as well as the snowpack over glaciated surfaces at multiple elevation columns. For this purpose, the initial conditions for all model components except CLM were taken from an existing preindustrial simulation with CESM (without the new land ice component and the new land units with multiple elevation classes over glaciated surfaces). To initialize with a spunup snowpack, the initial conditions coming from this coupled preindustrial run were blended with another simulation that was performed with CLM including glacier_mec land units but forced offline with atmospheric data from the same coupled CESM simulation. This offline CLM simulation had a length of ~120 years, sufficient for spinning up the snowpack. It was started with a snowpack of thickness H = Hmax/2 = 0.5 m over glacier_mec land units, so as to minimize the time needed for spinup, and a snow temperature of −10°C in all layers.
All of the simulations analyzed in this study used an updated dataset for topography and ice sheet coverage over Greenland (Bamber et al. 2001). For this reason, the glacier coverage used to generate the atmospheric forcing for the offline glacier_mec simulation with GLAC_MEC is inconsistent with the glacier coverage in the preindustrial run, and the first half of this simulation was discarded prior to for analysis since the model is still equilibrating.
The historical simulation is branched from year 100 of the preindustrial simulation. The transient external forcing consists of insolation, volcanic eruptions, aerosols, and long-lived greenhouse gas concentrations.
c. Data for model evaluation
For evaluation of the model, the period between years 1960 and 2005 is used. The model output is evaluated using results from a regional model, RACMO2, which has been thoroughly evaluated against available in situ observations (Ettema et al. 2009, 2010a), GRACE data (van den Broeke et al. 2009), and other satellite measurements (Fettweis et al. 2011).
RACMO2 and preceding model versions have been applied to several studies of the Greenland Ice Sheet (Ettema et al. 2009; van den Broeke et al. 2009). The version used in this study has an improved simulation of surface albedo (van Angelen et al. 2012). RACMO2 is forced at the boundaries with data from the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) (Uppala et al. 2005) up to 1979 and with the ECMWF Interim Re-Analysis (ERA-Interim) (Simmons et al. 2007) thereafter. The model is run with a horizontal resolution close to 11 km, which is the highest among available simulations from regional models for Greenland that include an interactive snow model.
The climate fields and terms of surface mass balance over glacier land units that are plotted and analyzed in this paper have been averaged over different elevation classes, weighted by the fractional area at each elevation computed from a high-resolution topographic dataset. Unless otherwise specified, these fields are plotted only over grid cells belonging to the ice sheet with glacier coverage of at least 1%.
Note that RACMO2 has been forced by reanalysis data, which assimilates data from observations. This limits the natural variability of the large-scale climate to a state that might be quite different from the single CESM realization analyzed in this paper.
3. Simulated near-surface climate of Greenland
We first focus on near-surface (2 m) atmospheric temperature and the surface energy balance. The latter determines melt rates and surface temperature. Special emphasis is given to the summer energy fluxes since they are most important for ice sheet ablation. The surface albedo is mentioned here, but it is analyzed in more detail in the following section because of its strong ties with melt, snowfall, and rainfall. The Greenland climate and energy fluxes as simulated by RACMO2 are discussed and evaluated in detail by Ettema et al. (2010b).
Analysis of temperature and surface fluxes is done with the data provided by CESM for each land grid cell. The mean values represent the area-weighted average of all surface types in each grid cell, including vegetated as well as glaciated areas. Therefore, these mean values contain information from areas outside the ice sheet. The surface mass balance scheme uses values for the incoming radiation that are shared by all surface types, but upward radiation and net turbulent fluxes are calculated specifically for each surface type. In the analysis of the correlation between CESM and RACMO2 fluxes, two cases will be considered in order to evaluate the influence of the model biases over noniced surfaces. For the first case, grid cells with partial ice sheet coverage (defined as higher than 20%) will be considered; for the second, we include only grid cells with total ice sheet coverage (defined as higher than 99%).
The area of CESM grid points over the ice sheet decreases with increasing latitude, ranging between ~3000 and ~6000 km2. The number of RACMO2 grid points for the Greenland Ice Sheet (GIS) is 13 962. The number of CESM (land model) grid points with more than 50% glacier cover is 431. The highest RACMO2 grid point is at 3227 m and at 3142 m in CESM. The median heights are 2098 and 1928 m, respectively, and the interquartile ranges are 998 (RACMO2) and 1111 m (CESM).
a. Near-surface temperature
The simulated mean winter, December–February (DJF), near-surface temperature over the GIS is −31.2°C, only 0.4°C higher than for RACMO2. The pattern agrees well with RACMO2 (the correlation coefficient is r = 0.95). CESM simulates higher temperatures over the northern interior of the ice sheet (Fig. 2), very likely because of the smoothed topography. Near-surface temperatures over the ocean around Greenland are lower in CESM. At the margins of the ice sheet and in the southern interior, temperatures are also lower. The RMSE is 2.7°C for grid cells with more than 20% of ice sheet coverage and 2.5°C for grid cells with total ice sheet coverage. The colder climate on the eastern side of Greenland compared to the western side at the same altitude and latitude (Ettema et al. 2010b) does not have a corresponding pattern in CESM.
The simulated mean summer, June–August (JJA), near-surface temperature over the GIS is −7.1°C, which is 0.2°C higher than for RACMO2. The correlation with RACMO2 data is higher than for winter (r = 0.98) and the RMSE is smaller (1.3°C). The spatial standard deviation (STD) is also smaller during summer (2.9° and 6.1°C). Both STDs are smaller than in RACMO2 (summer and winter values of 3.5° and 7.7°C). The mean summer bias is 0.65°C. The RMSE is only slightly reduced when considering grid cells with more than 99% of ice sheet coverage (1.1°C). As in winter, temperatures are higher in CESM over the northern interior and lower over the rest of the ice sheet. Ocean temperatures around Greenland are also lower in CESM. Over the tundra, there are major discrepancies between models in the north where CESM temperatures range between −4° and 0°C, while RACMO2 simulates temperatures in the range of 2°–4°C. These lower temperatures in CESM can be at least in part attributed to the presence of sea ice at these northern latitudes in CESM (as evidenced by negative temperatures over the ocean). Also, the large glacier fraction within the grid cells containing tundra surfaces (see Fig. 1) causes major cooling of the cell-averaged climate (e.g., incoming longwave radiation). In the western tundra these glacier fractions are smaller, and the temperature bias is less.
b. Surface fluxes
The simulated mean summer downwelling shortwave radiation SWd for years 1960–2005 is 268 W m−2, which is 11 W m−2 lower than for RACMO2 (Table 1). The SWd shows a south-to-north gradient (Fig. 3, top center) and higher values (>250 W m−2) over the ice sheet than in the ice-free areas, caused by clouds. The eastern part of the ice sheet receives more radiation than the western counterpart, also due to different cloud cover. Several minima of SWd are located in the northwest (NW) and in parts of the western ablation area, as in RACMO2 (Fig. 3, middle center). The correlation with RACMO2 is much higher when areas with partial ice sheet coverage are considered (0.77 versus 0.54 for areas totally ice covered). CESM simulates lower SWd over most of the ice sheet and the areas with partial tundra coverage in the west and higher SWd over the rest of the partially ice-covered regions. The RMSE is 16 W m−2 (17 W m−2 if only totally ice-covered regions are considered).
The mean summer albedo is 0.77, lower than RACMO2 by 0.008 (Table 1). The correlation coefficient is r = 0.78 (r = 0.70 for cells totally ice covered). The RMSE is 0.10. Over most of the ice sheet, CESM simulates lower albedos. This is in part explained through the high dependence of albedo on snow grain size and, thus, on snow temperature. In the previous subsection, it was shown that the near-surface temperatures in the interior of the ice sheet were higher in CESM. Another process that can explain the lower albedo in CESM is rainfall events. CESM simulates rainfall events in the interior of the ice sheet during summer, a process that is unrealistic because of the low atmospheric temperatures there. For most of the areas where RACMO2 albedos are less than 0.7, CESM albedos are higher. This is likely connected with the lower near-surface temperatures simulated in ice-free regions of Greenland (see previous subsection), which favor a more persistent snow cover during summer.
The mean absorbed (net) solar radiation SWnet is 61 W m−2, 3 W m−2 lower than in RACMO2 (Table 1). The SWnet is minimum in the northern interior of the ice sheet, where the albedo values are high (~0.825) and increases toward the margins, with maximum values ~100 W m−2 in the western ablation area where bare ice is exposed during long periods in the summer. The correlation (r = 0.77) is similar to that of the albedo and SWd and decreases (r = 0.73) if only totally ice-covered grid cells are considered. The biases are smaller in the interior of the ice sheet (mean bias and RMSE are −1.3 and 7 W m−2 for the areas totally ice covered and −5.3 and 20 W m−2 when considering areas with at least 20% ice cover). The biases are mostly negative for areas with partial ice sheet coverage, regardless of the positive SWd biases over most of the margins.
The incoming (downward) longwave radiation LWd varies over Greenland mainly owing to differences in height (higher elevations receive less radiation because of a thinner overlying atmosphere) and cloudiness, which is also often related to the former. The LWd is also a function of temperature and humidity of the lowest atmosphere layers and especially so of surface elevation. The pattern of LWd is similar to SWd but reversed (Figs. 4, top left, and 3, top center). Maximum LWd is found over the tundra regions in the west and south. Also, as in RACMO2, there is a zonal asymmetry for the northern half of the ice sheet with the minimum LWd over the high eastern part. Over most of the ice sheet, CESM simulates higher LWd than in RACMO2, with the largest differences (>20 W m−2) in the northern high elevation areas where CESM does not simulate values below 216 W m−2. The mean LWd (235 W m−2) is 7 W m−2 lower than in RACMO2. The correlation between models is much higher than for SWd (r = 0.92), and the RMSE is similar (14 W m−2). Most of the grid cells with partial ice coverage have lower LWd than in RACMO2.
Over the areas in the north and east that do not belong to the ice sheet, LWd is slightly lower in CESM, probably for the same reason that SWd is higher, that is, a relatively high topography that reduces cloud coverage and, additionally, atmospheric thickness. Also, the higher sea ice cover and lower near-surface temperature in CESM described in the previous subsection influence LWd via lower humidity and temperatures of the lowest atmosphere layers. The link of this relationship between LWd and SWd with height will be further explored in future studies, in particular to evaluate the need of a height correction for radiation fluxes for the downscaling to the different elevation classes. As explained in the methods section, radiation is currently kept constant among all elevation classes for CLM glacier units.
The net longwave radiation LWnet is the difference between the downwelling radiation LWd and the emitted radiation, which depends on surface temperatures. The LWd varies between −70 and −30 W m−2 in summer over the ice sheet with the lowest values in the northeast (NE), as in RACMO2 (Fig. 4, middle). In general, the heat loss is higher toward the margins of the ice sheet. Over the western ablation area, however, the heat loss is reduced at the lowest elevations owing to more frequent melting there during summer (van den Broeke et al. 2008a). The mean LWnet is −46 W m−2, 7 W m−2 higher than in RACMO2 (Table 1). The correlation is worse than for LWd (r = 0.51) but it improves dramatically when considering only fully glaciated grid cells (r = 0.71). The anomaly pattern is similar to that of LWd. Over the northern high areas in the ice sheet, CESM generates LWnet that are >10 W m−2 larger than in RACMO2, due to higher LWd as described above. Over the northern tundra regions, CESM does not simulate the low values simulated in RACMO2, due to the negative temperature bias.
The mean net radiation Rnet (i.e., the sum of the net shortwave and longwave surface radiation) is 15 W m−2 (Table 1). Lower incoming and net shortwave fluxes over the ice sheet than in RACMO2 added to higher longwave fluxes result in mean net radiation of 15 W m−2, which is 4 W m−2 higher than in RACMO. The correlation is high (r = 0.79). While the annual mean radiation is negative over the ice sheet (not shown), the summer SWnet exceeds LWnet over most of the ice sheet, except in a small area in the northern interior with small negative fluxes (Fig. 5, top left). Over the margins, it reaches values close to 70 W m−2. The values of Rnet agree well with RACMO2 over the ice sheet margins, but they are higher in the interior of the ice sheet, owing to higher LWd and LWnet fluxes.
The turbulent sensible (SHF) and latent (LHF) heat fluxes are shown in Fig. 5 (top center left and right panels). The mean SHF (6.6 W m−2) is 0.2 W m−2 higher than in RACMO2. SHF takes negative values only at the highest elevations over the ice sheet and over most of the ice-free regions. The highest positive values (10–24 W m−2) are found between 1000 and 2000 m. In RACMO2, these values are also found in the lower part of the ablation areas (Fig. 5, middle left). The discrepancy is due, in part, to partial tundra coverage of the corresponding grid cells in CESM (Fig. 1). The correlation with RACMO2 is low (r = 0.21) but improves dramatically when grid cells with partial tundra coverage are excluded (r = 0.76). In northern Greenland outside the ice sheet, SHF values are in the same range as in the interior of the ice sheet (0–10 W m−2), corresponding to the values that RACMO2 simulates over glaciers and ice caps in this region. The areas covered by tundra in the north reach very negative values in RACMO2 (~100 W m−2). CESM does not model comparable negative values in the northern tundra, as it does in the western part of the ice sheet and to some extent in the east. This bias will be evaluated together with other biases in surface fluxes in this area of Greenland that have been described above.
The mean summer LHF (−8.3 W m−2) is 3 W m−2 lower than in RACMO2. LHF is negative all over the ice sheet in summer, indicating sublimation, with the highest absolute values (10–25 W m−2) taking place below 2000 m. Over the western ablation area, sublimation decreases slightly toward the margin. This is because sublimation is higher where the snow surface has higher temperatures but is not melting. This pattern is more obvious in RACMO2 owing to increased resolution. The correlation with RACMO2 is moderate (r = 0.61), and is the worst of individual fluxes when only areas exclusively ice sheet covered are considered (r = 0.65). Outside the ice sheet, absolute LHF values are higher than over the ice sheet. In the west, they are higher (25–50 W m−2) than in the east (10–25 W m−2), with similar values as in RACMO2. In the north, LHF is less negative than in RACMO2.
The modeling of the subsurface heat flux GF in CESM is conditioned by the limited thickness of the snow cover, which does not exceed Hmax = 1-m LWE (section 2). CESM correctly simulates annual mean fluxes that are close to zero in the accumulation zone and represent a net gain for the surface in percolation and ablation areas (not shown) owing to the heat released by refreezing of meltwater in the inner snow layers. The annual mean GF is 1.5 ± 2.2 W m−2, only 0.2 W m−2 higher than in RACMO. The correlation between the models is high when only fully glaciated CESM cells are considered (r = 0.85 versus r = 0.45, if partially ice-covered cells are included). Summer GF is less consistent between the two models. Over the real ice sheet, summer GF is negative (heat is subtracted from the surface) in the accumulation zone since the surface is warmer than the inner snow layers. In CESM, GF is also negative over the high south, while positive in RACMO2 (Fig. 5). Over other high areas, CESM shows higher fluxes owing to the shallower snow cover. In the percolation area, where GF is positive because of the heat released by refreezing, CESM shows lower values, in connection with the lower refreezing capacity of its snow cover. The mean summer GF is −1.4 W m−2 (Table 1), −1.2 W m−2 lower than in RACMO2. The correlation between both models is acceptable (r = 0.62 if only fully glaciated cells are considered).
In summary, summer surface heat fluxes over the Greenland Ice Sheet agree well between CESM and RACMO2, with the major difference being higher net radiative flux in the interior in CESM mainly caused by higher incoming longwave radiation and, in the southern part, also lower surface albedo. Lower incoming solar radiation in CESM only partially offsets this difference.
4. Simulated 1960–2005 surface mass balance
The mean surface mass balance for the years 1960–2005 integrated over the whole ice sheet amounts to 359 Gt yr−1 (Table 2), with a temporal standard deviation of 125 Gt yr−1. Both values agree well with RACMO2 values for the same period (376 ± 117 Gt yr−1), but this is mainly attributed to compensation of higher precipitation (+143 Gt yr−1) with higher runoff (+151 Gt yr−1) and sublimation (+14 Gt yr−1) in CESM.
Figure 6 compares CESM simulated 1960–2005 SMB at 5-km resolution with data from 475 locations over the GIS (Cogley 2004; Bales et al. 2009; van de Wal et al. 2012; other sources referenced in supplementary material of Ettema et al. 2009). Most of the observational data corresponds to the accumulation zone, with measurement periods varying from location to location. The modeled SMB agrees reasonably well with the observations (r = 0.80). Accumulation rates are overestimated over the northern interior. CESM reproduces observed accumulation rates in the southern part, except in the high-accumulation area of the southeast (SE) margin, where it fails to reproduce the highest SMB values from this compilation that are close to 4000 kg m−2 yr−1 (data from ice core drilled by Jason Box). Figure 6c shows the comparison with 1990–2007 data from eight stations along the K transect (van de Wal et al. 2012) in the western ablation and percolation area around 67°N. The modeled equilibrium line altitude is close to ~1500 m, as in the observations. The SMB differences are relatively small over 1000 m. The SMB gradient is underestimated in CESM, which models too high SMB over the lowest three stations. Two main factors can explain this mismatch. First, these stations are located over a narrow fjord that is framed by tundra. The model resolution is too coarse to reproduce this small-scale feature. Second, these stations are located within the “dark zone” (Wientjes and Oerlemans 2010), which shows values for bare ice albedo that are among the lowest of the ice sheet [between 0.3 and 0.4; see, for instance, Fig. 1 in van Angelen et al. (2012), including MODIS data]. The fixed value of bare ice albedo in CESM (~0.5 when averaged over the visible and near-infrared spectra) is too high for this area.
Figure 7 illustrates the map of mean SMB for 1960–2005 from CLM (~1° resolution) and CISM (5-km resolution) outputs. The output from CLM for each grid cell represents area-weighted averages over the elevation classes. The ablation area (i.e., the area with negative SMB) represents 10% of the ice sheet in CESM and 12% in RACMO. The ablation areas are similarly distributed in RACMO2 and CESM: the main ablation areas along the western, northern, and northeastern margins have a similar width, while the southeastern margin lacks an ablation area in both models. Also, the small dimension of the ablation areas of the midwestern and northwestern margins, which present high topographic gradients and local maxima for precipitation rates, are similar in both models. Main discrepancies occur in areas of complex topography, as in the mideastern margin, where CESM does not simulate some of the ablation areas captured by RACMO2. The simulated equilibrium line is placed at ~1000 m in the north and northeast and at ~1500 m in the southern half of the western margin.
When CESM and RACMO2 data are both interpolated to the 5-km grid, a correlation of r = 0.79 is found (Fig. 7e; Table 2). The SMB range is substantially smaller in CESM (from −2337 to 1188 kg m−2 yr−1) than in RACMO (from −3143 to 4405 kg m−2 yr−1). When the high accumulation areas in RACMO, where the SMB > 1188 m−2 yr−1 (4% of the ice sheet), are excluded, the correlation improves substantially (r = 0.86), and the RMSE is reduced from 309 to 199 kg m−2 yr−1. CESM simulates a higher SMB over the areas where RACMO2 simulates the lowest rates (mean bias = +508 Gt yr−1 for the 2% of the ice sheet where SMBRACMO2 < 1000 kg m−2 yr−1).
In the following, precipitation and melt will be analyzed separately. Seasonality and albedo will be evaluated in the two following subsections.
The mean total precipitation over the reference period is 866 Gt yr−1 (STD = 88 Gt yr−1), which exceeds the range of values from regional models (Table 2), owing to the smoother topography that permits a higher inland penetration of moisture. The precipitation map (Fig. 8) shows two bands of high precipitation along the southeastern and western margins and a strong reduction of precipitation toward the northern interior of the ice sheet, where rates are between 100 and 200 kg m−2 yr−1. The model captures the local maxima in the northwest. Rates over 1000 kg m−2 yr−1 are modeled along the southeast and southwest margins, as in RACMO2, but miss the peak rates there. Along the northern margin, the model shows a local maximum of more than 500 kg m−2 yr−1 that is associated with orography but misses other maxima resolved by RACMO2. The overall pattern agrees well with RACMO2 (r = 0.91), with major biases being overestimation of rates in the interior of the ice sheet and underestimation in the southeast, attributed to the smoothed topography of CESM, which prevents the model from fully capturing the processes associated with orographic forcing of precipitation. Also, CESM overestimates precipitation in the west in most of the area between 1000-m and 2000-m height, with rates of more than 700 kg m−2 yr−1 at locations where RACMO2 simulates precipitation of 300 to 500 kg m−2 yr−1. The RMSE is 274 kg m−2 yr−1.
In CESM, rainfall represents a very large fraction (15%) of the precipitation over the ice sheet with a total of 127 ± 19 Gt yr−1. Rainfall is also simulated at high elevations of the ice sheet (not shown), which is unrealistic because the year-round low temperatures in this region. This overestimation of precipitation is due to a model bias in the atmosphere model, which simulates rainfall events at temperatures well below freezing. This model bias has been corrected in the next version of the atmospheric model [the Community Atmosphere Model version 5 (CAM5)].
b. Melt and runoff
The total melt rate amounts to 568 ± 112 Gt yr−1, which is in the high range of previous model studies (Table 2). Melt rates of more than 1 kg m−2 yr−1 occur all over the ice sheet except for the northern half of the area with elevations higher than 2000 m (Fig. 8). Below 1000 m the highest melt rates are simulated along the western and southeastern margins. This is also the case for the heights between 1000 and 2000 m, with melt rates of 700–1000 kg m−2 yr−1. The lowest melt rates between these heights occur in the north, as in RACMO2 (Fig. 8). The correlation coefficient is 0.73, and the RMSE is 489 kg m−2 yr−1 (Table 2).
CESM simulates runoff all over the ice sheet (Fig. 8). In reality, the low melt rates above 2000 m and the low temperatures of the snowpack should cause all the available liquid water (defined here as the sum of meltwater and rainfall) to refreeze. However, as explained above, in CESM rainfall over snow surfaces that reach the maximum thickness Hmax is prescribed to run off. This maximum thickness is reached during most of the time in the accumulation area. This explains the differences between the runoff maps of CESM and RACMO2 (Fig. 8) in this area. If this area is excluded, the mean bias and RMSE are +414 and 675 kg m−2 yr−1.
Sublimation/evaporation amounts to 54 ± 3 Gt yr−1, in the midrange of the values from the regional models shown in Table 2, and very close to RACMO2 values (40 ± 4 Gt yr−1). Net deposition takes place in the northernmost part over 2000 m (negative values, Fig. 8). Over most of the other high areas, sublimation rates are less than 20 Gt yr−1. Sublimation/evaporation increases toward the margins and decreases again over part of the ablation zones. The highest values are reached in the south and southwest margins. The agreement with RACMO2 is high (r = 0.81, RMSE = 17 kg m−2 yr−2).
d. Seasonality and comparison with GRACE data
Figure 9 (top) illustrates mean daily values of melt, refreezing, precipitation, and sublimation averaged over the ice sheet for the reference period (1960–2005). Snow starts melting in mid-April, reaching maximum melt rates in July and sharply decreasing during August. There is little melt in September and none after October. All snowmelt refreezes during the first month of melt. Maximum refreezing rates take place in the months of higher snowmelt, representing a bit more than one-third of the snowmelt, at the peak of the melt season. Some refreezing occurs in the fall, when the only source of available liquid water is rainfall. Bare ice starts melting in early June, with maximum rates in late July. Approximately 90% of the melt of bare ice takes place in the two months July and August. From September, melt rates of snow and bare ice are similar.
Snowfall is relatively constant throughout the year, with a maximum in late summer and early fall. Rainfall occurs between May and November, with most of it occurring in July and August. Sublimation takes place mostly between March and the end of September, with maximum daily rates of 180 Gt yr−1 in early July.
Figure 9 (lower) compares CESM SMB with data of the GIS mass from GRACE. We use the release 05 data from the Center for Space Research spanning years from 2003 to 2011, with monthly frequency. We applied all common corrections to the data, as in van den Broeke et al. (2009), including a correction for the displacement of the earth's center of mass (Swenson et al. 2008). The cumulative mass anomalies for the Greenland Ice Sheet are derived from the GRACE data using the basin method of Wouters et al. (2008). A year-by-year comparison with the CESM is not possible because of internal variability of the model data. For this reason, we compare the seasonal cycle from nine CESM years (1996–2004). It must be noted that the time interval for GRACE is shifted forward with respect to the years from CESM; therefore, seasonal variations could be more influenced by the global warming trend. GRACE data represents cumulative mass anomalies, with monthly frequency. The data has been detrended (quadratic fit). CESM cumulative (surface) mass anomalies have been constructed from monthly SMB fields and detrended (quadratic fit). For this comparison, we assume that ice flux to the ocean is approximately constant throughout the year, as suggested in van den Broeke et al. (2009). For both model and satellite data, monthly anomalies to the annual mean have been constructed for each year. Nine-year-mean anomalies (thick lines) and the maximum and minimum anomalies for each month (thin lines) have been plotted in Fig. 9. The mean cycles are well correlated (r = 0.88). The CESM mean cycle (thick blue line) shows its maximum in May, as in GRACE (thick black line), and its minimum in August, a month earlier than GRACE. The amplitude of this cycle is similar in the model and data, with a maximum to minimum anomaly of 278 Gt in CESM (305 Gt in GRACE). Interannual variability in the seasonal cycle is greater in CESM. The modeled mass recovery after the late-summer minimum is somehow stronger in CESM for the mean values, although within the variability range shown by GRACE data.
Figure 10 illustrates the spatial variability in timing and duration of total melt and bare ice melt as modeled by the CESM. Melt starts in early spring in the areas of southern Greenland below 1000 m and in late April between 1000 and 2000 m. In the northern half of the ice sheet, melt starts in late spring, with the latest inception for the height range between 1000 and 2000 m taking place in the northwest. Areas between 1000 and 2000 m begin to melt in mid to late July, roughly at the same time as the areas over 2000 m in the southern third of the ice sheet. In the other two thirds of the high elevations of the ice sheet melt does not exceed the daily mean threshold of 1 kg m−2 yr−1.
The melt season lasts five months in the lower part of the western ablation area and more than six months at some smaller areas in the low south. In eastern Greenland and at elevations between 1000 and 2000 m of the western ablation area, the melt season lasts almost four months. In the north the season lasts between two and three months, depending on the location. The locations over 2000 m show a clear south-to-north gradient, with the melt season ranging from almost three months in the southernmost areas to only one day in the northernmost grid point that experiences melt according to the chosen daily mean threshold of 1 kg m−2 yr−1.
Figure 11 shows the surface albedo at the beginning (April) and peak (July) of the melt season only for the glacier fraction of the grid cells that belong to the ice sheet. These albedo values represent the weighted average over different elevation classes.
The CESM mean albedo over the ice sheet changes from 0.84 in April to 0.76 in July. This is higher than RACMO2 by 0.02 in both months. At the beginning of the melt season in April, the surface albedo over the ice sheet is above 0.825 at the locations where no melt has happened yet and slightly lower at the locations where there has been some melt, in the western and eastern margins below 2000 m (Fig. 10a). RACMO2 shows a uniform pattern of albedo exceeding 0.825 (not shown). At the peak of the melt season, albedo has dropped most in the western ablation area, reaching values close to 0.5, which are typical for bare ice (Fig. 10b). With the exception of the northern margin, most of the areas below 1000 m show albedo values below 0.7. Above 2000 m the albedo is that of dry snow (>0.825) in the northern half and decreases southward until reaching values close to 0.7, which are typical for wet snow in the southernmost areas.
5. SMB variability and trends
The mean surface mass balance for the period 1850–2005 is 386 ± 107 Gt yr−1. The preindustrial value, corresponding to years 55 to 100 of the preindustrial simulation, is similar for the mean with a much smaller standard deviation: 390 ± 66 Gt yr−1. The map of the preindustrial SMB (Fig. 12) is similar to the map for 1960–2005 (Fig. 6) but with smaller ablation areas (8% of the ice sheet surface).
Figure 13a displays the time evolution of the SMB, precipitation, melt, and runoff during 1850–2005. The net surface mass balance shows a negative trend of −0.75 Gt yr−2, significant at the 95% level. Both runoff and precipitation show positive trends: +0.87 and +0.14 Gt yr−2, respectively, with only the runoff trend being significant. The SMB ranges between 17 and 710 Gt yr−1. The second and third lowest SMB values are both reached after 1990. For all individual years in the period 1998–2005, the modeled SMB is below the 1850–2005 mean. The trend in summer near-surface temperature averaged over the whole ice sheet (Fig. 13b) is 0.0086°C yr−1, which is significant. The SMB over Greenland increases from 1991 to 1992 in response to the Mt Pinatubo eruption. Temperature, runoff, and melt decrease between these two years, while precipitation increases.
The standard deviations of the different integrated terms of the SMB compare well with RACMO2 (Table 1), indicating a similar variability in both models.
Figure 14 shows a map of SMB trends in CESM and RACMO. During 1960–2005 CESM has negative trends of more than 10 kg m−2 yr−1 in the lower part of most of the ablation areas except in the northeast, which has positive trends of ~3 kg m−2 yr−1 (Fig. 14b). These trends in the ablation areas are similar to those modeled by RACMO2 (Fig. 14c). RACMO2 does not model the positive trend in the northeast. Instead, positive trends are modeled at some areas along the western part of the ice sheet. This mismatch between models is most likely attributable to the differences in snowfall related to natural variability (RACMO2 is forced with reanalysis data). The 1960–2005 negative trend in the ablation areas is also present when considering the whole simulation period (Fig. 14a). The trend in the northeast has the opposite sign.
6. Summary and conclusions
We have described the performance of the Community Earth System Model (CESM) in the calculation of the surface mass balance (SMB) of the Greenland Ice Sheet. The CESM includes a completely new ice sheet component. The SMB is calculated in the land model at its coarse horizontal resolution with the same snow physics that is used in the model for nonglaciated surfaces. Over glaciated surfaces, the SMB is calculated at several prescribed elevations and then downscaled to a finer grid. In this study, the finer grid has a horizontal resolution of 5 km. This resolution is required for an adequate modeling of the SMB gradients at the margins of the Greenland Ice Sheet (van den Broeke et al. 2008a).
The simulated near-surface climate of the ice sheet was evaluated against output from the regional climate model RACMO2. The near-surface temperature over the ice sheet is in good agreement with RACMO2 but with higher temperatures over 2500 m in CESM. The radiation and turbulent fluxes from CESM compare well with RACMO over the ice sheet but show large biases over the tundra regions.
The surface albedo over the ice sheet changes from nearly uniform values before the beginning of the melt season to a very heterogeneous map at the peak of the season in July, with maximum values of more than 0.825 in the northern half over 2000 m and minimum values close to the albedo of bare ice (~0.5) in the lowest parts of the western ablation area. Intermediate values close to 0.7, which are typical of wet snow, are modeled at heights between 1000 and 2000 m. Since the surface energy balance is dominated by radiation fluxes (van den Broeke et al. 2008a,b), changes in surface albedo have a profound impact on the energy available for melt. The tight relationship between mass balance processes (snowfall, rainfall, snowmelt, refreezing, and bare ice melt) and albedo is well modeled in the CESM, enabling a correct coupling between SMB and local climate.
Melt rates are greatest over the western, northeast, and southeast ablation areas and decrease strongly with height, with melt rates below 1 kg m−2 yr−1 in the northern half above 2000 m. Runoff peaks at the same locations as melt. Over 2000 m, runoff is unrealistically large owing to the model constraint of limiting the snow thickness to a prescribed value and the assumption that additional rain runs off.
Rainfall is unrealistically simulated in the model to occur over the interior of the ice sheet. In addition, the total amount of rainfall is too large. This rainfall excess has an impact in several surface mass balance processes. It represents an artificial source of mass in regions where the snowpack has a thickness H between 0 m and Hmax of 1-m LWE, such as the percolation area and in the ablation area early in the melt season before bare ice is exposed. At these locations, rainfall reduces the refreezing capacity of the snowpack, which warms when rainfall freezes. Also, rainfall reduces surface albedo, causing transitions from the high albedo for dry snow to low albedo for wet snow.
The simulated melt season starts in April and lasts until September. Northern areas start melting in late May. The amount of refreezing is roughly equal to the amount of melt during the first month of the melt season. Total refreezing represents 35% of the available liquid water, which is the sum of rainfall and melt.
A major limitation in the modeling of refreezing in the CESM is the lack of a true snowpack model over the ice sheet. In the CLM snowpack model, the maximum thickness of snow cannot exceed a prescribed value Hmax. The thickness Hmax cannot be set to high values since, otherwise, the model cannot initialize with a correct distribution of ablation areas because negative SMB is only modeled at locations where the spring snow layer is completely melted and bare ice is exposed. While the approach in CESM of having a snowpack model in common with nonglaciated surfaces benefits from the advantages of model consistency and a more efficient coupling between atmosphere and land, a true firn model would provide a more adequate treatment of refreezing.
In a simulation covering 1850–2005, the simulated SMB variability agrees well with regional models. The SMB responds adequately to external forcing, for instance to volcanic eruptions such as the Pinatubo eruption in 1991. During 1960–2005, the SMB shows negative trends over the western and northern ablation areas, as well as in the northeast margin, in agreement with RACMO2.
After evaluating the processes included in the model and comparing to regional climate models, we conclude that the CESM is the first global climate model that is capable of realistically simulating the surface mass balance of the Greenland Ice Sheet, in spite of limitations imposed by the coarse horizontal resolution. The good match between CESM and state-of-the-art regional models and observations is not attributed to tuning (as in temperature-based melt index methods) but to adequate physical modeling of surface mass and energy balance processes, including albedo evolution, and an effective downscaling scheme. As shown in this study, a good simulation of the surface climate over Greenland is also strongly tied to a correct simulation of surface mass balance processes. Models with a weaker coupling between surface climate and surface mass balance are likely to suffer from biases in their simulation of the SMB and sensitivity to climate change.
M. Vizcaíno was funded via NSF Grant ATM-0917755 to UC Berkeley and a Marie Curie International Incoming Fellowship within the 7th European Community Framework (EC FP7) Programme (FP7-PEOPLE-2010-IIF-272956). W. Lipscomb was supported by the Scientific Discovery through Advanced Computing (SciDAC) project funded by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Biological and Environmental Research. The Los Alamos National Laboratory is operated by the DOE National Nuclear Security Administration under Contract DE-AC52-06NA25396. Additional support was provided by NSF through Awards ANT-1103686 (for W. Sacks) and a Small Grant for Exploratory Research (for W. Lipscomb). B. Wouters is funded by a Marie Curie International Outgoing Fellowship within the EC FP7 Programme (FP7-PEOPLE-2011-IOF-301260). The CESM project is supported by the National Science Foundation and the Office of Science (BER) of the U.S. Department of Energy. Computing resources were provided by the Climate Simulation Laboratory at the NCAR Computational and Information Systems Laboratory (CISL), sponsored by the National Science Foundation and other agencies. The CMIP5 simulations used in this research were enabled by CISL compute and storage resources. Bluefire, a 4064-processor IBM Power6 resource with a peak of 77 TeraFlops, provided more than 7.5 million computing hours, the GLADE high-speed disk resources provided 0.4 PetaBytes of dedicated disk and the CISL 12-PB HPSS archive provided over 1 PetaByte of storage in support of this research project. In addition, this work used resources of the Oak Ridge Leadership Computing Facility, located in the National Center for Computational Sciences at Oak Ridge National Laboratory, which is supported by the Office of Science of the Department of Energy under Contract DE-AC05-00OR22725. Janneke Ettema provided RACMO2 data from a previous version for preliminary evaluation of CESM. GRACE data were provided by Bert Wouters. J. Ettema, W. J. van de Berg, J. Lenaerts, and X. Fettweis are thanked for insightful discussions about the surface mass balance of the Greenland Ice Sheet.
Current affiliation: Department of Geoscience and Remote Sensing, Delft Institute of Technology, Delft, Netherlands.
The National Center for Atmospheric Research is sponsored by the National Science Foundation.
This article is included in the CESM1 Special Collection.