Journal topic
The Cryosphere, 13, 3193–3208, 2019
https://doi.org/10.5194/tc-13-3193-2019
The Cryosphere, 13, 3193–3208, 2019
https://doi.org/10.5194/tc-13-3193-2019

Research article 04 Dec 2019

Research article | 04 Dec 2019

# Surface mass balance downscaling through elevation classes in an Earth system model: application to the Greenland ice sheet

Surface mass balance downscaling through elevation classes in an Earth system model: application to the Greenland ice sheet
Raymond Sellevold1, Leonardus van Kampenhout2, Jan T. M. Lenaerts3, Brice Noël2, William H. Lipscomb4, and Miren Vizcaino1 Raymond Sellevold et al.
• 1Geoscience and Remote Sensing, Delft University of Technology, Delft, the Netherlands
• 2Institute for Marine and Atmospheric Research Utrecht, Utrecht University, Utrecht, the Netherlands
• 3Department of Atmospheric and Oceanic Sciences, University of Colorado Boulder, Boulder, CO, USA
• 4Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, CO, USA

Correspondence: Raymond Sellevold (r.sellevold-1@tudelft.nl)

Abstract

The modeling of ice sheets in Earth system models (ESMs) is an active area of research with applications to future sea level rise projections and paleoclimate studies. A major challenge for surface mass balance (SMB) modeling with ESMs arises from their coarse resolution. This paper evaluates the elevation class (EC) method as an SMB downscaling alternative to the dynamical downscaling of regional climate models. To this end, we compare EC-simulated elevation-dependent surface energy and mass balance gradients from the Community Earth System Model 1.0 (CESM1.0) with those from the regional climate model RACMO2.3. The EC implementation in CESM1.0 combines prognostic snow albedo, a multilayer snow model, and elevation corrections for two atmospheric forcing variables: temperature and humidity. Despite making no corrections for incoming radiation and precipitation, we find that the EC method in CESM1.0 yields similar SMB gradients to RACMO2.3, in part due to compensating biases in snowfall, surface melt, and refreezing gradients. We discuss the sensitivity of the results to the lapse rate used for the temperature correction. We also evaluate the impact of the EC method on the climate simulated by the ESM and find minor cooling over the Greenland ice sheet and Barents and Greenland seas, which compensates for a warm bias in the ESM due to topographic smoothing. Based on our diagnostic procedure to evaluate the EC method, we make several recommendations for future implementations.

1 Introduction

During the 20th century, the Arctic warmed much faster than the rest of the world due to shrinking sea ice cover , associated positive albedo–temperature feedbacks , and increased moisture and heat transport from the midlatitudes . The Greenland ice sheet (GrIS) lies within this fragile and rapidly changing environment. The GrIS is the world's second largest ice sheet, after the Antarctic ice sheet, and has an estimated volume of 2.96×106 km3 of ice, which would lead to an increase in global mean sea level by 7.36 m if it were all melted . Since the 1990s, the GrIS has lost mass at an accelerated rate . This mass loss is projected to be sustained and contribute a 0.04–0.21 m sea level rise by the end of the 21st century, depending on the climate scenario . The range of uncertainty is due to uncertainties in climate scenarios, climate sensitivity, and simulated mass balance of the GrIS by ice sheet models (ISMs). This latter uncertainty is currently being targeted by the Ice Sheet Model Intercomparison for CMIP6 (ISMIP6; Nowicki et al.2016), a major international effort to investigate future ice sheet evolution, constrain estimates of future global mean sea level, and explore ice sheet sensitivity to climate forcing.

State-of-the-art Earth system models (ESMs; coupled climate models capable of simulating the Earth’s chemical and biological processes, in addition to the physical processes, Flato2011) typically operate at a resolution of 1 (∼100 km), which poses a challenge for studies with a regional interest, such as GrIS surface mass balance (SMB). For instance, the extent of GrIS ablation areas may be underestimated . Also, there is a significant disparity between different model estimates of GrIS SMB even for models with higher resolution (Fettweis2018). Downscaling techniques are likely required to capture realistically the sharp gradients of SMB with elevation in the GrIS ablation zone . The most common downscaling techniques for the GrIS SMB are listed as follows.

1. Dynamical downscaling. This is done in regional climate models (RCMs, e.g., Box and Rinke2003; Noël et al.2018; Fettweis et al.2017) and recently as regional grid refinement within ESMs . This type of downscaling allows for explicit modeling at relatively high resolution for a region of interest. Physical parameterizations need to be readjusted over the fine grid , and in some cases, the model physics can be better tuned for this region. A major disadvantage of this downscaling method is the computational cost and the dependency on another model for lateral forcing in the case of RCMs.

2. Statistical downscaling . This uses elevation corrections on either SMB or components of SMB (e.g., runoff). This type of downscaling is successful when realistic topographic gradients of SMB or melt are captured in the model . However, in an ESM these gradients are typically not well-captured , making this technique unsuitable.

3. Hybrid downscaling. This is where elevation corrections are applied to components of SMB or surface energy balance (SEB), and the full SEB and/or SMB are explicitly calculated offline at a higher resolution. This method was used by to construct an SMB field from a global climate model for coupling to an ice sheet model.

A variant of the hybrid approach with “online” (that is, within the ESM) implementation has been developed recently. This method is based on the use of elevation classes (ECs) . It simulates the SEB and SMB over glaciated surfaces, with specific albedo and snowpack evolution for each EC. A benefit of this online approach is that it is able to capture feedbacks between the downscaled surface simulation and the atmospheric component of the ESM. This method has been successfully applied to the simulation of historical and RCP8.5-scenario projections of the GrIS SMB and mass balance evolution with the Community Earth System Model version 1.0 (CESM1.0). However, the EC downscaling in itself and its effects on the downscaled SMB and SEB components in CESM1.0 or other models have not been analyzed or evaluated before. Our study aims to fill this gap in three steps. First, we compare the simulated EC gradients of SMB and SEB components with gradients simulated by an RCM. Second, we investigate the sensitivity of the GrIS surface mass balance simulation to the main EC downscaling parameter, i.e., the temperature forcing lapse rate (it must be noted that our model does not downscale precipitation). Third, as the downscaling of SMB in the ECs takes place online within the climate model, we investigate how the EC implementation impacts regional climate.

Although we analyze the particular EC implementation in a specific ESM (CESM1.0), we aim to provide an evaluation and diagnostic framework to guide the future implementation of EC downscaling in other climate models, for offline SMB estimates and/or forcing of ice sheet models.

The paper is structured as follows: Sect. 2 describes the modeling setup as well as the regional model used for evaluation. In Sect. 3 we present the results. The discussion (Sect. 4) addresses the strengths and limitations of the EC implementation in CESM1.0. Section 5 gives the main conclusions and outlook.

2 Methods

## 2.1 CESM1.0 and EC downscaling scheme

The model used for this study is the Community Earth System Model 1.0.5 (CESM1.0) with all components active. The atmospheric model is the Community Atmosphere Model 4 (CAM4; Neale et al.2013), which is run at a horizontal resolution of 0.9× 1.25 and has a finite volume dynamical core. The land model is the Community Land Model 4.0 (CLM4.0; Lawrence et al.2011), which is run at the same horizontal resolution as CAM4. Within a CLM4.0 grid cell, different land cover types can exist. The grid cell average passed to the atmosphere is calculated with an area-weighted average of the fluxes. The ocean is simulated with the Parallel Ocean Program 2 (POP2; Smith et al.2010) with a nominal resolution of 1. The ocean model grid is a dipole with its northern pole centered over Greenland to prevent numerical instabilities, implying a higher effective resolution around Greenland. Sea ice is modeled with the Los Alamos Sea Ice Model 4 (CICE4; Hunke et al.2010; Jahn et al.2012), which runs on the same grid as the ocean. The ice sheet model in CESM1.0 is the Glimmer Community Ice Sheet Model 1.0 (CISM1.0; Rutt et al.2009; Lipscomb et al.2013), with a default resolution of 5 km. For the simulations performed in this study, the GrIS ice thickness and extent do not evolve. A static ice sheet surface that corresponds to present-day observations is used to downscale SMB, energy fluxes, and other quantities at the land–atmosphere interface through the EC scheme.

