1 Introduction

Global climate models often assume a stable, non-responding, Greenland Ice Sheet (GrIS) or a very simple treatment of ice sheet processes because the ice sheet time-scales (millennial and beyond) are assumed to be much longer than the time-scales of the atmosphere and the ocean (seasonal, annual, decadal). Although studies with ice sheet models have been a rapidly expanding area of research during the last two decades (e.g., Goelzer et al. 2017; Hanna et al. 2020), the interactive coupling between an ice sheet and climate in Earth System Models is still in its infancy (e.g., Vizcaino 2014; Nowicki et al. 2016; Fyke et al. 2018; Muntjewerf et al. 2020a, b).

Observations during the last three decades show a highly variable and increasing mass loss rate from the GrIS (e.g., Bamber et al. 2018; Cazenave et al. 2018; Bevis et al. 2019; The IMBIE Team 2020; Velicogna et al. 2020; Sasgen et al. 2020). A 46-year long record shows both temporal and spatial variability in mass loss with 66 ± 8% of the loss attributed to ice dynamics and 34 ± 8% to the surface mass balance (SMB) for the period 1972–2018 (Mouginot et al. 2019). A reconciled estimate of mass balance from numerous sources for the period 1992–2018 shows that 49.7% of the mass loss was due to increased dynamical imbalance and 50.3% of mass loss was due to increased meltwater run-off (The IMBIE Team 2020) indicating an increased role of surface melt for mass loss. Indirect observations of mass changes over the twentieth century (Box and Colgan 2013; Kjeldsen et al. 2015; Khan et al. 2015, 2020) indicate that the GrIS mass loss has been rapidly increasing after 2000. These rapid changes include a widening of the surface area experiencing melt, a prolonged melting season, and an increase of the interannual variability (e.g., Mote 2007; Tedesco et al. 2013; Tedesco and Fettweis 2020). The temporally evolving rate and location of the mass loss indicate a highly sensitive system (King et al. 2020; Mankoff et al. 2020; Moon et al. 2015, 2020; Mouginot et al. 2019; The IMBIE Team 2020). The rapidly responding GrIS challenges the assumption of fixed ice sheets in global climate models. Instead, it reveals the need to include comprehensive process-based ice sheet models to explore ice sheet-climate interactions under climate change conditions to improve the reliability of climate change projections (Fyke et al. 2018; Hanna et al. 2020).

The interaction between the climate and an ice sheet can introduce significant feedback mechanisms that alter both the climate system and ice sheets (e.g., Fyke et al. 2018). On short timescales (diurnal, seasonal or interannual), the melt-albedo feedback, i.e., variations in surface albedo as a result of surface melting or refreezing, is a major feedback mechanism (Tedesco et al. 2011; Box et al. 2012; Ryan et al. 2019), while the elevation-SMB feedback (e.g., van den Broeke et al. 2017) can have a significant effect on intermediate, multi-decadal, time scales (Edwards et al. 2014; Le clec’h et al. 2019). On longer time scales, ice sheet elevation changes affect precipitation patterns by modifying the atmospheric circulation and moisture transport across the ice sheet (e.g., Solgaard and Langen 2012; Solgaard et al. 2011).

The interaction between the ice sheet and the ocean is important as the increased amounts of freshwater from the melting ice sheet may strengthen the vertical upper ocean stratification, which ultimately weakens the Atlantic Meridional Overturning Circulation (AMOC; Swingedouw et al. 2012, 2014; Bakker et al. 2016; Caesar et al. 2021). This weakened circulation reduces the oceanic heat supply to the Northern North Atlantic region, and eventually dampens Arctic warming (e.g., Gierz et al. 2015; Rahmstorf et al. 2015; Hátún et al. 2005).

During wintertime, seaward of the sea ice edge, a pronounced temperature difference between ocean and atmosphere triggers open-ocean deep water formation via buoyancy loss in cyclonic gyres (Watson et al. 1999; Killworth 1983). In climate simulations, a warming sea surface accelerates sea-ice melting, accompanied by a poleward retreating sea ice edge in the Arctic. In a warming climate, the ocean surface warming dominates, compared to an enhanced freshening, in the amplified stratification of the upper ocean (Gregory et al. 2005). Since additional freshwater dampens the AMOC and reduces northward heat transports, it leads to colder atmospheric conditions in the northern Atlantic and its surrounding regions (e.g., Ivanovic et al. 2018; Lenaerts et al. 2015; Swingedouw et al. 2012, 2014). As the CMIP models assume fixed ice sheets (Eyring et al. 2016), they do not include enhanced calving and run-off from the ice sheet, and therefore potentially underestimate freshwater fluxes to the surrounding ocean. In addition, climate models tend to be biassed towards a state stabilizing the AMOC (Liu et al. 2017).

The integration of an ice sheet model in a global model is challenging because of their different spatial and temporal scales. Consequently, relatively few climate models include a two-way coupling to a dynamical ice sheet model (Ridley et al. 2005; Mikolajewicz et al. 2007a, b; Vizcaino et al. 2008, 2013; Fyke et al. 2011; Gierz et al. 2015; Muntjewerf et al. 2020a; Smith et al. 2021). In the first climate model experiments with an interactive ice sheet (Ridley et al. 2005; Mikolajewicz et al. 2007a; Vizcaino et al. 2008), the climate model supplied near-surface air temperature and precipitation, while the ice sheet model computed the SMB using a positive degree-day model (Reeh 1991). Other studies have applied an energy balance model as an interface between climate and ice sheet model to calculate the SMB at the resolution of the ice sheet model after downscaling the atmospheric forcing from the climate model (Mikolajewicz et al. 2007b; Vizcaino et al. 2010, 2015; Ziemen et al. 2014).

As resolution is well-known to have a large influence on the modelled SMB (e.g., Lucas-Picher et al. 2012; Noël et al. 2016; Langen et al. 2015), a number of methods for downscaling the SMB in global climate models have recently been explored (e.g., Vizcaino et al. 2013; van Kampenhout et al. 2019; Sellevold et al. 2019). Vizcaino et al. (2013) introduced the concept of elevation classes to calculate melt directly in the GCM while still accounting for the height differences between the ice sheet model and the coarser resolved GCM. Other models have used an increased grid resolution for Greenland (van Kampenhout et al. 2019).

Some Earth System Models of Intermediate Complexity (EMICs; e.g., Ganapolski et al. 2010; Golledge et al. 2019) also include ice sheet models. The resolution in these models is much coarser, and the ice sheet forcing is often based on present-day observations in combination with modelled anomalies.