The steps for of the EC calculation in an ESM are roughly as follows.

1. A set of elevation classes are defined for each (partially) glaciated grid cell in the land model.

2. A selected set of atmospheric variables are downscaled by applying simple elevation corrections (typically, prescribed lapse rates).

3. The land model calculates the SEB and SMB per EC.

4. EC outputs are area-averaged per grid cell, and these averages are coupled to the atmospheric component.

In the following, the EC calculation is described in more detail. SMB calculations in CESM1.0 are done in CLM4.0 through ECs using the CLM4.0 snowpack mass balance scheme. EC downscaling accounts for subgrid elevation variability. SMB is explicitly calculated at multiple surface elevations to force the higher-resolution ice sheet model. The EC calculation is activated in the glaciated fraction of any grid cell with total or partial glacier coverage within a pre-defined region of interest (e.g., Greenland for the present study).

The EC method takes subgrid surface elevation from the ice sheet model and bins them into n ECs. In this study, n is 10 and the n+1 boundaries are fixed at 0, 200, 400, 700, 1000, 1300, 1600, 2000, 2500, 3000, and 10 000 m elevation a.s.l. The choice of n=10 was motivated by a compromise between computing time and increased (vertical) resolution. The offline test showed this number to be appropriate and is the default for CESM1.0. After this binning, CLM4.0 calculates the relative weight of each EC within a given grid cell, as well as the mean topography for each EC. The weight of each EC within a grid cell is determined by the area of the high-resolution topography dataset that lies within an EC. These weights are used to calculate the grid cell average that will be the output of CLM4.0 and coupled to CAM4, as well as for the interpolation of SMB and ice sheet surface temperature (which is equivalent to the temperature at the bottom snow/ice layer in CLM), which are standard forcings for ice sheet models .

Through the coupling with the atmosphere model, CLM4.0 receives surface incoming shortwave and longwave radiation, precipitation, 10 m wind, relative and specific humidity, surface pressure, and 2 m air temperature. Incoming radiation, precipitation, and wind are kept constant across all ECs within a grid cell. In contrast, the method downscales near-surface (2 m) air temperature to the ECs with a default lapse rate of 6 K km−1, and specific humidity is downscaled by assuming the relative humidity to be constant with elevation . At each EC, an energy balance model is used to calculate the surface energy balance every 30 min (SEB; W m−2) as

$\begin{array}{}\text{(1)}& M={\mathrm{SW}}_{\mathrm{in}}\left(\mathrm{1}-\mathit{\alpha }\right)+{\mathrm{LW}}_{\mathrm{in}}-\mathit{ϵ}\mathit{\sigma }{T}^{\mathrm{4}}+\mathrm{SHF}+\mathrm{LHF}+\mathrm{GHF},\end{array}$

where M is the melt energy (W m−2), SWin is the incoming solar radiation (W m−2), α is the surface albedo (–), LWin is incoming longwave radiation (W m−2), ϵ is surface emissivity (–), σ is the Stefan–Boltzmann constant (W m−2 K−4), T is the surface temperature (K), SHF is the sensible heat flux (W m−2), LHF is the latent heat flux (W m−2), and GHF is the subsurface heat flux into the snow or ice (W m−2). For these surface fluxes, positive values indicate energy transfer from the atmosphere to the land surface and from the subsurface to the surface for GHF. Snow albedo is calculated based on snow grain size, depth, density, and other properties . The first term on the right-hand side of Eq. (1) is the net solar radiation, and the sum of the second and third term on the right-hand side is the net longwave radiation. As a result of the SEB calculation, CLM4.0 calculates prognostic temperature, wind, relative humidity, and other quantities, taking into account the simulated exchanges of heat and moisture and surface roughness.

Additionally, the SMB (millimeter water equivalent per year, referred to as mm yr−1 in this paper) is calculated at each EC, with the same frequency as the SEB calculation, as

$\begin{array}{}\text{(2)}& \mathrm{SMB}=\mathrm{SNOW}+\mathrm{REFR}-\mathrm{MELT}-\mathrm{SUBL},\end{array}$

where SNOW refers to the snowfall rate, REFR is the refreezing rate of snowmelt and rainfall, MELT is the sum of snow and ice melt rates, and SUBL is the rate of sublimation/evaporation minus deposition/condensation. Rain and meltwater that do not refreeze are routed to runoff. For further details on the calculation of SEB and SMB, see . Total snow mass is limited to 1 m water equivalent.

The resulting SMB is linearly interpolated onto the ice sheet grid in two steps: first, with a bilinear horizontal interpolation per EC; and second with a vertical linear interpolation between two ECs (above and below), based on the ice sheet model high-resolution topography.

## 2.2 Simulation design

We perform four CESM1.0 simulations with an identical setup, except for a different temperature lapse rate forcing to each of the n=10 ECs. These four lapse rates are 1, 4, 6 (default), and 9.8 K km−1, and we refer to the corresponding simulations as EC-1K, EC-4K, EC-6K, and EC-9.8K, respectively. EC-1K is chosen to represent minimal activation of the EC calculation. EC-4K is chosen as a lapse rate forcing between EC-1K and EC-6K that is close to the summer lapse rate over the Greenland ice sheet as estimated from observations (e.g., Fausto et al.2009). As the upper limit of the magnitude of the lapse rate, 9.8 K km−1 (dry adiabatic lapse rate) is used.

All simulations start in 1955 from a CMIP5 historical run that is evaluated in detail in (which also describes the spinup procedure and the setup for the EC-6K) and run to 2005. All CESM1.0 model components are allowed to vary freely. The first 10 years are used for model adjustment to the new lapse rate, leaving the period 1965–2005 for analysis.

## 2.3 RACMO2.3 and evaluation procedure

For evaluation of the EC downscaled simulation of SEB and SMB, we compare with the dynamical downscaling in the Regional Atmospheric Climate Model version 2.3 (RACMO2.3; Noël et al.2015) with a horizontal resolution of ∼11 km and forced by the ERA-Interim reanalysis . We analyze the period between 1965 and 2005 for both RACMO2.3 and CESM1.0. As we are only comparing CESM1.0 simulations with identical initial conditions, we are likely to sample a different realization of climate variability than the reanalysis-forced RACMO2.3. RACMO2.3 has been successfully evaluated in multiple studies by comparison with in situ and remote sensing observations . Version 2.3 includes updates in cloud microphysics, surface and boundary layer microphysics, radiation, and precipitation . For the latter, precipitation falls exclusively as snow when near-surface temperatures are between 7 and −1C.

For the comparison, we use SEB and SMB components simulated at each EC with those simulated at the native grid of RACMO2.3. For CESM1.0, this results in between 1 and 10 values per CLM4.0 grid cell, depending on subgrid elevation heterogeneity. We subtract each EC value of SEB or SMB component from the grid cell average, as well as the corresponding EC topographic height from the CLM4.0 mean height. We subtract these averages to only capture gradients within each grid cell, and to reduce the effect of internal climate variability. With these differences, we calculate a gradient or linear function with elevation. To generate these gradients for RACMO2.3, we first cluster RACMO2.3 model output from the 11 km native grid onto the CLM4.0 grid (∼100 km). We then calculate averages for each RACMO2.3 SEB/SMB component and surface elevation over the coarse CLM4.0 grid cells. We subtract these averages from the native original values, and we construct the gradients via a linear fit. In this way, up to 56 RACMO2.3 grid cells are mapped into each CLM4.0 grid cell, giving a total of 13 311 points for evaluation. For CLM4.0, the resulting number of points is 1551.

For comparison of the overall downscaled SMB in CESM1.0 to a previous RACMO version (2.1), and an evaluation of the simulation at the mean elevation, see .

3 Results

## 3.1 Process-based comparison of EC and dynamical downscaling

We use CESM1.0 output from a simulation using the default lapse rate forcing of 6 K km−1 (EC-6K). Figure 1 illustrates the comparison of the downscaled SEB component gradients for the CESM1.0 ECs and RACMO2.3 RCM. Regression slopes, m (gradients) and r values (correlation with elevation), are given in Table 1.

Figure 1Comparison of EC downscaling (simulation EC-6K, blue) versus dynamical downscaling in a RCM (RACMO2.3, black) for several summer (JJA) SEB components and near-surface climate, (a) albedo, (b) net solar radiation (W m−2), (c) net longwave radiation (W m−2), (d) sensible heat flux (W m−2), (e) latent heat flux (W m−2), and (f) melt energy (W m−2). The x values show the deviation of surface elevation (m) from the coarse grid cell (∼100 km) mean, and the y values show the deviation of the physical quantity from the grid cell mean. In plots (b) through (f), positive y values indicate more energy available for melting. Melt energy (f) is the sum of the radiation and turbulent flux terms in (b) through (e) plus the ground heat flux (not shown). The lines represent least-squares linear regression. The annotated m is the least-squares linear regression gradient (W m−2 km−1 or km−1 for albedo); r is the correlation coefficient.

Table 1Gradients (m) and correlation with elevation (r; unitless) of surface energy and mass balance components as simulated through EC downscaling in CESM1.0. Values correspond to JJA (energy) and annual (mass) averages for the period 1965–2005. Melt energy is the sum of the net shortwave and longwave radiation and the heat fluxes. Surface mass balance is the sum of snowfall and refreezing, minus melt and sublimation.

Figure 2Same as Fig. 1 but for annual SMB components from EC-6K (blue) and RACMO2.3 (black). (a) Snowfall (mm yr−1), (b) surface melt (mm yr−1), (c) refreezing (mm yr−1), and (d) surface mass balance (mm yr−1). Surface mass balance is the sum of snowfall (a) and refreezing (c) minus the surface melt (b) and sublimation (not shown). The lines represent least-squares linear regressions. The annotated m is the least-squares linear regression gradient (mm yr−1 km−1); r is the correlation coefficient.

The downscaled net longwave radiation (difference between incoming and outgoing longwave radiation, Eq. 1) in CESM1.0 has an opposite elevation gradient (8.9 W m−2 km−1) compared to RACMO2.3 (−3.1 W m−2 km−1) as shown in Fig. 1c. That is, the net longwave energy available for melting increases with lower elevation for RACMO2.3 but decreases with lower elevation for CESM1.0. The reason for this difference is that CESM1.0 does not downscale the incoming longwave radiation, while RACMO2.3 simulates a gradient of −17.6 W m−2 km−1 with a relatively high correlation with elevation ($r=-\mathrm{0.81}$, Fig. S1c, Table 1). This negative correlation in RACMO2.3 is caused by thicker clouds as well as higher water vapor and atmospheric temperatures at lower elevations . As the outgoing thermal radiation depends on the surface temperature, both models simulate negative gradients (Fig. S1d). The result is a positive gradient for the net longwave in CESM1.0. In RACMO2.3, the magnitude of the outgoing longwave gradient is smaller than the incoming longwave gradient, resulting in a net negative gradient. Due to the complex relationship between the different components of the longwave radiation, the net longwave has a low correlation with elevation in RACMO2.3 ($r=-\mathrm{0.30}$). In contrast, CESM1.0 simulates a high correlation (r=0.76) as the surface temperature gradient directly controls the net longwave gradient. The net radiation gradient in CESM1.0 is 5.4 W m−2 km−1 and in RACMO2.3 it is −22.6 W m−2 km−1 (Table 1).

In summary, biases in the downscaling of net radiation in CESM1.0 are due to null gradients of incoming radiation in the model, as well as weaker albedo gradients. As a result, the gradient is dominated by the outgoing longwave gradient in CESM1.0 and by the albedo and incoming longwave gradients in RACMO2.3.

Next, turbulent fluxes of latent and sensible heat are examined, as well as their contribution to the available melt energy with respect to radiation. The gradients of sensible and latent heat fluxes are negative in both models (Table 1); more energy is available for melting at lower elevation. The sensible heat flux gradient is stronger than the latent heat flux gradient and shows a larger spread of values (Fig. 1d, e). In CESM1.0, this is a result of the elevation correction applied to the near-surface temperature (lapse rate). This correction increases atmospheric temperature and specific humidity at lower ECs and decreases them at higher ECs within each coarse grid cell. In RACMO2.3, these heat flux gradients are smaller and less correlated with elevation ($r=-\mathrm{0.42}$ and $r=-\mathrm{0.02}$, for sensible and latent heat fluxes, respectively) than in CESM1.0 ($r=-\mathrm{0.77}$ and $r=-\mathrm{0.76}$). Stronger sensible and latent heat gradients in CESM1.0 appear to compensate for most of the underestimation of the radiation gradients (Fig. 1c, d, e.), resulting in a melt energy gradient (−16.0 W m−2 km−1) which is similar in magnitude and sign to RACMO2.3 (−26.1 W m−2 km−1; Fig. 1f, Table 1).

Figure 2 compares snowfall, surface melt, refreezing, and SMB gradients between the two models. While CESM1.0 does not downscale snowfall, RACMO2.3 simulates an elevation gradient of −218 mm yr−1 km−1 that has little correlation with elevation (r=0.26), possibly due to the competition of the dominant effect of height desertification (less snowfall at higher elevations due to colder and drier air), orographic forcing of snowfall, and small-scale atmospheric circulation features . Consistent with the melt energy gradients, the surface melt gradient in RACMO2.3 is −717 mm yr−1 km−1 while for CESM1.0 it is −425 mm yr−1 km−1 (Table 1).

On the other hand, the CESM1.0 refreezing gradient (62 mm yr−1 km−1) is in disagreement with RACMO2.3 (−129 mm yr−1 km−1 and Fig. 2c). CESM1.0 simulates a positive gradient, implying increasing refreezing at higher ECs despite reduced melt rates. We hypothesize that at low-elevation ECs this is due to limited refreezing capacity in CLM4.0, as a result of the limited snow depth (Sect. 2.1). On the contrary, at the higher ECs, where the melt is lower, refreezing is favored due to lower snow temperatures, more available pore space, and thicker snowpacks. The overestimation of rainfall at higher elevation may also be an important factor. In contrast to CESM1.0, RACMO2.3 simulates a negative gradient of −129 mm yr−1 km−1 (Table 1), suggesting a dominant control from the increased melting at lower elevation. As the refreezing gradient results from the combination of opposite gradients, i.e., available meltwater and available refreezing capacity, the correlation with elevation is low in RACMO2.3 ($r=-\mathrm{0.45}$, Table 1). It is similarly low in CESM1.0, in part due to lower correlation for the melt gradient than in RACMO2.3.

Regardless of substantial differences in melt gradients in both models, the SMB gradient is relatively close (Fig. 2d; CESM1.0: 439 mm yr−1 km−1; RACMO2.3: 369 mm yr−1 km−1; Table 1). CESM1.0 compensates for underestimation of the melt gradient with the snowfall and refreezing gradients (in order of importance; see Table 1). In addition to the snowfall contribution of +218 mm yr−1 km−1 to the CESM1.0 SMB gradient difference with RACMO2.3, the difference in the refreezing gradient contributes with +191 mm yr−1 km−1. The higher-elevation correlation of SMB with elevation in CESM1.0 (r=0.58) compared to RACMO2.3 (r=0.27) is due to the null precipitation gradient in CESM1.0.

In summary, the EC method in CESM1.0 with the default lapse rate of 6 K km−1 (EC-6K) is approximately reproducing SMB gradients of RCM RACMO2.3. The EC method partially compensates for the biases in radiation downscaling with an overestimated turbulent heat flux gradient. The resulting melt energy gradients, however, are still lower than in RACMO2.3. However, the EC method compensates for this in the net SMB gradient due to lack of snowfall downscaling (leading to a more positive gradient relative to RACMO) and a positive bias in the refreezing gradient.

## 3.2 EC downscaling sensitivity to lapse rate of temperature forcing