This study presents the EC-Earth-PISM model system with two-way coupling between the global climate model EC-Earth v2.3 (Hazeleger et al. 2012) and the Parallel Ice Sheet Model PISM (Bueler and Brown 2009; Aschwanden et al. 2019; http://www.pism-docs.org; accessed 2021–12-15). The new model system is used to explore the role of an interactive GrIS in the climate response to increased levels of greenhouse gases, i.e., the warm conditions of idealized 4xCO2 scenarios. Section 2 describes the EC-Earth-PISM model system, including the initialization and characteristics of the pre-industrial ice sheet. Results from the model experiments are presented in Sect. 3 and further discussed in Sect. 4.

2 The model system and experiments

In this section, we describe the components that make up the coupled model system. EC-Earth v2.3 is described in Sect. 2.1 and Sect. 2.2 details how EC-Earth was adapted for the coupling. The PISM ice sheet model v0.5 is described in Sect. 2.3 and the details of the coupling scheme are provided in Sect. 2.4. The initialization procedure of the coupled model system is described in Sect. 2.5, followed by a description of the experiments performed (Sect. 2.6). Finally, in Sect. 2.7, the characteristics of the pre-industrial ice sheet are described.

2.1 EC-Earth v2.3

EC-Earth v2.3 (Hazeleger et al. 2012) was the model used by the EC-Earth consortium to contribute to the Coupled Model Intercomparison Project phase 5 (CMIP5). It includes physical processes in the atmosphere, land surface, ocean, and sea-ice. EC-Earth uses the Integrated Forecasting System (IFS, cycle 31r1) of the European Centre for Medium Range Weather Forecasts (ECMWF) for the atmosphere. The land surface component is the Hydrology Tiled ECMWF Scheme of Surface Exchanges over Land (HTESSEL; Balsamo et al. 2009), with an improved snow scheme (Dutra et al. 2010). The atmosphere and land surface components are run on the same grid. The ocean component is NEMO v2.2 (Madec 2008), which embedded the Louvain la Neuve sea-ice model LIM2 (Fichefet and Morales Maqueda 1997; Bouillon et al. 2009).

In the standard configuration, IFS and the HTESSEL are run on a linear N80 reduced Gaussian grid at a spectral horizontal resolution of T159 and 62 vertical layers, corresponding to approximately 125 km horizontal resolution globally. NEMO is run on the ORCA1 grid which has a 1° horizontal resolution with refinement to 1/3° around the equator and 42 vertical layers. The snowpack in EC-Earth v2.3 is represented by a single snow layer, for which the snow temperature is calculated according to the snow energy budget (Dutra et al. 2010):

$$\left[ {\left( {\rho C} \right)_{sn} D_{sn} + L_{f} S_{l}^{c} \frac{{\delta f\left( {T_{sn} } \right)}}{{\delta T_{sn} }}} \right]\frac{{\delta \left( {T_{sn} } \right)}}{\delta t} = R_{sn}^{N} - L_{s} E_{sn} - H_{sn} - G_{sn}^{B} - L_{f} M_{sn}$$
(1)

Here, \({\left(\rho C\right)}_{sn}\) is the snow volumetric heat capacity, \({D}_{sn}\) is the snow depth, \({L}_{f}\) is the latent heat of fusion, \({S}_{l}^{c}\) is the snow liquid water capacity, \({R}_{sn}^{N}\) is the net radiative fluxes, \({L}_{s}\) is the heat of sublimation, \({E}_{sn}\) is the snow evaporation, \({H}_{sn}\) is the sensible heat flux, \({G}_{sn}^{B}\) is the heat flux at the bottom of the snowpack and \({M}_{sn}\) is the snowmelt. Water and ice/snow can co-exist in the snowpack. The second term in the square bracket represents an additional snow heat capacity associated with internal phase changes of the snowpack water (Dutra et al. 2010). For deep snowpacks, a maximum snow depth of 7 cm is used in Eq. (1) to allow heating and melting of surface snow, despite the single-layer snow scheme.

The EC-Earth snow albedo scheme distinguishes between seasonal (less than 1 m water equivalent (m.w.e.)) and permanent snow. For seasonal snow, the albedo decreases due to aging (linear decay) and melting (exponential decay) and varies between 0.5 for wet/old snow and 0.85 for fresh snow. For permanent snow, a fixed value of 0.8 is used.

The ice sheets in the original EC-Earth v2.3 are represented by a permanent layer of snow. A few constraints are applied: (1) for snow depths exceeding 10 m.w.e., snow melts whenever the heat flux is towards the surface, i.e., even at temperatures below the melting point, and the meltwater is returned to the hydrological cycle as run-off. This constraint is meant to represent calving as it redistributes snow from the ice sheet to the ocean; (2) for snow depths exceeding 1 m.w.e., fixed values of 0.8 and 300 kg/m3 are used for snow albedo and density, respectively.

2.2 Adapting the EC-Earth model for coupling to an ice sheet model

In EC-Earth-PISM, the special treatment of permanent snow is replaced by the coupling to the PISM ice sheet model which simulates the ice dynamical and thermo-dynamical processes (Fig. 1). As part of the coupling, the EC-Earth surface scheme has been modified to account for the presence of the ice sheet in the surface energy and mass balance calculations.

Fig. 1
figure 1

The EC-Earth-PISM coupling scheme. Surface mass balance and the temperature at the ice surface (below the snow layer) is calculated in EC-Earth and used to force the PISM ice sheet model. The ice sheet’s topography and extent are calculated in PISM and given back to EC-Earth along with the heat and freshwater fluxes associated with mass loss related to ice discharge and basal melt

Full size image

2.2.1 Land ice in the land surface scheme

The EC-Earth land surface scheme (HTESSEL) differentiates between eight surface types. Each grid box is divided into fractions (tiles) according to the surface type. The land surface tiles include bare ground, low vegetation, high vegetation, intercepted water, exposed snow (on bare land or low vegetation) and snow under high vegetation (Balsamo et al. 2009). Each tile has its own surface properties (albedo, long-wave emissivity and roughness lengths) and for each tile, a separate solution for the surface energy balance equation is found. These solutions are then combined in a weighted average to obtain values for each grid box. Some grid boxes have only one surface type, e.g., snow on bare ground, while other grid boxes may have a combination of several surface types.

As there is no land ice tile, we introduce an ice mask in EC-Earth-PISM. The ice mask is based on the ice sheet thickness in PISM (see details in Sect. 2.4). A grid-box in the EC-Earth land surface scheme is considered glacierized if at least half the grid box intersects with PISM grid boxes defined as ice sheet. As an entire EC-Earth land surface grid box is either seen as glazierized or not, we note that this is a rough approach, especially when the resolution difference is large. For glacierized grid boxes in the EC-Earth land surface scheme, ice surface properties are adopted for all snow-free tiles: i.e., we use an albedo of 0.6, a long-wave emissivity of 0.98, and roughness lengths of 1 mm. Note that the surface characteristics are not affected for snow tiles, i.e., glacierized grid boxes may be partly covered by snow.

Below the snow, the vertical profiles of soil temperature and moisture are represented by four layers extending to a depth of 2.89 m (Viterbo and Beljaars 1995). In EC-Earth-PISM, these layers are treated as ice for glacierized grid boxes, i.e., the ice heat conductivity (2.2 W/mK) and volumetric heat capacity (2.05 J/m3K) replace the soil values otherwise used. In this way, the presence of the ice sheet modifies the vertical distribution of heat and the exchange of energy with the overlying snow or atmosphere.

2.2.2 Snow on ice sheets

For glacierized grid boxes, the fixed albedo used in EC-Earth is replaced by a simple albedo scheme where the albedo decays exponentially under wet conditions with a minimum value of 0.6; for non-melting conditions, the snow albedo is 0.80 as in the standard model. The decrease in albedo under wet conditions allows for the influence of the melt-albedo feedback. The value of 0.6 used in EC-Earth-PISM for the albedo of bare ice and melting snow is a rather conservative estimate. For comparison, Wehrlé et al. (2021) find a bare ice albedo of 0.565 at ice-ablation onset at Greenland weather stations, and Ryan et al. (2019) estimated a bare ice albedo in the range 0.45–0.57. The snowmelt is sensitive to the parameters in the snow albedo scheme, and more elaborate parameterizations may improve the estimated SMB (Helsen et al. 2017).

In EC-Earth-PISM, the snow energy balance calculation takes glacierized grid boxes into account by using ground properties of ice instead of soil. In the snow energy balance calculation (Eq. 1), the basal heat flux at the bottom of the snowpack is calculated as

$$G_{sn}^{B} = \frac{{T_{sn} - T_{soil,ice} }}{{r_{sn} }},$$
(2)

where \({r}_{sn}\) is the thermal resistance between the middle of the snowpack and the middle of the upper soil/ice layer. For glacierized grid points, this resistance is calculated from the conductivity of ice and the presence of the underneath ice sheet is thereby considered.

2.2.3 Ice sheet surface melt

In the snow-free part of the grid box, the exposed ice surface melts when the skin temperature, i.e., the temperature of the surface in radiative equilibrium, reaches 273.16 K, and energy is available at the surface. The surface energy balance equation assumes a zero heat capacity of the skin layer and is solved separately for each tile

$$\left(1-{\alpha }_{i}\right){{R}_{s}}^{\downarrow }+{\varepsilon }_{i}\left({{R}_{T}}^{\downarrow }-\sigma {T}_{sk,i}^{4}\right)+{H}_{i}+{L}_{v,s}{E}_{i}={\Lambda }_{sk,i}\left({T}_{sk,i}-{T}_{s}\right)$$
(3)

where the subscript, i denotes the tile index. \({\alpha }_{i}\) is the albedo, \({{R}_{s}}^{\downarrow }\) and \({{R}_{T}}^{\downarrow }\) are the downward short-wave and longwave radiation, respectively, \({\varepsilon }_{i}\) is the surface emissivity, σ is the Stefan–Boltzman constant, \({T}_{sk,i}\), is the skin temperature, \({H}_{i}\) the sensible heat flux, \({L}_{v,s}\) the latent heat of vaporization/sublimation, \({E}_{i}\) the evaporation from the skin layer, \({\Lambda }_{sk,i}\) the skin conductivity for tile i and \({T}_{s}\) is the temperature of the uppermost ice layer (Viterbo and Beljaars 1995). The right-hand side represents the conductive heat flux: \({G}_{i}\). If the calculated skin temperature reaches the melting point temperature, Tsk is set to 273.16 K, the conductive heat flux is recalculated, and the excessive energy is the energy available for melting the ice, which is calculated for each surface tile as

$${Q}_{Melt,i}=\left(1-{\alpha }_{i}\right){{R}_{s}}^{\downarrow }+{\varepsilon }_{i}\left({{R}_{T}}^{\downarrow }-\sigma {T}_{sk,i}^{4}\right)+{H}_{i}+{L}_{v,s}{E}_{i}-{G}_{i}$$
(4)

The total available energy is summed over the snow-free land surface tiles and multiplied by the ice mask (see Sect. 2.4), as ice sheet surface melt should only occur in the part of the grid box which is ice sheet in PISM. Surface ice melt is then calculated, dividing by the latent heat of fusion, and added to the snowmelt.

Meltwater from snow and ice is assumed to run off immediately following the EC-Earth routing scheme, i.e., run-off is collected in pre-defined river basins and instantaneously distributed into ocean points near outlets of major rivers (Hazeleger et al. 2012). Refreezing is not implemented in the present version of EC-Earth-PISM.

2.3 The PISM ice sheet model

The ice sheet model PISM (version 0.5) is a 3D hybrid stress balance model that combines the solutions of the shallow shelf and shallow ice approximations (SSA and SIA, respectively) to compute the gravitational flow and horizontal stretching of ice (Bueler and Brown 2009). A pseudo-plastic power law (Schoof and Hindmarsh, 2010) relates bed-parallel shear stress and the sliding velocity (Aschwanden et al. 2016, 2019). The thermodynamic equation considers an enthalpy rather than a temperature formulation to represent polythermal ice (Aschwanden et al. 2012). PISM uses a modified version of the Lingle-Clark bedrock deformation model (Lingle and Clark 1985) assuming an elastic lithosphere and a resistant asthenosphere with viscous flow in the half-space below the elastic plate (Bueler et al. 2007). At the base of grounded ice, basal melting is driven by both ice-internal shear stress and geothermal heat flux, while the pressure dependence of the local melting temperature accounts for the gravitational ice load. The ice sheet model requires bedrock topography and geothermal heat flux at the bottom of the ice sheets as boundary conditions besides SMB and ice surface temperature. Ocean thermal forcing is assumed to be constant in space and time. The sub-shelf ice temperature is set to the pressure melting point which ultimately depends on the ice geometry. The sub-shelf melt rate is proportional to the heat flux that is driven by the temperature difference between the ice and the pressure melting point. We assume a constant proportional coefficient of 0.5 W/m2 between heat flux and temperature difference. The ocean thermal forcing contributes to basal melting of floating ice but does not directly impact calving.

Around Greenland, current observations reveal that ice does not advance into the deep ocean. Therefore, ice that flows across the continental shelf break calves at this shelf break in all our simulations. The continental shelf break is assumed to follow the 100 m depth contour in the SeaRISE bathymetry data set (Bindschadler et al. 2013). Hence, calving does not occur directly at the coastline.

PISM has been used for modelling Greenland (Aschwanden et al. 2013, 2019; Aðalgeirsdóttir et al. 2014; Goelzer et al. 2018) and Antarctic (e.g., Winkelmann et al. 2011; Rodehacke et al. 2020) ice sheets in several studies. Results have shown that the model realistically simulates the past evolutions and present state of the Greenland and Antarctic ice sheets. Therefore, PISM is applicable for coupling with EC-Earth to study the future changes of the ice sheets.

2.4 The coupling scheme

The coupling of EC-Earth and PISM occurs at the script level with exchange of files between the two models once per simulation year (Fig. 1). EC-Earth provides the forcing for the ice sheet model as monthly fields of SMB and subsurface temperature. The temperature of the lowermost ice layer in EC-Earth, which is located a few meters below the surface, is interpreted as the ice surface temperature in PISM. The temperature is interpolated to the PISM grid (bilinear interpolation) while considering a uniform lapse rate correction of − 6.8 K/km to account for elevation differences between the two models (Fausto et al. 2009). Precipitation, evaporation (including sublimation), and surface run-off determine the SMB:

$$SMB=Precipitation-evaporation-runoff$$

The SMB is interpolated from the IFS grid to the PISM grid using conservative remapping without any downscaling or corrections of elevation differences between the two models. This simple procedure ensures that the SMB is consistent with the atmospheric model, while sub-grid scale changes in the ablation zone are not accounted for.

EC-Earth receives ice thickness and surface topography (via ice sheet thickness changes and Glacial Isostatic Adjustment) from PISM. On the finer PISM grid we apply an ice thickness threshold of 1 m (water equivalent) to distinguish between seasonal and perennial ice (to determine an ice mask with values of zero and one). After regridding of this ice mask from the higher resolved ice sheet grid to the coarser EC-Earth land surface grid, we have a fractional ice mask in EC-Earth. In the EC-Earth land surface scheme, grid boxes which are at least 50% ice sheet are treated as glacierized, i.e., ice properties are used in the energy balance equations (Sect. 2.2.1). The fractional ice mask is applied in the calculation of ice sheet surface melt which is proportional to the ice sheet fraction (Sect. 2.2.3).

EC-Earth also receives calving and basal melt fluxes from PISM. Calving produces ice bergs which release freshwater into the surface ocean, while the energy required for melting the related ice berg is subtracted from the surface ocean. Once the surface water reaches the freezing temperature, further energy subtraction eventually forms sea ice. Since the PISM model grid generally has a much higher horizontal resolution than the ocean and atmosphere model, fjords resolved in the ice sheet model as ocean points usually correspond to land points in the ocean and atmosphere grid. Therefore, heat and freshwater fluxes due to calving are remapped to the nearest grid box in the ocean model NEMO. These are located along Greenland’s coast. Accordingly, any freshwater flux from basal melting of grounded and floating ice is mapped to the nearest EC-Earth land point, from where it is redistributed to the ocean by the EC-Earth routing scheme which also redistributes freshwater from surface run-off.

2.5 Initialization of EC-Earth-PISM

For the experiments with EC-Earth-PISM, the ice sheet is assumed to have been in equilibrium with the climate under the pre-industrial conditions.

To achieve this, we start with the state of the GrIS obtained from a bootstrapping method as described in the PISM user’s manual (The PISM authors 2015) using present-day observed data of bedrock topography and ice sheet surface altitude (SeaRISE, Bindschadler et al. 2013). This ice sheet state is then run through a pre-spinup simulation using the stand-alone PISM, for which a monthly climatology of the pre-industrial EC-Earth simulation is applied as a constant forcing. After 100 ka, a quasi-stationary ice sheet state is achieved, and this pre-spin-up/stand-alone state of the ice sheet is used as the ice sheet initial state in the ensuing spin-up of the fully coupled EC-Earth-PISM system.

The initial state of the atmosphere, ocean and sea ice for the spin-up of the coupled model system is taken from an EC-Earth initial state, which has gone through a multi-century spin-up and control simulation and reached a quasi-equilibrium state under pre-industrial conditions. The EC-Earth-PISM model system is then spun-up from this model state together with the pre-spinup ice sheet described above and under pre-industrial conditions. The interactive ice sheet appears to have no significant effects on the overall heat balance at the top and the bottom of the atmosphere, and thus no striking changes occur in the modeled global climate. After 130 years of coupled-model spin-up, the emerging model state is considered stable and is used as the initial state for the three specific EC-Earth-PISM experiments described below. Fig. S1 demonstrates the stability of the global mean temperature in the pre-industrial experiment, continuing from the coupled spin-up.

The modelled ice sheet shape is, to a great extent, similar to the present-day observed ice sheet with the maximum ice thickness over the Greenland Summit and a second maximum on the South Dome (Fig. 2a, initial elevation). The modelled ice sheet is larger than the observed ice sheet. It covers a larger area (2.3·106 km2 versus 1.8·106 km2 in Kargel et al. 2012) and extends to the coast in most regions, such as the Northern part of Greenland that is currently ice-free. The ice sheet volume is too large compared to the present-day observed state (3.84 ·106 km3 versus 2.99 ± 0.02 ·106 km3 in Morlighem et al. 2017), mostly due to the thicker ice along the margins of the ice sheet. It is a common problem that modelled ice sheets are too large compared to observations, when the SMB forcing is very different from the observed (e.g., Vizcaino et al. 2010; Nowicki et al. 2013; Goelzer et al. 2018). In our case, where the SMB is calculated at the coarse resolution of the atmosphere model that does not resolve the small-scale topography of the ablation zone, this would be the main contribution to the large ice sheet. The inherent bias of the model climate and the simplified representation of processes in the ice sheet model itself (e.g., stress regime at the margin, iceberg calving) and the resolution of the ice sheet model (i.e., 20 km) in EC-Earth-PISM may contribute to the large ice sheet as well. Since small ice caps and glaciers along the ice sheet’s margin cannot be resolved sufficiently by the ice sheet model resolution, these ice caps merge with the main ice sheet and contribute to its large size.

Fig. 2
figure 2

Evolution of the surface elevation over time in a piControl and b Abrupt4xCO2. The initial surface elevation is shown as the gray-scale background field, while color-coded contour lines of the 1000 m, 2000 m and 3000 m show the temporal evolution (see legend for the color-year relation). The lower row compares c satellite-based velocity estimates from Joughin et al. 2018 with d the simulated mean surface velocities of the last thirty years of piControl

Full size image

2.6 Experiments

The EC-Earth-PISM system is configured to use the EC-Earth standard resolution for the atmosphere and ocean. The atmosphere has a spectral resolution of T159 horizontally and 62 vertical layers on a linear N80 reduced Gaussian grid corresponding to approximately 125 km horizontal resolution globally. The ocean uses the ORCA1 grid which has a 1° horizontal resolution with refinement to 1/3° around the equator and 42 vertical layers. The ice sheet model utilizes a regular polar-stereographic grid at a 20 km resolution.

To evaluate the performance of the EC-Earth-PISM coupled model and to explore the influence of the interactive ice sheet under warm climate conditions, one control and two climate warming experiments are performed with both EC-Earth-PISM and EC-Earth. These three experiments are designed as:

  1. 1.

    piControl: a control simulation with all forcings (i.e., greenhouse gas concentrations, aerosol loadings, and the solar forcing) held constant at the pre-industrial level (CO2 concentration of 285 ppm).

  2. 2.

    Abrupt4xCO2: a simulation with the pre-industrial atmospheric CO2 concentration abruptly quadrupled to 1140 ppm, following the CMIP5 protocol.

  3. 3.

    1pctCO2: a simulation with the atmospheric CO2 concentration gradually increased by 1% per year until it reaches four times the pre-industrial level (after 140 years) and stabilizes after that.

In the EC-Earth experiments without the interactive ice sheet, a fixed present-day topography was used for Greenland. All experiments are run for 350 model years.

2.7 Characteristics of the pre-industrial climate and ice sheet

The climate of EC-Earth in pre-industrial and present-day control experiments was evaluated in Hazeleger et al. (2012) and Sterl et al. (2012). Overall, the large scale characteristics of the atmosphere are well represented, but the atmosphere in general has a cold bias (Hazeleger et al. 2012). The northern hemispheric ocean has a cold and fresh bias, while sea-ice extent and concentrations are simulated realistically. The AMOC is relatively low compared to observational estimates (Sterl et al. 2012). In the twentieth century, the Arctic is found to have a cold bias of 2 K, while sea-ice extent and concentration is overestimated (Koenigk et al. 2013). As compared to the CMIP5 model ensemble, the EC-Earth bias in near surface air temperature (precipitation) is about the same as (~ 10% smaller than) the median bias of all CMIP5 models (Flato et al. 2013).

For the ice sheet, Fig. 2a shows the temporal evolution of its surface elevation in piControl starting from the initial state. Throughout the 350 years of experiment, only small inland movements of the 1000 and 2000 m contour lines occur at some locations in southern Greenland. The modeled velocity at the end of the piControl experiment is shown in Fig. 2d. It compares relatively well with the MEaSUREs velocity (Fig. 2c, Jougin et al. 2018), the fast flowing regions are well simulated considering the spatial resolution of 20 km. Since the simulated ice sheet is larger and thicker, this may cause the expanded region of low velocities near the ice divides (Fig. 2c and d). The ice stream feeding the 70 N glaciers (Nioghalvfjerdsbrae) is underrepresented in our model, probably because the special basal conditions, such as the network of supraglacial streams, are only rudimentarily represented at this model resolution.

The annual mean surface mass balance, the ice surface temperature and the total ocean flux (i.e., calving and thermal interactions) calculated in PISM is illustrated in Fig. 3a, c and e as an average over the last 30 years of the piControl experiment. The SMB pattern reflects the resolution of the atmospheric model. It captures the high accumulation zone in southeastern Greenland and the negative SMB in Southern Greenland and along the central part of the west coast. The annual mean SMB is positive in the south-western part of the ablation zone, indicating a low amount of surface run-off. Regions of underestimated surface ablation are characterized by an ice margin that progresses until the coastline, where ice discharge compensates for the insufficient surface run-off (Fig. 3e). The ice sheet elevation has a strong control on the ice surface temperature (Fig. 3c) such that higher elevated areas are generally colder. The modeled pre-industrial ice sheet is in quasi-equilibrium with the EC-Earth model climate. During the 350 years of the piControl simulation, it has a small drift of -0.06 mm SLE/year (Fig. 4c), which integrates to 2.1 cm of potential sea-level fall. This drift is small compared to the estimated sea-level contribution over the last century of 0.6 mm SLE/year (Mitrovica et al. 2001) and is not subtracted from the results presented in the following section.

Fig. 3
figure 3

The subpanels depict 30-year mean values of surface mass balance, ice surface temperatures, and total ocean flux as simulated by EC-Earth-PISM at the end of the experiments (year 321–350). a and b show the surface mass balance calculated by the atmospheric module. The surface mass balance is regridded to the ice sheet model grid by conserving the total flux. c and d show the ice surface temperature below the snow layer as seen by PISM. e and f show the total ocean flux due to iceberg calving and ice terminus melting

Full size image
Fig. 4
figure 4

Time series of a five-year running means of SMB (dashed lines) and calving (solid lines), b five-year running means of total freshwater flux into the ocean c annual mean ice sheet area and d ice sheet volume given in meters of potential sea-level rise in the three EC-Earth-PISM experiments. The same initial state is used in all experiments. A total ocean area of 3.619 million km2 (Eakins and Sharman 2010) is used to convert from volume in m3 to m/SLE

Full size image

3 Results

This section discusses the climate response to increased levels of CO2 in order to assess the influence of an interactive GrIS on the response of the climate system. The GrIS mass balance changes are described in Sect. 3.1, and the influence of these changes on the Arctic and North Atlantic regions are discussed in Sects. 3.2–3.4. For the comparison between EC-Earth and EC-Earth-PISM, the average of the last 50 years of each experiment (years 301–350) is considered, i.e. after a stable climate state has been achieved in the two climate warming experiments. In each grid point, differences are considered significant, if the difference between the 50-year averages (301–350) is larger than one standard deviation of the same 50 years in the reference experiment.

3.1 Mass balance of the GrIS

The components contributing to the mass balance of the GrIS are shown in Table 1 for all experiments as averages over the model years 301–350. All surface fluxes are calculated on the EC-Earth land surface grid and integrated over the GrIS using the EC-Earth-PISM fractional ice mask. For the three EC-Earth-PISM experiments, the GrIS integrated time evolution of SMB, calving and freshwater fluxes are shown in Fig. 4a and b. The evolution of total area and volume are shown in Fig. 4c and d, and maps of SMB, ice surface temperature and total ocean flux, as seen by the ice sheet model, are shown in Fig. 3 for piControl and Abrupt4xCO2. Note that the total ocean flux is dominated by calving but includes an insignificant contribution from ocean thermal forcing.

Table 1 The mass balance components (Gt/year) integrated over the GrIS as simulated by EC-Earth-PISM and EC-Earth and averaged over model years 301–350
Full size table

3.1.1 SMB components

Under pre-industrial conditions, the SMB of EC-Earth-PISM is 531 ± 57 Gt/year (Table 1). The amounts of precipitation and surface run-off vary substantially from year to year, and the resulting SMB has a large interannual variability. For comparison, van Kampenhout et al. (2020) estimate an SMB of 508 ± 73 Gt/years (1961–1990) using the CESM2 model, and Muntjewerf et al. (2020b) report an SMB of 544 ± 103 Gt/year under pre-industrial conditions using the CESM2.1-CISM2.1 coupled ice sheet climate model. Fettweis et al. (2020) estimate an SMB of 338 ± 111 Gt/year for 1980–2012 using an ensemble of models, including both physically based climate models (RCMs and GCMs), energy balance models and positive degree day (PDD) models. As described below, the relatively high SMB in EC-Earth-PISM under pre-industrial conditions results from a relatively low estimate of surface melt and run-off, while the amount of precipitation compares well with other studies. It should also be noted, that the EC-Earth-PISM total SMB is estimated using an ice mask which is larger than the mask used in Fettweis et al. (2020) and Muntjewerf et al (2020b).

In piControl, precipitation on the GrIS accounts for 743 ± 53 Gt/year in EC-Earth-PISM. About 85% (634 Gt/year) of the precipitation falls as snow, and this amount compares well with the present-day (1980–2012) estimate of 642 ± 59 Gt/year in Fetttweis et al. (2020). The total amount of precipitation on the GrIS is very similar in EC-Earth and EC-Earth-PISM (Table 1; Fig. 5a and b). The differences on the GrIS are patchy and very small (less than ± 0.6 mm/day), even though significant, as can be seen in Fig. 5c. These small differences in precipitation are likely due to the application of the initial ice sheet topography in EC-Earth-PISM that results in higher and steeper topography on the margin of the ice sheet than the topography in EC-Earth (Fig. S2). Outside Greenland, there are no significant differences between the precipitation patterns in EC-Earth and EC-Earth-PISM.

Fig. 5
figure 5

Left: Precipitation (P, in mm/day) in the a EC-Earth-PISM and b EC-Earth piControl experiments, and c the difference between (a) and (b). Right: Change in precipitation in Abrupt4xCO2 with respect to piControl in d EC-Earth-PISM, e EC-Earth and f the difference between (d) and (e). All values are averaged over the years 301–350 (grey shading in Fig. 6). In each subpanel, the hatching marks regions where the difference exceeds one standard deviation of the respective reference simulation

Full size image

In piControl, the surface run-off in EC-Earth-PISM is 154 ± 31 Gt/year. This amount is less than half of the estimated present-day (1980–2012) surface run-off of 331 ± 102 Gt/year in Fettweis et al. (2020). The relatively low estimate of surface run-off is in agreement with EC-Earth’s cold bias in the Arctic (Koenigk et al. 2013) and the model system’s horizontal resolution which is not sufficient to resolve the highly varying conditions in the ablation zone. The ice surface temperature as shown in Fig. 3c shows the cold conditions over the ice sheet. In EC-Earth, the surface run-off is 576 ± 45 Gt/year, more than three times larger than in EC-Earth-PISM. This very large contribution from surface run-off is due to the parameterized melting of permanent snow (> 10 m.w.e.), which is meant to substitute calving in the model without the interactive ice sheet. The EC-Earth surface run-off thus effectively represents mass loss from both surface run-off and calving.

In the climate warming experiments, the hydrological cycle is enhanced and the SMB increases over central Greenland and decreases along the margin (Table 1 and Fig. 4a). The SMB has a large interannual variability, reflecting the pronounced year to year variations in ablation (e.g., Fyke et al. 2014). The strong forcing results in significant increases in precipitation over the whole Arctic region for both models with the largest increase appearing in the high accumulation zone in the southeastern part of Greenland and the smallest in central and northern Greenland (Fig. 5d and e). However, the differences in the response of precipitation to the warming between EC-Earth-PISM and EC-Earth are small and only appear in the south of Greenland and the surrounding ocean. These differences are not statistically significant (Table 1 and Fig. 5f). Due to the increase in precipitation, the ice sheet thickens in the interior of Greenland and other regions of high elevation (Fig. 2b). The precipitation integrated over the GrIS increases 67% in 1pctCO2 (averaged over year 301–350 compared to piControl) which is large compared to the precipitation changes of 32% found in Muntjewerf et al (2020b) for a similar scenario. This reflects the large model uncertainty in the regional precipitation response to warming (Collins et al. 2013).

In the warming experiments with EC-Earth-PISM, the surface run-off increases by a factor of eight compared to the pre-industrial (Table 1). Surface melt increases along the ice margin and the SMB turns negative in most coastal grid boxes outside the high-accumulation zone. Calving occurs at the same locations as in piControl but the calving rate reduces in areas with increased surface melting (Fig. 3e and f). The melt-albedo feedback drives increased snowmelt which leaves a larger part of the ice sheet snow-free and exposed to ice surface melt. In EC-Earth, the surface run-off increases by less than half the amount in EC-Earth-PISM, as the fixed albedo of permanent snow does not allow for the melt-albedo feedback to be activated.

For comparison, Muntjewerf et al. (2020b) used the CESM2.1-CISM2.1 coupled ice sheet climate model to simulate the mass loss of the Greenland Ice Sheet in a 1pctCO2 scenario. For this model, the surface melt at the end of the 350-year long experiment was 3804 Gt/year which is 3–4 times the amount simulated in EC-Earth-PISM.

3.1.2 Freshwater flux from the GrIS

In EC-Earth-PISM, the freshwater flux to the ocean is the sum of surface run-off, calving, and basal melt (see Table 1). The evolution of the freshwater flux in the EC-Earth-PISM experiments is shown in Fig. 4b. In EC-Earth, there is no calving/basal melt; instead, the surface run-off effectively equals the freshwater flux to the ocean, as it includes the contribution from melting of snow exceeding 10 m.w.e.

In piControl, the total freshwater flux to the ocean from the GrIS is 682 ± 49 Gt/year in EC-Earth-PISM. Calving and basal melting contribute 511 ± 27 and 17 ± 1 Gt/year, respectively, and surface run-off the remaining 154 ± 31 Gt/year. The modelled calving and basal melt compares well with recent estimates of Greenland discharge (Mouginot et al. 2019; King et al. 2020; Mankoff et al. 2020), e.g. the discharge increased after 2000 from 458 ± 10 Gt/year (1972–1980) to 537 ± 5 Gt/year (2010–2018; Mouginot et al. 2019). In EC-Earth-PISM, the simulated ice sheet reaches a quasi-equilibrium because the too intense calving rate compensates for the underestimated run-off flux. However, the partitioning between surface run-off and calving is not well captured in EC-Earth-PISM. As a result of the modified albedo scheme and the additional freshwater sources (calving and basal melt from PISM), the freshwater flux from the GrIS to the ocean is 18% larger in EC-Earth-PISM than in EC-Earth during the last 50 years of the piControl (682 ± 49 Gt/year) vs (575 ± 45 Gt/year).

Compared to the pre-industrial period, the freshwater flux from the ice sheet into the ocean more than doubles in the climate warming experiments with EC-Earth-PISM (Table 1 and Fig. 4d). The small ablation from basal melting of grounded ice doubles (from 17 ± 1 Gt/year to 29 ± 2 and 39 ± 3 Gt/years, in 1pctCO2 and Abrupt4xCO2, respectively). Calving reduces by about 10% (from 511 ± 27 to 468 ± 33 and 442 ± 75 Gt/year, in 1pctCO2 and Abrupt4xCO2, respectively) as the ice sheet is thinning and starts to retreat from the coast (Fig. 3e and f). The freshwater flux from the GrIS increases about 55% more in EC-Earth-PISM than in EC-Earth as a result of the interactive ice sheet that takes the melt-albedo feedback into account.

3.1.3 Total mass balance

The pre-industrial ice sheet is not fully stable but gains mass by + 10 Gt/year (Table 1) throughout the 350 years of experiment as estimated on the IFS grid. One should note that estimating on the finer grid of PISM results in a slightly bigger mass gain of ~ 23 Gt/year or a contribution of − 0.06 mmSLE/year during the entire pre-industrial experiment (Fig. 4d). The coupling between EC-Earth and PISM is thus not fully mass-conserving. Due to the implementation of the ice mask, mass conservation is not guaranteed in grid boxes which are partially ice sheet, and the mass received on the finer PISM grid is slightly higher than that on the original coarser grid of the IFS model.

The ice sheet initially loses mass at very different rates in Abrupt4xCO2 and 1pctCO2, consistent with the warming rates of the respective experiment. Averaged over the whole experiment, the ice sheet loses mass at a rate of 504 and 360 Gt/year in Abrupt4xCO2 and 1pctCO2, respectively. However, after the first 140 years, the mass loss rate is very similar in the two experiments (cmp., Fig. 4d). Once the CO2 level has reached its maximum, the GrIS loses mass at a nearly constant rate of 1.4 mm SLE/year in both warming experiments (Fig. 4d). Averaged over the last 50 years of the experiments, the net mass loss is 499 Gt/year in Abrupt4xCO2 and 532 Gt/year in 1pctCO2 (Table 1). At the end of the simulation, the ice loss under Abrupt4xCO2 amounts 4.9% of its initial mass.

In the EC-Earth climate warming experiments, there is a net increase of mass over the ice sheet domain, as the amount of precipitation increases more than the surface run-off. As the melt-albedo feedback is not active in the permanent snow regions, snowmelt increases significantly less than in EC-Earth-PISM (Table 1). Furthermore, there is no calving in EC-Earth and the cut-off of snow above 10 m works less well in the climate change scenarios with increased precipitation. In piControl, snow accumulation above 10 m.w.e. occurred in a few grid boxes; this occurs more frequently and to a larger extent in the warm scenarios. These results show that EC-Earth cannot simulate a realistic Greenland mass balance under climate change conditions and thus demonstrates the importance of the interactive ice sheet.

3.2 Atmospheric response

As mentioned in Sect. 2.5, the interactive ice sheet does not change the global feature of the climate (Fig. S1). However large impacts can be seen in Northern high latitudes, in particular in the climate warming experiments.

Under pre-industrial conditions, the mean Arctic near-surface air temperature (north of 60° N) is stable and very similar in the EC-Earth and EC-Earth-PISM experiments in both the mean values and the variability (Fig. 6). Figure 7a and b show maps of near-surface air temperature in the EC-Earth-PISM and EC-Earth piControl experiments. The local pattern of near-surface air temperature over Greenland highlights the ice sheet topography differences between the EC-Earth and EC-Earth-PISM experiments (Fig. 7c). The GrIS in EC-Earth-PISM is broader, accompanied by a higher elevation along the ice sheet’s periphery (Fig. S2), and therefore lower temperatures occur at the ice sheet’s margin.

Fig. 6
figure 6

Time series of annual mean near-surface air temperature (TAS, in °C) for the EC-Earth and EC-Earth-PISM experiments. Values are averaged over all grid cells north of 60° N. Thin lines indicate annual means, and thick lines are 11-year running means. The grey shading illustrates the 50-year period shown in Table 1 and Figs. 5, 7, 8b-d, 9 and 10

Full size image
Fig. 7
figure 7

Left: Near-surface air temperature (TAS, in °C) in the a EC-Earth-PISM and b EC-Earth piControl experiments and c the difference between (a) and (b). Right: Change in near-surface air temperature in Abrupt4xCO2 compared to piControl in d EC-Earth-PISM and e EC-Earth and f the difference between (d) and (e). All values are averaged over the years 301–350 (grey shading in Fig. 6). In each subpanel, the hatching marks regions where the difference exceeds one standard deviation of the respective reference simulation. In subpanels (d) and (e) the entire shown area exceeded one standard deviation so the hatching has been removed to improve the visibility

Full size image

In the climate warming experiments, the near-surface air temperature initially evolves similarly in experiments with and without the interactive ice sheet (Fig. 6). In Abrupt4xCO2, the rapid increase of the radiative forcing results in an immediate increase of about 3 K during the first year and about 10 K in the first 20 years (Fig. 6). The increase gradually slows down, and the temperature stabilizes after 200 years. In 1pctCO2, the near-surface air temperature rises steadily following the pace of the amplified forcing. The warming rate reduces once the final atmospheric carbon concentrations have been reached (year 140) and, eventually, the temperature stabilizes at a level similar to the warming in Abrupt4xCO2. During the first 200 years there is only a small difference between the near-surface air temperature in EC-Earth-PISM and EC-Earth but the stabilized temperature averaged over the high latitudes is about 1 ºC colder in EC-Earth-PISM than in EC-Earth (Fig. 6).

Geographically the most substantial warming relative to the piControl in Abrupt4xCO2 occurs over the Arctic Ocean, especially over the central Arctic Ocean and the Greenland Sea (Fig. 7d and e). Across the GrIS, the temperature differences between the EC-Earth and EC-Earth-PISM experiments are predominantly determined by the contrasts in the initial ice sheet elevation. However, the difference in near-surface temperature change in the EC-Earth and EC-Earth-PISM experiments is significant over most of the central Arctic Ocean and around the Faroe Islands (hatched areas in Fig. 7f), and the reduced warming caused by the interactive ice sheet affects most of the Arctic region (Fig. 7f). North of 55 ºN, the atmospheric winter temperature is colder up to a height of 400 hPa in EC-Earth-PISM, and the Arctic stratosphere is up to 1.5 °C warmer (Fig. S3).

3.3 Arctic sea ice

The changes in Arctic near-surface air temperature are strongly related to the changes in sea ice. In winter, sea ice limits the heat exchange between the cold Arctic atmosphere and the much warmer ocean. Figure 8a shows the evolution of the annual mean northern hemispheric sea ice area, Fig. 8b its seasonal variation, and Fig. 8c and d show the mean sea ice extent in March for the last 50 years of the experiments (301–350). In piControl, the sea ice area is stable and evolves very similarly in EC-Earth and EC-Earth-PISM; it reaches its maximum in February, and the minimum occurs in August.

Fig. 8
figure 8

a Time series of mean annual NH sea ice area for all experiments. The thin lines indicate the annual means and thick lines the 11-year running means. b Mean monthly NH sea ice area (mean and one standard deviation) for years 301–350 (grey shaded period highlighted in Fig. 6). Mean March NH sea ice extent for the same 50-year period for c EC-Earth and d EC-Earth-PISM. Red indicates sea-ice coverage in both piControl and Abrupt4xCO2 while blue areas have sea-ice only in piControl

Full size image

In Abrupt4xCO2, there is a rapid reduction in sea ice area during the first years of simulation in both EC-Earth and EC-Earth-PISM (Fig. 8a). The sea ice reduction occurs gradually in 1pctCO2 but after 200 years, the sea ice area remains small in both warming experiments, and the Arctic is seasonally ice-free. However, the colder Arctic and the more freshwater into the Arctic in EC-Earth-PISM do influence the existence of sea-ice in the warm scenarios: the two EC-Earth-PISM experiments have about two times more sea ice in the winter months than the EC-Earth experiments and the sea ice recovers significantly faster in winter (Fig. 8b). Sea ice still exists across the Central Arctic from the Beaufort Sea to the Laptev Sea and the East Siberian Sea, while these areas remain ice-free all year round in EC-Earth. The sea ice in the Baffin Bay also extends further south in EC-Earth-PISM than in EC-Earth (Fig. 8c and d).

3.4 The oceanic response

Sea surface temperature (SST) and sea surface salinity (SSS) for the piControl and Abrupt4xCO2 experiments are shown in Figs. 9 and 10. In piControl, both SST and SSS are comparable in EC-Earth and EC-Earth-PISM with only small differences (Figs. 9c and 10c).

Fig. 9
figure 9

Left: Sea surface temperature (SST, in °C) in the a EC-Earth-PISM and b EC-Earth piControl experiments and c the difference between (a) and (b). Right: Change in SST in Abrupt4xCO2 compared to piControl in d EC-Earth-PISM and e EC-Earth and f the difference between (d) and (e). All values are averaged over the years 301–350. In each subpanel, the hatching marks regions where the difference exceeds one standard deviation of the respective reference simulation. The boxes indicate the four regions used to estimate Deep Mixed Volume in Fig. S5 (starting from bottom left): Labrador Sea, Iceland-Scotland, the Greenland-Iceland-Norway Seas (GIN), and an area north of Barents Sea (Nansen)

Full size image
Fig. 10
figure 10

Left: Sea surface salinity (SSS, in PSU) in the a EC-Earth-PISM and b EC-Earth piControl experiments and c the difference between (a) and (b). Right: Change in SSS in Abrupt4xCO2 compared to piControl in d EC-Earth-PISM, e EC-Earth and f the difference between (d) and (e). All values are averaged over the years 301–350. In each subpanel, the hatching marks regions where the difference exceeds one standard deviation of the respective reference simulation. Note that the salinity in the Caspian Sea is piling up due to unbalanced precipitation and evaporation in the region, which is a closed sea in the ocean model. As the EC-Earth-PISM initial state has gone through a longer spin-up, EC-Earth-PISM has a higher salinity in this region than EC-Earth

Full size image

In Abrupt4xCO2, SST increases in the whole Northern Hemisphere. The most extensive warming (locally up to 12 K) occurs in a region covering the Barents and Greenland Seas (9d and e). This region is located close to the sea ice edge (Fig. 8c and d) and also has an increase in salinity (Fig. 10d and e). Regions initially covered with sea ice, experience the highest warming, as the darker ocean surface absorbs more radiation, which drives the warming and locally amplifies sea ice melting. The strongest warming occurs in winter, close to the sea ice edge, while in summer, the melting of sea ice keeps the SST close to the freezing point (not shown).

Overall, the northern hemispheric ocean surface warms less in EC-Earth-PISM than in EC-Earth, especially in the Arctic Ocean and around Greenland (Fig. 9f). The surface waters are also fresher in EC-Earth-PISM (Fig. 10f), especially around Greenland. The colder and fresher surface water in the Arctic Ocean in Abrupt4xCO2 is a result of the larger freshwater flux from GrIS melt. The Arctic Ocean stays colder and fresher down to a depth of at least 2000 m in EC-Earth-PISM compared to EC-Earth (Fig. S4).

The change in SST and SSS observed in the warming experiments may affect the North Atlantic deep water formation and thereby influence the thermohaline circulation. The ventilation of mixed layer volume below a certain threshold is a proxy of deep convection. Here, the definition of Deep Mixed Volume (DMV) in Brodeau and Koenigk (2016) is used to assess the deep water formation strength. The mean March mixed layer volume below a basin-specific threshold is analyzed in four regions (see Fig. 9a for the locations of regions): the Labrador Sea, Iceland-Scotland (including the Irminger Sea, the Iceland basin, and an area west of Scotland), the Greenland-Iceland-Norway Seas (GIN), and an area north of the Barents and Kara Seas (Nansen). The depth threshold is 1000 m for the first two regions and 750 m for the latter two (Brodeau and Koenigk 2016).

In the piControl experiments, deep convection occurs in the Labrador and GIN Seas, and, to a smaller extent, in the Iceland-Scotland region (Fig S5). In all basins, deep convection is weaker in EC-Earth-PISM because the surface water is fresher. There is no deep convection in the Nansen region. In Abrupt4xCO2, deep convection disappears almost immediately from the Labrador and GIN Seas; instead, it shifts into the northern Nansen region (Fig. S5). This result agrees with the northward movement of deep water formation reported in Brodeau and Konigk (2016) for the RCP8.5 scenario. The deep water formation sites follow the sea ice edge rather closely, shifting the deep water formation towards the North Pole and the Kara Sea. The deep convection is weaker in the experiments with the interactive ice sheet (Fig. S5).

Along with the northward shift of the DMV region in the warmer climate, the North Atlantic meridional overturning circulation (AMOC) slows down. Figure 11 shows the AMOC strength, defined as the maximum stream function at 30° N, in the different experiments. In piControl, the AMOC is slightly weaker in EC-Earth-PISM (15.6 ± 0.8 Sv) than in EC-Earth (16.3 ± 0.9 Sv) and shows a reduced decadal variability. In the warming experiments, the AMOC strength reduces and reaches a minimum approximately 50 years after the maximum CO2 concentration is reached (year 50 in the Abrupt4xCO2 experiment and year 190 in the 1pctCO2 experiment). After the initial reduction, the AMOC recovers gradually. In the first 100 years of the experiments, the AMOC behaves similarly in EC-Earth and EC-Earth-PISM, but the recovery is slower in EC-Earth-PISM; for years 301–350, the average strength is 12.1 ± 1 and 12.2 ± 2 for the Abrupt4xCO2 and 1pctCO2 experiments, respectively, compared to 13.8 ± 0.8 and 13.0 ± 0.8 in EC-Earth for the Abrupt4xCO2 and 1pctCO2 experiments, respectively. The increased amount of freshwater input to the ocean from the GrIS melt in the coupled system reduces the deep water formation, delays the recovery of the AMOC and eventually results in a reduced warming in the Arctic and North Atlantic regions.

Fig. 11
figure 11

Time series of the Atlantic Meridional Overturning Circulation (AMOC) maximum strength (in Sv) at 30° N for all experiments. The thin lines indicate the annual means and the thick lines the 11-year running means. The grey shading illustrates the 50-year period shown in Table 1 and Figs. 5, 7, 8b-d, 9 and 10

Full size image

Reduction of AMOC due to freshwater input from Greenland is also observed in idealized experiments, i.e., common freshwater hosing experiments (e.g., Haskins et al. 2020; Swingedouw et al. 2012, 2014). These experiments are usually performed with a much stronger, unrealistic freshwater input to kick the AMOC response. For example, Swingedouw et al. (2014) applies a fixed hosing rate of 100 mSv (= 0.1 Sv) and find the AMOC strength reduced about 1.1 ± 0.6 Sv, representing a 27 ± 14% of the weakening. In EC-Earth-PISM, we find a similar reduction in the AMOC recovery, even though the additional freshwater release rate from Greenland to the ocean is only about 30 mSv at the end of the Abrupt4xCO2 and 1pctCO2 experiments (Table 1 and Fig. 4b). The lower freshwater flux in the EC-Earth-PISM simulations appear therefore to be sufficient to create the effect of a reduced recovery of the AMOC.

4 Discussion and conclusions

In this study we present the climate model system EC-Earth-PISM, with a 2-way coupled, interactive ice sheet model for Greenland. The climate model simulates ice surface processes and surface mass balance, and the ice sheet model accounts for elevation changes, calving, and basal melt rates. The work demonstrates the first step towards the development of a fully coupled climate-ice sheet model that includes an interactive GrIS.

Multi-century experiments for stable, pre-industrial conditions and two idealized warming scenarios with four times the pre-industrial level of CO2 are performed using this system and compared with similar experiments performed using the AOGCM model EC-Earth. The pre-industrial experiment with EC-Earth-PISM shows a stable GrIS and a climatology comparable to that of the standard EC-Earth model. In the relatively cold pre-industrial climate, the interactive ice sheet does not affect the climate equilibrium. The main differences between the experiments with and without the interactive GrIS are mostly due to differences in the Greenland topography which influences surface temperature and precipitation over Greenland. The GrIS freshwater flux is about 18% larger in EC-Earth-PISM than in EC-Earth but the resulting differences in salinity and temperature are small and patching.

However, inclusion of an interactive GrIS in the climate model may result in significant differences in the response of GrIS and the climate to increases in the CO2 concentration. Under warming climate conditions, in both abrupt and gradual 4xCO2 scenarios, the melt-albedo feedback plays an important role, causing the surface melt to increase by a factor of eight while the calving rate reduces by about 10% as the ice sheet retreats. The increase in the total freshwater flux from the GrIS into the ocean is about 55% larger in the experiments with EC-Earth-PISM and this influences the climate of the Arctic. In comparison to the model system without an interactive GrIS, the EC-Earth-PISM results show lower near-surface air temperatures for the northern high latitudes and particularly over the Arctic Ocean. It also has substantially more Arctic sea ice in winter, and a delayed recovery of a reduced AMOC as shown in Figs. 8 and 11.

The new snow albedo parameterization in EC-Earth-PISM is an important part of the coupling as it allows the melt-albedo feedback to be activated. The surface snow-physics of EC-Earth-PISM could be further improved by adopting a more sophisticated snow albedo parameterization (Helsen et al. 2017) as well as a multi-layer snow model, including compaction of snow and refreezing of meltwater.

Since the SMB in EC-Earth-PISM is calculated at the resolution of the atmospheric model, the ablation zone is not well resolved. Another important limitation of the current coupling approach is that mass-conservation is not guaranteed for grid boxes which are partially ice sheet. This is due to the implementation of the ice mask in EC-Earth and the large differences in resolution between EC-Earth and PISM. The results presented here show that the model underestimates the interactions and feedbacks between the climate and the ice sheet model components in the ablation zone, resulting in a modelled ice mass loss in the 1pctCO2 experiment that is only about one third of that simulated in a model that calculates SMB using elevation classes (Muntjewerf et al 2020b). This suggests the need for applying a more sophisticated downscaling of the SMB to the resolution of the ice sheet model in order to better represent the ablation zone in the climate-ice sheet coupled model system (e.g., Vizcaino et al. 2015; Sellevold et al. 2019; van Kampenhout et al. 2019).

The limited understanding and representation of physical processes, such as calving and ice-ocean interaction, in large scale ice sheet models has been identified as a major source of uncertainties for future projections of the GrIS (Goelzer et al. 2020). In EC-Earth-PISM, the ice sheet model uses a simple geometric calving criterion. The introduction of a physical calving model based on the stress balance divergence estimates within the narrow fjords (e.g., Nick et al. 2013; Benn et al. 2007, 2017; Mercenier et al. 2018) would be needed to better simulate the ice dynamics, which would further require a better resolved atmospheric forcing and an explicit representation of the ocean forcing. Furthermore, a higher horizontal spatial resolution of the ice sheet model would improve the representation of the topographically steered ice streams that transport ice towards the numerous narrow fjord systems surrounding Greenland (e.g., Aschwanden et al. 2016).

Despite the limitations of the current coupling approach, the results presented here stress that an interactive GrIS is important for understanding and assessing climate change in the Arctic. As most GCMs do not include interactive ice sheets, it is likely that these models underestimate the mass loss of the ice sheets under increased levels of CO2 and hence overestimate the northern hemispheric warming. As global warming continues, the melting of the GrIS intensifies, and its influence on the climate in the Arctic as well as the northern hemisphere increases. Climate models with a representation of climate-ice sheet interactions are thus necessary for reliable climate change projections.