Figure 3 shows how the most relevant energy fluxes and SMB respond to different lapse rate forcings. With a larger lapse rate forcing, the simulated sensible heat flux gradient is stronger, from −3.2 W m−2 km−1 in EC-1K to −20.0 W m−2 km−1 in EC-9.8K (Fig. 3a–d). This implies that the stronger the lapse rate forcing, the more heat is redistributed from upper to lower elevations. The correlation with elevation only increases marginally when increasing the lapse rate forcing (Fig. 3a–d).

Figure 3Comparison of 1965–2005 summer (JJA) downscaled energy fluxes among four simulations with different elevation corrections for the atmospheric temperature forcing. The first column corresponds to EC-1K, the second to EC-4K, the third to EC-6K, and the last to EC-9.8K for (a)(d) sensible heat flux (W m−2), (e)(h) albedo (-), (i)(l) melt energy (W m−2), and (m)(p) surface mass balance (mm yr−1). The x and y values represent the deviation of surface elevation and energy component for each data point with respect to the climate model grid (∼100 km) mean. The lines represent least-squares linear regressions. The annotated m is the least-squares linear regression gradient (W m−2 km−1 or km−1 for albedo); r is the correlation coefficient.

Albedo gradients are sensitive to lapse rate forcing, from close to zero gradients in EC-1K to 0.029 km−1 in EC-9.8K (Fig. 3e–h). Even with the maximum lapse rate forcing, CESM1.0 is only able to produce an albedo gradient that is 35 % of the RACMO2.3 gradient. Albedo gradients are triggered by surface temperature and melt gradients resulting from turbulent heat flux gradients. In the case of EC-1K, the turbulent heat flux gradient is not sufficient to trigger substantial albedo–melt feedback. Downscaled albedos have a variation range of similar magnitude in EC-4K and EC-6K; however, more points in EC-6K have non-null variations.

The combined effects of the turbulent heat flux gradients and the associated albedo gradients result in higher melt energy gradients with higher lapse rate forcing (Fig. 3i–l). The melt energy gradient in EC-1K is −3.5 W m−2 km−1, which is very similar to the sensible heat flux gradient (−3.2 W m−2). With higher lapse rate forcings, the difference between melt energy and sensible heat gradients becomes larger, which is interpreted as an effect of the albedo–melt feedback.

The melt energy gradient as simulated by RACMO2.3 is best matched with EC-9.8K (Fig. 3, Table 1). However, EC-6K matches the SMB gradient best (SMB gradients for EC-1K, EC-4K, and EC-9.8K are 110, 310, and 711 mm yr−1 km−1; Fig. 3; compare with Table 1). This is explained by compensation from the snowfall and refreezing gradients.

Figure 4 compares the downscaled SMB maps on the ice sheet model grid (5 km resolution) for the four lapse rates and RACMO2.3 (11 km resolution). Spatially, the largest responses to a varying lapse rate occur along the margin of the ice sheet and close to the equilibrium line (Fig. 4c, d, e). At the margins, a low lapse rate leads to a higher SMB with respect to EC-6K in a very narrow band of only 10–20 km, due to the aforementioned relatively low turbulent fluxes and weak albedo–temperature feedbacks. In the EC-9.8K, this effect becomes opposite, resulting in a similarly narrow band of lowered SMB (blue rim). Further inland, this extreme lapse rate leads to larger areas with higher SMB, as higher melt energy gradients reduce melt at high-elevation ECs.

Figure 4Climatological (1965–2005) SMB for RACMO2.3 (a), CESM1.0 downscaled to 5 km (b), and SMB anomalies (c, d, e) (mm yr−1) using lapse rates (c) 1 K km−1, (d) 4 K km−1, and (e) 9.8 K km−1. Anomalies are with respect to the default lapse rate of 6 K km−1. Solid black contour shows the ice sheet margin. Elevation contours (dashed) are plotted every 500 m. The black line shows the ice sheet margin. Black dots show where differences are significant at the 95 % level according to a Student t test.

Larger lapse rates result in reduced ablation area, from 16.4 % of the GrIS in EC-1K to 13.0 % in EC-9.8K (Table 2). This reduction is due to an enhanced melt gradient (Fig. 3i–l), reducing melt at higher ECs and resulting in a lower equilibrium line altitude (ELA, where SMB equals zero), and it reduces interannual variability (although only mildly, from 4.0 % to 3.0 %). Due to this expansion of the accumulation area with higher lapse rates, the total SMB of the accumulation area increases (Table 2), although within the standard deviations. For the SMB of the ablation area, the area reduction is partially compensated for with higher specific (local) ablation rates for higher lapse rates, resulting in the most negative SMB in the ablation area for EC-4K. The total SMB is the sum of the SMB for ablation and accumulation areas, and it is maximum for EC-6K. The SMB for EC-6K is at the same time the closest to RACMO2.3, also for the standard deviation. However, the range of variation of the mean total SMB across the four simulations is not large and is within the standard deviations. As an additional note of caution, the values in Table 2 result from four simulations with independent atmospheric simulation, perhaps sampling different segments of, e.g., multidecadal precipitation variability , and they therefore reflect more than just the effect of the lapse rate choice.

Table 2Simulated whole-ice-sheet SMB, ablation area, total SMB in the ablation and accumulation areas, and prognostic near-surface temperature gradients for the four simulations performed in this study with varying lapse rates and for the reference regional model RACMO2.3. Values correspond to the climatological (1965–2005) average with the standard deviation in parentheses.

To summarize, lapse rates lower than EC-6K result in larger ablation areas and lower integrated SMB. These results indicate a dominant effect on the CESM1.0 ELA simulation of higher melt rates at high-elevation ECs versus reduced melt rates at low-elevation ECs.

To complete this sensitivity investigation, we compare prognostic near-surface temperature gradients across the four simulations (Table 2). This prognostic temperature is calculated per EC within each CLM4.0 time step and is a result of heat and moisture exchange between the surface and atmosphere. Therefore it differs from the prescribed lapse rate forcing. The prognostic temperature gradients are lower in magnitude than the respective lapse rate forcing for all CESM1.0 simulations. The magnitude of the June–August (JJA) gradient is also less than for December–February (DJF) and is approximately half of the forcing lapse rate. The former is also the case for RACMO2.3. The simulation EC-9.8K gives the prognostic temperature gradient closest to RACMO2.3, which is in between the EC-6K and EC-9.8K gradients. It is remarkable that the simulation EC-4K with the lapse rate forcing that is closest to the observational summer gradient (4.7 K km−1, ), and RACMO2.3 (4.3 K km−1) is not the simulation with the closest prognostic gradient.

## 3.3 Impact of the EC calculation on regional climate simulation

Next, we examine how the EC calculation in the land component (CLM4.0) affects the simulation of Arctic climate in CESM1.0. If the EC method is active in CLM4.0, subgrid gradients in the ice sheet surface budget are coupled to the atmosphere model (and via the atmosphere to other components) during runtime. We compare two simulations for this analysis. The EC-1K simulation serves as the control as it represents the simulation closest to non-active EC downscaling, which is the standard for most CMIP5 ESMs. The EC-6K is used to assess the climatic effect of using the EC method. Figure 5 shows differences in selected climate variables between EC-6K and EC-1K.

Figure 5Annual climatology (1965–2005) of EC-1K (left column) and anomalies of EC-6K with respect to EC-1K, which approximate the EC imprint (right column). (a, b) Near-surface air temperature (K), (c, d) turbulent (sensible + latent) heat fluxes (W m−2), (e, f) net longwave radiation (W m−2), and (g, h) sea ice concentration (–). Black dots indicate significance at the 95 % level according to a Student t test. Positive signs for (a–f) indicate energy transfer from atmosphere to the surface.

Near-surface temperatures decrease over large parts of the GrIS and on average by 0.9 K in EC-6K with respect to EC-1K (Fig. 5a,b and Table 3). This relative cooling in EC-6K is due to two factors. First, because the atmospheric topography is more smoothed than the topography in the ice-sheet-covered land grid cell, the atmospheric mean elevation per grid cell is lower than the land model mean elevation per grid cell. This gives higher ECs a higher areal weight per grid cell. Second, the characteristic quasi-parabolic shape of the ice sheet contributes to this areal effect. This results in the dominance of the net (negative) energy anomalies from high-elevation ECs. Maximal cooling coincides with areas of rapid change in slope in the SE and NW. Downwind advection of colder air masses from the eastern side of the ice sheet causes mild cooling in the Greenland and the Barents Sea, which is amplified by the growth of sea ice (Fig. 5h).

Table 3Simulated annual (ANN) and summer (JJA) GrIS-averaged components of the surface energy with the standard deviation in parentheses. The period considered is 1965 to 2005. Closest values to RACMO2.3 are given in bold.

Turbulent heat fluxes respond most strongly over the Greenland ice sheet, over the Labrador Sea, and along the sea ice edges in the Greenland and Barents Sea (Fig. 5c, d, and Table 3). Significant differences over the Greenland ice sheet are collocated with areas showing a significant decrease in air temperature. In these simulations, the atmosphere transfers turbulent heat to the surface on average (Fig. 5c). The reduction in air temperature, and consequently air humidity (not shown), results in decreased turbulent heat transfer. Over the Barents Sea, larger sea-ice-covered areas cause a reduction in the heat transfer from the ocean to the atmosphere.

Net surface longwave radiation increases over the Greenland ice sheet where the near-surface temperature decreases (Fig. 5f). Over these areas, incoming longwave radiation decreases; however, this is overcompensated for by a reduction in emitted longwave radiation due to surface cooling.

Figure S2 compares near-surface temperature, turbulent heat fluxes, net longwave radiation, and sea ice extent in EC-1K and EC-6K with ERA-Interim over the entire area in Fig. 5, with the tentative goal of assessing whether the EC method improves or deteriorates the climate simulation. However, the differences between EC-1K and EC-6K are small compared to the difference between these simulations and ERA-Interim, likely due to different realizations of internal climate variability. This precludes a robust conclusion. For Greenland, on the other hand, an assessment is more reliable as the differences between the EC-1K and EC-6K simulations are of the same magnitude as differences with RACMO2.3. The simulation of the GrIS-averaged annual and summer near-surface air temperature is improved in EC-6K, using RACMO2.3 as a reference, as well as the net longwave radiation, melt energy, and (only annual) turbulent heat flux (see bold values in Table 3). The simulated cooling partially counteracts the temperature overestimation in the ESM due to topographic smoothing, resulting in a close fit to RACMO2.3.

4 Discussion

This study has evaluated for the first time the EC method for SMB downscaling from a global climate model of ∼100 km resolution to the much higher resolution (5 km) of an ice sheet model. Other studies (e.g., Alexander et al.2019) have evaluated the effect of implementing ECs on the coarse grid cell but not at the subgrid resolution as done here. This evaluation uses gradients of SEB and SMB components as a primary metric. These gradients are obtained by linear regression of the components on subgrid elevations in all GrIS grid cells. While this provides a systematic framework of comparison, it does not account for relevant non-linear relationships for SMB gradients (e.g., Helsen et al.2012; Noël et al.2016) and SMB components (e.g., precipitation) or for heterogeneity arising from different Greenland climate subregions, local influences on climate (e.g., proximity of tundra, valleys, fjords), or proximity to the ELA.

We justify our comparison with the RCM as dynamical downscaling is the most advanced downscaling technique as shown in numerous evaluations (e.g., Ettema et al.2010; Noël et al.2015). However, one of the limitations of comparing with an RCM is that, unlike an ESM, the RCM is laterally forced with reanalysis. Also, there are fundamental differences in the physical schemes and simulated climate components between the ESM and RCM compared here. Additionally, RACMO2.3 has some well-documented biases, e.g., an underestimation of net longwave radiation, which is compensated for by the sensible heat flux . Further, the RACMO2.3 model was forced at its lateral boundaries by the ERA-Interim reanalysis , which limits the intrinsic or natural climate variability compared to an ESM. Therefore, a more systematic comparison could be made by forcing an RCM with the same ESM where the EC method is implemented.

As a result of the combination of EC downscaling and advanced snow physics , CESM1.0 shows high skill in simulating GrIS climate compared to same-generation global climate models and Earth system models . The ability to realistically represent GrIS SMB in ESMs has been utilized for projections of future SMB change , without an RCM for additional dynamical downscaling. Reliable simulation of the GrIS surface climate at ESM resolution enables exploration of the interaction between the high-resolution surface simulation and other climate components (e.g., atmosphere, ocean, sea ice).

While the EC method in CESM1.0 realistically simulates SMB gradients, we have shown here major deficiencies in the simulation of individual gradients of surface energy and mass balance components compared to the RACMO2.3 RCM. This is an important caveat for modelers who may need to calculate the SMB from individual components of the energy or mass balance, e.g., to perform corrections for one atmospheric forcing field. It also limits the possibility to investigate individual processes at a higher resolution. In the following, we discuss the relative importance and possible fixes of the biases in these individual processes as identified for CESM1.0.

1. CESM1.0 does not capture low enough albedo values due to the use of a single fixed ice albedo, while bare ice has a broader range of albedos . We recommend therefore the use of spatially varying ice albedos, e.g., to simulate the impacts of impurities on ice darkening .

2. The EC scheme in CESM1.0 does not downscale incoming radiation, despite the fact that it varies over small scales at the GrIS surface . The lack of downward longwave downscaling leads to an underestimation of net radiative energy at low-elevation ECs and an overestimation at high-elevation ECs. We recommend downscaling of incoming radiation to reduce overcompensation from the turbulent heat flux gradients and more realistically capture radiation–snow–ice interactions such as shortwave-generated subsurface snowmelt.

3. Since snowfall has no elevation corrections in CESM1.0, small-scale orographically induced precipitation, height-desertification effects, and small-scale variations in the rain-to-precipitation ratio are not captured. Designing realistic and effective elevation corrections for precipitation is a challenging task as the precipitation's correlation with elevation is spatially highly variable over the GrIS . To account for fine-scale variations in the rain-to-precipitation ratio with a simple parameterization, we propose the implementation of a scheme relating the phase of precipitation with atmospheric near-surface temperature, similarly as in .

4. CESM1.0 does not realistically simulate the refreezing gradient, mainly due to limited snow mass in the CLM4.0 snowpack and biased high rainfall rates at high elevations. A realistic simulation of refreezing is key in modeling the response time of an ice sheet to a changing climate as it acts as a buffer for meltwater to run off the ice sheet surface. A more physically based treatment of snow could be used with a snow densification scheme that does not impose a maximum allowed snow depth. An intermediate approach is using relatively large snow and firn depths. As an example along this line, the maximum snow depth can be increased, as in the version 5.0 of CLM, with respect to CLM4.0 due to the further development of the snow scheme to allow for realistic firn simulation .

Assessing the optimal choice of lapse rate forcing proves challenging. In this study, the EC-1K results in the turbulent heat flux gradients closest to RACMO2.3 (Fig. 3a) but almost null melt energy and SMB gradients. EC-4K does not stand out in any way. EC-6K results in the most realistic SMB gradients, despite EC-9.8K comparing the best with RACMO2.3 for the melt gradient. This discrepancy is because CESM1.0 does not downscale snowfall which has an opposite slope to the melt gradient. For the downscaled SMB, EC-6K and EC-9.8K give fairly similar results, making it hard to distinguish one or the other as the best choice. Further improvements of the physical representation of SMB processes at the EC scale might allow for a better identification of an observationally constrained optimal lapse rate.

Global climate models often have warm biases over high areas like the ice sheets, due to topographic smoothing. Here we showed that the EC implementation in CESM1.0 results in moderate cooling over Greenland, which fully compensates for the warm bias with respect to the RCM. The cooling pattern from the EC method is similar to that of , who explored the sensitivity of the simulated GrIS surface climate to horizontal resolution with an RCM.

5 Conclusions

The EC downscaling as implemented in CESM1.0 results in realistic GrIS SMB gradients as shown through comparison with a state-of-the-art RCM. In CESM1.0, high turbulent heat flux gradients compensate for the absence of incoming radiation downscaling. Explicit simulation of snow albedo enables the albedo–melt feedback which is shown to contribute to realistic melt gradients and consequently realistic SMB gradients. Therefore, we conclude that the EC classes method in CESM1.0 efficiently generates a realistic downscaled SMB, despite the fact that only temperature and humidity are downscaled.

Our sensitivity experiments show that a larger lapse rate for the temperature correction results in higher melt energy gradients, as expected. As a consequence of these gradients, ablation areas narrow in CESM1.0, although this result may be different for other models or ice sheet topographies. In turn, this leads to a general cooling downwind of Greenland and an increase in sea ice cover over the Greenland Sea and the Barents Sea. For future implementations of the EC classes within ESMs, we recommend evaluation of the effects on regional climate simulation.

Future improvements of the EC method could be headed towards realistic downscaling of the individual surface energy and mass budget components. Some concrete examples include (1) a lower and/or spatially varying albedo; (2) downscaling of incoming radiation; (3) downscaling of precipitation phase; and (4) development of more adequate snowpack parametrizations for realistic representation of, e.g., snow compaction, firn, and refreezing, fit for polar conditions.

This study aims to guide the future implementation of the EC method, providing diagnostic metrics and evaluation methodology. We recommend in any case that these metrics are adapted to the particular targets of scientific research to be conducted with each model.

Code and data availability
Code and data availability.

The model CESM1.0.5 can be downloaded from http://www.cesm.ucar.edu/models/cesm1.0/ (last access: 21 June 2019). The output from the CESM1.0 simulations, together with code for processing the data and creating the figures for this paper, is available at https://doi.org/10.5281/zenodo.3479410 .

Supplement
Supplement.

Author contributions
Author contributions.

The idea of the study and simulation design came from RS and MV. RS carried out the model simulations, data analysis, and writing of the manuscript, under the supervision of MV. LvK contributed to the development of the analysis software and BN provided RACMO2.3 data. All authors read and commented on the manuscript.

Competing interests
Competing interests.

The authors declare no competing interests.

Acknowledgements
Acknowledgements.

Computing and data storage resources, including the Cheyenne supercomputer (https://doi.org/10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at the National Center for Atmospheric Research (NCAR). The material is based upon work supported by NCAR, which is a major facility sponsored by the National Science Foundation under cooperative agreement no. 1852977. The CESM project is supported primarily by the National Science Foundation. Brice Noël acknowledges funding from the Polar Program of NWO and NESSC. We thank the editor Xavier Fettweis and three anonymous reviewers, whose comments helped improve the manuscript.

Financial support
Financial support.

This research has been supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (grant no. ALWOP.2015.096), the European Research Council (grant no. ERC-StG-678145-CoupledIceClim), and the Netherlands Earth System Science Centre (OCW, grant no. 024.002.001).

Review statement
Review statement.

This paper was edited by Xavier Fettweis and reviewed by three anonymous referees.

References

Alexander, P. M., Tedesco, M., Fettweis, X., van de Wal, R. S. W., Smeets, C. J. P. P., and van den Broeke, M. R.: Assessing spatio-temporal variability and trends in modelled and measured Greenland Ice Sheet albedo (2000–2013), The Cryosphere, 8, 2293–2312, https://doi.org/10.5194/tc-8-2293-2014, 2014. a

Alexander, P. M., LeGrande, A. N., Fischer, E., Tedesco, M., Fettweis, X., Kelley, M., Nowicki, S. M. J., and Schmidt, G. A.: Simulated Greenland Surface Mass Balance in the GISS ModelE2 GCM: Role of the Ice Sheet Surface, J. Geophys. Res.-Earth Surf., 124, 750–765, https://doi.org/10.1029/2018JF004772, 2019. a, b

Bamber, J. L., Griggs, J. A., Hurkmans, R. T. W. L., Dowdeswell, J. A., Gogineni, S. P., Howat, I., Mouginot, J., Paden, J., Palmer, S., Rignot, E., and Steinhage, D.: A new bed elevation dataset for Greenland, The Cryosphere, 7, 499–510, https://doi.org/10.5194/tc-7-499-2013, 2013. a, b

Bamber, J. L., Westaway, R. M., Marzeion, B., and Wouters, B.: The land ice contribution to sea level during the satellite era, Environ. Res. Lett., 13, 063008, https://doi.org/10.1088/1748-9326/aac2f0, 2018. a

Box, J. E. and Rinke, A.: Evaluation of Greenland ice sheet surface climate in the HIRHAM regional climate model using automatic weather station data, J. Climate, 16, 1302–1319, https://doi.org/10.1175/1520-0442-16.9.1302, 2003. a

Bromwich, D. H., Chen, Q.-S., Bai, L.-S., Cassano, E. N., and Li, Y.: Modeled precipitation variability over the Greenland Ice Sheet, J. Geophys. Res.-Atmos., 106, 33891–33908, https://doi.org/10.1029/2001JD900251, 2001. a

Church, J. A., Clark, P. U., Cazenave, A., Gregory, J. M., Jevrejeva, S.,Levermann, A., Merrifield, M. A., Milne, G. A., Nerem, R. S., Nunn, P. D., Payne, A. J., Pfeffer, W. T., Stammer, D. and Unnikrishnan, A. S.: Sea level change, Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, 1137–1216, 2013. a

Cullather, R. I., Nowicki, S. M., Zhao, B., and Suarez, M. J.: Evaluation of the Surface Representation of the Greenland Ice Sheet in a General Circulation Model, J. Climate, 27, 4835–4856, https://doi.org/10.1175/JCLI-D-13-00635.1, 2014. a, b, c

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteorol. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011. a, b

Ettema, J., Van Den Broeke, M. R., Van Meijgaard, E., Van De Berg, W. J., Bamber, J. L., Box, J. E., and Bales, R. C.: Higher surface mass balance of the Greenland ice sheet revealed by high-resolution climate modeling, Geophys. Res. Lett., 36, 1–5, https://doi.org/10.1029/2009GL038110, 2009. a, b

Ettema, J., van den Broeke, M. R., van Meijgaard, E., and van de Berg, W. J.: Climate of the Greenland ice sheet using a high-resolution climate model – Part 2: Near-surface climate and energy balance, The Cryosphere, 4, 529–544, https://doi.org/10.5194/tc-4-529-2010, 2010. a, b, c, d, e

Fausto, R. S., Ahlstrøm, A. P., Van As, D., Bøggild, C. E., and Johnsen, S. J.: A new present-day temperature parameterization for Greenland, J. Glaciol., 55, 95–105, https://doi.org/10.3189/002214309788608985, 2009. a, b

Fettweis, X.: The SMB Model Intercomparison (SMBMIP) over Greenland: first results, available at: http://hdl.handle.net/2268/232923 (last access: 27 November 2019), 2018. a

Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033, https://doi.org/10.5194/tc-11-1015-2017, 2017. a

Fischer, R., Nowicki, S., Kelley, M., and Schmidt, G. A.: A system of conservative regridding for ice-atmosphere coupling in a General Circulation Model (GCM), Geosci. Model Dev., 7, 883–907, https://doi.org/10.5194/gmd-7-883-2014, 2014. a

Flanner, M. G. and Zender, C. S.: Linking snowpack microphysics and albedo evolution, J. Geophys. Res.-Atmos., 111, D12208, https://doi.org/10.1029/2005JD006834, 2006. a, b

Flato, G. M.: Earth system models: an overview, Wiley Interdisciplinary Reviews: Clim. Change, 2, 783–800, https://doi.org/10.1002/wcc.148, 2011. a

Franco, B., Fettweis, X., Lang, C., and Erpicum, M.: Impact of spatial resolution on the modelling of the Greenland ice sheet surface mass balance between 1990–2010, using the regional climate model MAR, The Cryosphere, 6, 695–711, https://doi.org/10.5194/tc-6-695-2012, 2012. a

Fyke, J. G., Weaver, A. J., Pollard, D., Eby, M., Carter, L., and Mackintosh, A.: A new coupled ice sheet/climate model: description and sensitivity to model physics under Eemian, Last Glacial Maximum, late Holocene and modern climate conditions, Geosci. Model Dev., 4, 117–136, https://doi.org/10.5194/gmd-4-117-2011, 2011. a

Fyke, J. G., Vizcaíno, M., Lipscomb, W., and Price, S.: Future climate warming increases Greenland ice sheet surface mass balance variability, Geophys. Res. Lett., 41, 470–475, https://doi.org/10.1002/2013GL058172, 2014a. a, b

Fyke, J. G., Vizcaíno, M., and Lipscomb, W. H.: The pattern of anthropogenic signal emergence in Greenland Ice Sheet surface mass balance, Geophys. Res. Lett., 41, 6002–6008, https://doi.org/10.1002/2014GL060735, 2014b. a, b

Goelzer, H., Huybrechts, P., Fürst, J., Nick, F., Andersen, M., Edwards, T., Fettweis, X., Payne, A., and Shannon, S.: Sensitivity of Greenland Ice Sheet Projections to Model Formulations, J. Glaciol., 59, 733–749, https://doi.org/10.3189/2013JoG12J182, 2013. a

Hanna, E., Huybrechts, P., Janssens, I., Cappelen, J., Steffen, K., and Stephens, A.: Runoff and mass balance of the Greenland ice sheet: 1958–2003, J. Geophys. Res.-Atmos., 110, https://doi.org/10.1029/2004JD005641, 2005. a

Hanna, E., Huybrechts, P., Cappelen, J., Steffen, K., Bales, R. C., Burgess, E., McConnell, J. R., Peder Steffensen, J., Van den Broeke, M., Wake, L., Bigg, G., Griffiths, M., and Savas, D.: Greenland Ice Sheet surface mass balance 1870 to 2010 based on Twentieth Century Reanalysis, and links with global climate forcing, J. Geophys. Res.-Atmos., 116, D24121, https://doi.org/10.1029/2011JD016387, 2011. a

Hanna, E., Navarro, F. J., Pattyn, F., Domingues, C. M., Fettweis, X., Ivins, E. R., Nicholls, R. J., Ritz, C., Smith, B., Tulaczyk, S., Whitehouse, P. L., and Zwally, H. J.: Ice-sheet mass balance and climate change, Nature, 498, 51–59, https://doi.org/10.1038/nature12238, 2013. a

Hartmann, D. L., Klein Tank, A. M., Rusticucci, M., Alexander, L. V., Brönnimann, S., Charabi, Y. A. R., Dentener, F. J., Dlugokencky, E. J., Easterling, D. R., Kaplan, A., Soden, B. J., Thorne, P. W., Wild, M., and Zhai, P.: Observations: Atmosphere and surface, Climate Change 2013 the Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, 9781107057, 159–254, https://doi.org/10.1017/CBO9781107415324.008, 2013. a

Helsen, M. M., van de Wal, R. S. W., van den Broeke, M. R., van de Berg, W. J., and Oerlemans, J.: Coupling of climate models and ice sheet models by surface mass balance gradients: application to the Greenland Ice Sheet, The Cryosphere, 6, 255–272, https://doi.org/10.5194/tc-6-255-2012, 2012. a, b

Hourdin, F., Mauritsen, T., Gettelman, A., Golaz, J.-C., Balaji, V., Duan, Q., Folini, D., Ji, D., Klocke, D., Qian, Y., Rauser, F., Rio, C., Tomassini, L., Watanabe, M., and Williamson, D.: The Art and Science of Climate Model Tuning, B. Am. Meteorol. Soc., 98, 589–602, https://doi.org/10.1175/BAMS-D-15-00135.1, 2017. a

Hunke, E. C., Lipscomb, W. H., Turner, A. K., Jeffery, N., and Elliott, S.: CICE: the Los Alamos Sea Ice Model Documentation and Software User's Manual Version 4.1 LA-CC-06-012, T-3 Fluid Dynamics Group, Los Alamos National Laboratory, 675, 2010. a

Hurrell, J. W., Holland, M. M., Gent, P. R., Ghan, S., Kay, J. E., Kushner, P. J., Lamarque, J. F., Large, W. G., Lawrence, D., Lindsay, K., Lipscomb, W. H., Long, M. C., Mahowald, N., Marsh, D. R., Neale, R. B., Rasch, P., Vavrus, S., Vertenstein, M., Bader, D., Collins, W. D., Hack, J. J., Kiehl, J., and Marshall, S.: The community earth system model: A framework for collaborative research, B. Am. Meteorol. Soc., 94, 1339–1360, https://doi.org/10.1175/BAMS-D-12-00121.1, 2013. a

Jahn, A., Sterling, K., Holland, M. M., Kay, J. E., Maslanik, J. A., Bitz, C. M., Bailey, D. A., Stroeve, J., Hunke, E. C., and Lipscomb, W. H.: Late-twentieth-century simulation of Arctic sea ice and ocean properties in the CCSM4, J. Climate, 25, 1431–1452, 2012. a

Kjeldsen, K. K., Korsgaard, N. J., Bjørk, A. A., Khan, S. A., Box, J. E., Funder, S., Larsen, N. K., Bamber, J. L., Colgan, W., van den Broeke, M., Siggaard-Andersen, M.-L., Nuth, C., Schomacker, A., Andresen, C. S., Willerslev, E., and Kjær, K. H.: Spatial and temporal distribution of mass loss from the Greenland Ice Sheet since AD 1900, Nature, 528, 396–400, 2015. a

Lawrence, D. M., Oleson, K. W., Flanner, M. G., Thornton, P. E., Swenson, S. C., Lawrence, P. J., Zeng, X., Yang, Z.-L., Levis, S., Sakaguchi, K., Bonan, G. B., and Slater, A. G.: Parameterization improvements and functional and structural advances in Version 4 of the Community Land Model, J. Adv. Model. Earth Syst., 3, 1, https://doi.org/10.1029/2011MS00045, 2011. a

Lenaerts, J. T. M., Medley, B., van den Broeke, M. R., and Wouters, B.: Observing and modeling ice sheet surface mass balance, Rev. Geophys., 57, 376–420, https://doi.org/10.1029/2018RG000622, 2019. a

Lipscomb, W. H., Fyke, J. G., Vizcaíno, M., Sacks, W. J., Wolfe, J., Vertenstein, M., Craig, A., Kluzek, E., and Lawrence, D. M.: Implementation and initial evaluation of the glimmer community ice sheet model in the community earth system model, J. Climate, 26, 7352–7371, https://doi.org/10.1175/JCLI-D-12-00557.1, 2013. a, b, c, d, e

Mouginot, J., Rignot, E., Bjørk, A. A., van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.: Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, P. Natl. Acad. Sci. USA, 116, 9239–9244, https://doi.org/10.1073/pnas.1904242116, 2019. a

Neale, R. B., Richter, J., Park, S., Lauritzen, P. H., Vavrus, S. J., Rasch, P. J., and Zhang, M.: The Mean Climate of the Community Atmosphere Model (CAM4) in Forced SST and Fully Coupled Experiments, J. Climate, 26, 5150–5168, https://doi.org/10.1175/JCLI-D-12-00236.1, 2013. a

Noël, B., van de Berg, W. J., van Meijgaard, E., Kuipers Munneke, P., van de Wal, R. S. W., and van den Broeke, M. R.: Evaluation of the updated regional climate model RACMO2.3: summer snowfall impact on the Greenland Ice Sheet, The Cryosphere, 9, 1831–1844, https://doi.org/10.5194/tc-9-1831-2015, 2015. a, b, c, d, e, f

Noël, B., van de Berg, W. J., Machguth, H., Lhermitte, S., Howat, I., Fettweis, X., and van den Broeke, M. R.: A daily, 1 km resolution data set of downscaled Greenland ice sheet surface mass balance (1958–2015), The Cryosphere, 10, 2361–2377, https://doi.org/10.5194/tc-10-2361-2016, 2016. a, b, c

Noël, B., van de Berg, W. J., van Wessem, J. M., van Meijgaard, E., van As, D., Lenaerts, J. T. M., Lhermitte, S., Kuipers Munneke, P., Smeets, C. J. P. P., van Ulft, L. H., van de Wal, R. S. W., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 1: Greenland (1958–2016), The Cryosphere, 12, 811–831, https://doi.org/10.5194/tc-12-811-2018, 2018. a

Nowicki, S. M. J., Payne, A., Larour, E., Seroussi, H., Goelzer, H., Lipscomb, W., Gregory, J., Abe-Ouchi, A., and Shepherd, A.: Ice Sheet Model Intercomparison Project (ISMIP6) contribution to CMIP6, Geosci. Model Dev., 9, 4521–4545, https://doi.org/10.5194/gmd-9-4521-2016, 2016. a

Overland, J., Hanna, E., Hanssen-Bauer, I., Kim, S.-J., Walsh, J. E., Wang, M., Bhatt, U. S., and Thoman, R. L.: Surface air temperature, Arctic report card, 2018. a

Pithan, F. and Mauritsen, T.: Arctic amplification dominated by temperature feedbacks in contemporary climate models, Nat. Geosci., 7, 181–184, 2014. a

Ran, J., Vizcaino, M., Ditmar, P., van den Broeke, M. R., Moon, T., Steger, C. R., Enderlin, E. M., Wouters, B., Noël, B., Reijmer, C. H., Klees, R., Zhong, M., Liu, L., and Fettweis, X.: Seasonal mass variations show timing and magnitude of meltwater storage in the Greenland Ice Sheet, The Cryosphere, 12, 2981–2999, https://doi.org/10.5194/tc-12-2981-2018, 2018. a

Rutt, I. C., Hagdorn, M., Hulton, N. R. J., and Payne, A. J.: The Glimmer community ice sheet model, J. Geophys. Res.-Earth Surf., 114, F2, https://doi.org/10.1029/2008JF001015, 2009. a

Ryan, J. C., Hubbard, A., Stibal, M., Irvine-Fynn, T. D., Cook, J., Smith, L. C., Cameron, K., and Box, J.: Dark zone of the Greenland Ice Sheet controlled by distributed biologically-active impurities, Nat. Commun., 9, 1065, 2018. a

Schmidt, G. A., Bader, D., Donner, L. J., Elsaesser, G. S., Golaz, J.-C., Hannay, C., Molod, A., Neale, R. B., and Saha, S.: Practice and philosophy of climate model tuning across six US modeling centers, Geosci. Model Dev., 10, 3207–3223, https://doi.org/10.5194/gmd-10-3207-2017, 2017. a

Screen, J. A. and Simmonds, I.: The central role of diminishing sea ice in recent Arctic temperature amplification, Nature, 464, 1334–1337, https://doi.org/10.1038/nature09051, 2010. a

Screen, J. A., Deser, C., and Simmonds, I.: Local and remote controls on observed Arctic warming, Geophys. Res. Lett., 39, L10709, https://doi.org/10.1029/2012GL051598, 2012. a

Sellevold, R.: Data and scripts for Surface mass balance downscaling through elevation classes in an Earth System Model: application to the Greenland ice sheet, https://doi.org/10.5281/zenodo.3479410, 2019. a

Serreze, M. C. and Francis, J. A.: The arctic amplification debate, Clim. Change, 76, 241–264, https://doi.org/10.1007/s10584-005-9017-y, 2006. a

Serreze, M. C. and Stroeve, J.: Arctic sea ice trends , variability and implications for seasonal ice forecasting Subject Areas: Author for correspondence, Phil. Trans. R. Soc., 373, 1–16, 2015. a

Shepherd, A., Ivins, E. R., A, G., Barletta, V. R., Bentley, M. J., Bettadpur, S., Briggs, K. H., Bromwich, D. H., Forsberg, R., Galin, N., Horwath, M., Jacobs, S., Joughin, I., King, M. A., Lenaerts, J. T. M., Li, J., Ligtenberg, S. R. M., Luckman, A., Luthcke, S. B., McMillan, M., Meister, R., Milne, G., Mouginot, J., Muir, A., Nicolas, J. P., Paden, J., Payne, A. J., Pritchard, H., Rignot, E., Rott, H., Sorensen, L. S., Scambos, T. A., Scheuchl, B., Schrama, E. J. O., Smith, B., Sundal, A. V., van Angelen, J. H., van de Berg, W. J., van den Broeke, M. R., Vaughan, D. G., Velicogna, I., Wahr, J., Whitehouse, P. L., Wingham, D. J., Yi, D., Young, D., and Zwally, H. J.: A Reconciled Estimate of Ice-Sheet Mass Balance, Science, 338, 1183–1189, https://doi.org/10.1126/science.1228102, 2012. a

Smith, R., Jones, P., Briegleb, B., Bryan, F., Danabasoglu, G., Dennis, J., Dukowicz, J., Eden, C., Fox-Kemper, B., and Gent, P.: The parallel ocean program (POP) reference manual ocean component of the community climate system model (CCSM) and community earth system model (CESM), Rep. LAUR-01853, 141, 1–140, 2010. a

van Angelen, J. H., van den Broeke, M. R., Wouters, B., and Lenaerts, J. T.: Contemporary (1960–2012) Evolution of the Climate and Surface Mass Balance of the Greenland Ice Sheet, Surv. Geophys., 35, 1155–1174, https://doi.org/10.1007/s10712-013-9261-z, 2014.  a

Van den Broeke, M., Smeets, P., Ettema, J., and Munneke, P. K.: Surface radiation balance in the ablation zone of the west Greenland ice sheet, J. Geophys. Res.-Atmos., 113, 1–14, https://doi.org/10.1029/2007JD009283, 2008. a, b, c

van Kampenhout, L., Lenaerts, J. T. M., Lipscomb, W. H., Sacks, W. J., Lawrence, D. M., Slater, A. G., and van den Broeke, M. R.: Improving the Representation of Polar Snow and Firn in the Community Earth System Model, J. Adv. Model. Earth Syst., 9, 2583–2600, https://doi.org/10.1002/2017MS000988, 2017. a

van Kampenhout, L., Rhoades, A. M., Herrington, A. R., Zarzycki, C. M., Lenaerts, J. T. M., Sacks, W. J., and van den Broeke, M. R.: Regional grid refinement in an Earth system model: impacts on the simulated Greenland surface mass balance, The Cryosphere, 13, 1547–1564, https://doi.org/10.5194/tc-13-1547-2019, 2019. a

Van Tricht, K., Lhermitte, S., Gorodetskaya, I. V., and van Lipzig, N. P. M.: Improving satellite-retrieved surface radiative fluxes in polar regions using a smart sampling approach, The Cryosphere, 10, 2379–2397, https://doi.org/10.5194/tc-10-2379-2016, 2016. a

Vizcaíno, M., Mikolajewicz, U., Jungclaus, J., and Schurgers, G.: Climate modification by future ice sheet changes and consequences for ice sheet mass balance, Clim. Dynam., 34, 301–324, https://doi.org/10.1007/s00382-009-0591-y, 2010. a

Vizcaíno, M., Lipscomb, W. H., Sacks, W. J., Van Angelen, J. H., Wouters, B., and Van Den Broeke, M. R.: Greenland surface mass balance as simulated by the community earth system model. Part I: Model evaluation and 1850–2005 results, J. Climate, 26, 7793–7812, https://doi.org/10.1175/JCLI-D-12-00615.1, 2013. a, b, c, d, e

Vizcaíno, M., Lipscomb, W. H., Sacks, W. J., and van den Broeke, M.: Greenland Surface Mass Balance as Simulated by the Community Earth System Model. Part II: Twenty-First-Century Changes, J. Climate, 27, 215–226, https://doi.org/10.1175/JCLI-D-12-00588.1, 2014. a, b

Wientjes, I. G. M., Van de Wal, R. S. W., Reichart, G. J., Sluijs, A., and Oerlemans, J.: Dust from the dark region in the western ablation zone of the Greenland ice sheet, The Cryosphere, 5, 589–601, https://doi.org/10.5194/tc-5-589-2011, 2011. a

Wilton, D. J., Jowett, A., Hanna, E., Bigg, G. R., van den Broeke, M. R., Fettweis, X., and Huybrechts, P.: High resolution (1 km) positive degree-day modelling of Greenland ice sheet surface mass balance, 1870–2012 using reanalysis data, J. Glaciol., 63, 176–193, https://doi.org/10.1017/jog.2016.133, 2017. a