Improved Understanding of Multicentury Greenland Ice Sheet Response to Strong Warming in the Coupled CESM2‐CISM2 With Regional Grid Refinement

Abstract

The simulation of ice sheet‐climate interactions, such as surface mass balance fluxes, is sensitive to model grid resolution. Here we simulate the multi‐century evolution of the Greenland Ice Sheet (GrIS) and its interaction with the climate using the Community Earth System Model version 2.2 (CESM2.2) including an interactive GrIS component (the Community Ice Sheet Model v2.1 [CISM2.1]) under an idealized warming scenario (atmospheric $${\text{CO}}_{2}$$ increases by 1% $${\text{yr}}^{-1}$$ until quadrupling the pre‐industrial level and then is held fixed). A variable‐resolution (VR) grid with 1/$$4{}^{\circ}$$ regional refinement over the broader Arctic and $$1{}^{\circ}$$ resolution elsewhere is applied to the atmosphere and land components, and the results are compared with conventional $$1{}^{\circ}$$ lat‐lon grid simulations to investigate the impact of grid refinement. Compared with the $$1{}^{\circ}$$ runs, the VR run features a slower rate of surface melt, especially over the western and northern GrIS, where the ice surface slopes gently toward the periphery. This difference pattern originates primarily from higher snow albedo and, thus, weaker albedo feedback in the VR run. The VR grid better captures the CISM ice sheet topography by reducing elevation discrepancies between CAM and CISM and is, therefore, less reliant on the downscaling algorithm, which is known to underestimate albedo gradients. The sea level rise contribution from the GrIS in the VR run is 53 mm by year 150 and 831 mm by year 350, approximately 40% and 20% less than that of the $$1{}^{\circ}$$ runs, respectively.

Introduction

Recent data reveal an acceleration in mass loss from the Greenland Ice Sheet (GrIS), averaging 257 Gt $${\text{yr}}^{-1}$$ between 2017 and 2020 — a sevenfold increase compared to the early 1990s (Otosaka et al., 2023). GrIS mass loss is driven by both atmospheric warming (Hanna et al., 2021), which increases surface melt and meltwater runoff (Trusel et al., 2018), and oceanic warming, which leads to glacier acceleration and enhanced ice discharge (ID) (Straneo & Heimbach, 2013). Although increased ID played a stronger role in GrIS mass loss between 1992 and 2018 (66 $$\pm $$ 8%), over the past two decades, surface mass balance (SMB) decrease has become the dominant contributor, driven by increased surface melt (Enderlin et al., 2014; Mouginot et al., 2019). Interactions between the ice sheet, atmosphere, and ocean can trigger feedback mechanisms, further amplifying or dampening mass imbalance signals. One important positive feedback is the albedo/melt feedback. As snow or ice melts, surfaces with lower albedo, such as warmer snow, firn, bare ice, or ground, are exposed, leading to increased absorption of shortwave radiation and further accelerating melt in the affected and surrounding regions. Other feedbacks, such as the geometry/SMB feedbacks (Fyke et al., 2018) can enhance or restrain GrIS mass loss, underlining the need for coupled models that capture bidirectional interactions and feedbacks.

The accuracy of simulated SMB is sensitive to model grid resolution, especially in regions with steep and complex terrain. Conventional grid resolutions ($$1{}^{\circ}$$ to $$2{}^{\circ}$$) in global climate models (GCMs) are too coarse to capture the steep topographic gradients in the mountainous Greenland margins. As a result, these models fail to resolve key processes such as orographic precipitation and frequently allow excess moisture intrusion into the ice sheet interior (Pollard & Groups, 2000). Research has shown that higher horizontal resolution allows orographic precipitation to be more accurately resolved, thereby reducing positive precipitation biases (Herrington et al., 2022; van Kampenhout et al., 2019). In addition, the ablation zone around the GrIS margins, where the majority of summer melt occurs, may narrow to tens of kilometers, which cannot be resolved by $$1{}^{\circ}$$ to $$2{}^{\circ}$$ grids. It is therefore crucial to use finer resolution grids for accurately representing GrIS SMB processes.

Modeling with a variable-resolution grid offers several key advantages. Though rapidly advancing, widespread use of global-uniform high-resolution climate models (e.g., models participating in the High-Resolution Model Intercomparison Project (HighResMIP; Haarsma et al., 2016)) remains impractical due to current computational resource limitations. Regional climate models (RCMs), typically operating in one-way nesting mode, provide high-resolution simulations for specific regions at relatively lower computational cost. However, RCMs require boundary conditions from GCMs or reanalysis, thereby preventing two-way interactions across the boundaries. Moreover, boundary conditions derived from a separate host model can lead to inconsistencies between the host model and the RCM. Variable-resolution modeling addresses some of these challenges by using a unified framework that captures two-way interactions between the regional and large scales with better computational efficiency.

The application of regional grid refinement in GCMs dates back to the early use of stretched grids in the late 1970s (Schmidt, 1977; Staniforth & Mitchell, 1978) and is now implemented in many state-of-the-art GCMs (Golaz et al., 2019; Harris et al., 2016; Sakaguchi et al., 2023; Tang et al., 2023; Zängl et al., 2022). In the Community Earth System Model, version2 (CESM2; Danabasoglu et al., 2020), regional grid refinement is supported by the spectral-element (SE; Lauritzen et al., 2018) dynamical core of the atmospheric component. Studies have demonstrated its consistency in modeling global circulation and climatology (Gettelman et al., 2018; Zarzycki et al., 2015), fidelity in representing tropical and extra-tropical cyclones (Zarzycki, 2016; Zarzycki & Jablonowski, 2014; Zarzycki et al., 2014) and regional climate, especially in mountainous or steep terrain regions (Bambach et al., 2022; Huang & Ullrich, 2017; Huang et al., 2016; Rahimi et al., 2019; Rhoades et al., 2016, 2018; Wijngaard et al., 2023; Wu et al., 2017). The variable-resolution CESM2 (VR-CESM2) has also been applied to the polar regions. van Kampenhout et al. (2019) demonstrated significant improvement in simulating GrIS SMB in the accumulation zone by using two regionally refined grids over the GrIS at 1/$$2{}^{\circ}$$ and 1/$$4{}^{\circ}$$. In addition to improvements over the GrIS, the simulated clouds and precipitation in the Arctic are significantly improved with two Arctic-refined meshes, one at 1/$$4{}^{\circ}$$ and another with an additional 1/$$8{}^{\circ}$$ patch of refinement over Greenland (Herrington et al., 2022). For the Antarctic, 1/$$4{}^{\circ}$$ regional refinement over the Antarctic Ice Sheet and the surrounding Southern Ocean shows improvements, primarily in temperature and wind fields, and some degradations related to surface melt over the ice sheet, compared to $$1{}^{\circ}$$ CESM2 (Datta et al., 2023). The VR-CESM2 in the studies mentioned above was run in coupled land-atmosphere mode following the Atmospheric Model Intercomparison Project protocols (AMIP; Gates, 1992).

This study analyzes results from a set of simulations using the fully coupled configuration of CESM2 with a dynamic GrIS under an idealized strong warming scenario. A variable-resolution grid featuring 1/$$4{}^{\circ}$$ regional refinement over the broader Arctic region and $$1{}^{\circ}$$ horizontal resolution elsewhere is used in the atmosphere and land components of CESM2. Unlike prior VR-CESM2 studies, we incorporate coupling to a dynamic ocean model, similar to Tang et al. (2023). This work aims to: first, investigate the multicentury sensitivity of GrIS evolution to a changing climate, and second, compare the variable-resolution run with global $$1{}^{\circ}$$ resolution runs, to assess the added value of regional refinement. Section 2 documents the model, grids, experimental design, and analytical methods. Section 3 presents the results and simulation comparisons. Finally, in Section 4, the discussion and conclusions are provided.

Methods

Model Description

CESM2 is an Earth System Model (ESM) maintained by the National Center for Atmospheric Research (NCAR), which consists of atmosphere, ocean, land, sea ice, and land ice components and can be run in configurations of varying complexity. The ocean component, Parallel Ocean Program version 2 (POP2; Smith et al., 2010), runs on a nominal $$1{}^{\circ}$$ displaced-pole grid with 60 vertical levels. Sea ice is represented by the Community Ice CodE for sea ice version 5 (CICE5; Hunke et al., 2015), which uses the same horizontal grid as POP2. Land processes are simulated by the Community Land Model version 5 (CLM5; Lawrence et al., 2019), sharing the horizontal grid with the atmosphere model. CLM5 also embeds the Model for Scale Adaptive River Transport (MOSART; Li et al., 2013) to handle land surface runoff based on topographic gradients.

The GrIS is simulated using the Community Ice Sheet Model, version 2.1 (CISM2.1; Lipscomb et al., 2019), on a 4-km rectangular grid with 11 terrain-following vertical levels. To simulate ice flow, a depth-integrated higher-order approximation (Goldberg, 2011) of the Stokes equations is employed in the velocity solver. The basal sliding parameterization utilizes a pseudo-plastic sliding law and a simple basal hydrology model, following the approach described by Aschwanden et al. (2016). In this parameterization, the yield stress is determined by the till friction angle and the effective pressure, where the former is influenced by bedrock elevation through a fixed piecewise linear relationship. Bedrock evolution due to Glacial Isostatic Adjustment is governed by an Elastic plate Lithosphere plus Relaxing Asthenosphere (ELRA) model (e.g., Rutt et al., 2009). This study accounts for calving processes through a flotation criterion, where floating ice is discharged immediately to the ocean.

Two versions of CESM2 are used for the simulations in this study: CESM2.1 and CESM2.2. In CESM2.1, the atmosphere is simulated with the Community Atmosphere Model version 6 (CAM6; Gettelman et al., 2019), using the Finite-Volume (FV; Lin, 2004) dynamical core, with 32 vertical hybrid pressure-sigma levels. The CAM6 physical parameterization package is described in detail in Gettelman et al. (2019). CESM2.1 is one of the ESMs that contribute to the Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al., 2016) and the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6; Nowicki et al., 2016).

CESM2.2 uses the same CAM physics parameterizations and vertical grid but contains enhanced functionality for the SE dynamical core, including the capability for running VR grids (Herrington et al., 2022). The CMIP6 CESM2.1 simulations using CAM-FV are not reproducible in CESM2.2 due to code base updates, including regular maintenance and minor scientific changes that do not cause large differences. Therefore, two model versions are required to compare the VR grid with the CMIP6 workhorse configuration.

Surface Mass Balance

The GrIS SMB simulated in CESM2 is the sum of ice accumulation and ice ablation. The SMB processes are calculated in CLM5, which includes up to 10 vertical snowpack layers with a maximum total depth of 10-m water equivalent. Only snow accumulated over the 10-m threshold contributes to ice accumulation. Ice ablation incorporates surface ice melt as well as sublimation. Meltwater, together with rain, can penetrate the snow layers and refreeze, providing an additional source of ice. Liquid water that reaches the snow-ice interface or exceeds the holding capacity of the snowpack runs off to the ocean. Melt energy is calculated as the sum of net surface radiation, latent and sensible turbulent heat fluxes, and ground heat fluxes at the atmosphere-snow interface. To account for sub-grid variability, each glaciated grid cell in CLM5 is subdivided into 10 elevation classes (ECs) with predefined elevation ranges (Lipscomb et al., 2013; Sellevold et al., 2019). The area fractions of the ECs are calculated from the higher-resolution CISM topography. For each EC, surface energy fluxes and SMB are calculated independently by downscaling atmospheric variables. Near-surface temperature is downscaled with a fixed lapse rate of 6 K $${\text{km}}^{-1}$$. Relative humidity is assumed uniform vertically. Due to a CAM6 model bias leading to excessive rainfall over the GrIS, precipitation is repartitioned based on near-surface temperature thresholds: it falls as snow when the temperature is below −$$2{}^{\circ}$$C, as rain when the temperature is above $$0{}^{\circ}$$C, and as a linear combination of snow and rain for temperatures between −$$2{}^{\circ}$$C and $$0{}^{\circ}$$C.

Despite biases such as overestimated precipitation, CESM2.1 at $$1{}^{\circ}$$ resolution with a fixed GrIS geometry simulates a realistic historical GrIS SMB (van Kampenhout et al., 2020). When it is coupled to CISM2, higher GrIS SMB and interannual variability are reported, partially attributed to differences in ice sheet topography (Muntjewerf, Petrini, et al., 2020).

Coupling Scheme

The GrIS is interactively coupled to the other Earth system components within the model framework. CISM2 receives the CLM5 SMB for each EC, which is downscaled by the coupler using a trilinear remapping scheme (bilinear horizontally and linear vertically), with corrections of accumulation and ablation to conserve water mass during downscaling. As the ice sheet evolves, the coupler updates the EC fractional glacier coverage in each CLM5 grid cell based on the CISM2 ice sheet extent on an annual basis. The mean surface elevation in CAM6 is manually updated based on the CISM2 topography every 20 model years, using the CESM topography software (Lauritzen et al., 2015). Surface runoff from CLM5, along with basal melt and ID from CISM2, constitute the freshwater fluxes inputted into the ocean, which are supplied as salinity anomalies. Thus, the coupling is unidirectional, from ice to ocean, since floating ice is not permitted. A more detailed description of the coupling scheme can be found in Muntjewerf et al. (2021).

Grids

The variable-resolution grid used in this study, referred to as the ARCTIC grid (Figure 1a), is a $$1{}^{\circ}$$ SE grid with 1/$$4{}^{\circ}$$ regional refinement over the broader Arctic region (Herrington et al., 2022). It is generated using the software package SQuadgen (). The global $$1{}^{\circ}$$ resolution runs use the latitude-longitude $$1{}^{\circ}$$ grid, referred to as f09, supported by the FV dynamical core. Figures 1b and 1c present a snapshot of the surface topography over the GrIS before the start of the warming scenario for the f09 and ARCTIC grids, respectively. The ARCTIC grid represents more detailed Greenland topography (e.g., the south dome and ice sheet periphery) with a more accurate ice sheet mask. Surface elevation differences between the two grids can exceed 700 m (Figure 1d), partly due to the initial ice sheet volume differences. The physics time step of the ARCTIC simulations is set to 450 s, a 4$${\times} $$ reduction relative to the default 1800 s time step of the $$1{}^{\circ}$$ grid, to avoid large time truncation errors (Herrington & Reed, 2018).

[IMAGE OMITTED. SEE PDF]

Several reasons motivated the choice of f09 as the reference grid instead of the quasi-uniform $$1{}^{\circ}$$ SE grid (ne30). First, the coarser ne30 GrIS exacerbates biases in melt and precipitation relative to f09. Previous studies have demonstrated that f09 outperforms ne30 in representing historical GrIS SMB (Herrington et al., 2022). Hence, the f09 grid model, the CESM2 CMIP6 workhorse, presents a more challenging benchmark for the ARCTIC grid model to outperform than the ne30 grid. Furthermore, as mentioned in the following subsection, a pre-industrial simulation is required to bring the GrIS back to a near-equilibrium state after a grid change from the long spun-up state in the f09 grid. This process results in different initial ice sheet volumes at the start of the warming scenario for the SE grid simulations. This inevitable discrepancy also applies when comparing simulations using the ARCTIC grid and the f09 grid. Concurrent differences in the dynamical core and physics time step further complicate this comparison. Yet, we will explain in Section 4 that the differences between the simulations are primarily driven by varying the horizontal resolution rather than the other factors.

Experiment Design

First, a pre-industrial simulation was branched off from the last segment of the spun-up pre-industrial Earth system/ice sheet state (Lofverstrom et al., 2020). A series of experiments was conducted to bring the top of atmosphere radiative forcing into balance, using common tuning parameters such as the width of the sub-grid distribution of vertical velocity $$\mathtt{c}\mathtt{l}\mathtt{u}\mathtt{b}\mathtt{b}\mathtt{\_}\mathtt{g}\mathtt{a}\mathtt{m}\mathtt{m}\mathtt{a}$$ (Guo et al., 2015). Due to computational limitations, the tuned pre-industrial control was run for 180 years until the GrIS achieved a near-equilibrium state with a small residual drift (Table 1). An idealized warming scenario was then initialized, with atmospheric $${\text{CO}}_{2}$$ concentration increasing by 1% per year until reaching 4$${\times} $$ the pre-industrial value after 140 years, followed by a 210-year simulation with fixed 4$${\times} $$ pre-industrial $${\text{CO}}_{2}$$ concentration (Figure 2a). This simulation using the ARCTIC grid (hereafter ARCTIC) was compared to two simulations under the same forcing but on the f09 grid, run in the older CESM2.1 code base. One is from Muntjewerf, Sellevold, et al. (2020) (hereafter F09M). In this code base, to reduce the excessively high SMB over portions of the GrIS in CESM2 coupled runs, cold rain (<−$$2{}^{\circ}$$C) produced by CAM immediately runs off to the ocean instead of being converted to snow by CLM. Additionally, the sub-grid roughness over Greenland is artificially increased to facilitate low-level convergence and precipitation near the coasts. To limit the differences between models, another f09 simulation was run using CESM2.1 but without these two adjustments (hereafter F09).

Table 1 Annual Rate of Mass Loss (mm $${\text{yr}}^{\mathit{-}\mathit{1}}$$), Cumulative Greenland Ice Sheet (GrIS) Mass Loss (mm), Mass Balance Components (Gt $${\text{yr}}^{\mathit{-}\mathit{1}}$$), GrIS Area ($${\mathit{10}}^{\mathit{6}}$$ $${\text{km}}^{\mathit{2}}$$), and GrIS Volume ($${\mathit{10}}^{\mathit{6}}$$ $${\text{km}}^{\mathit{3}}$$)

	Last 20 years of CTRL	Years 131–150	Years 231–250	Years 331–350
$$\mathtt{F}\mathtt{09}\mathtt{M}$$	$$\mathtt{F}\mathtt{09}$$	$$\mathtt{A}\mathtt{R}\mathtt{C}\mathtt{T}\mathtt{I}\mathtt{C}$$	$$\mathtt{F}\mathtt{09}\mathtt{M}$$	$$\mathtt{F}\mathtt{09}$$	$$\mathtt{A}\mathtt{R}\mathtt{C}\mathtt{T}\mathtt{I}\mathtt{C}$$	$$\mathtt{F}\mathtt{09}\mathtt{M}$$	$$\mathtt{F}\mathtt{09}$$	$$\mathtt{A}\mathtt{R}\mathtt{C}\mathtt{T}\mathtt{I}\mathtt{C}$$	$$\mathtt{F}\mathtt{09}\mathtt{M}$$	$$\mathtt{F}\mathtt{09}$$	$$\mathtt{A}\mathtt{R}\mathtt{C}\mathtt{T}\mathtt{I}\mathtt{C}$$
Annual mass loss	−0.04	−0.10	−0.05	2.08	2.06	1.48	5.28	4.49	3.50	6.36	5.93	5.40
Cumulative mass loss	−0.8	−2.0	−1.1	97	84	53	501	447	344	1,098	976	831
MB	19	41	23	−776	−761	−542	−1,974	−1,669	−1,285	−2,376	−2,195	−2,001
SMB	616	723	685	−380	−319	−72	−1,797	−1,463	−1,081	−2,284	−2,097	−1,909
ID	573	654	636	376	420	448	161	187	187	78	81	77
BMB	−24	−27	−25	−20	−22	−22	−16	−18	−18	−14	−16	−16
GrIS area	1.97	2.00	2.02	1.92	1.96	1.99	1.77	1.80	1.83	1.60	1.64	1.66
GrIS volume	3.23	3.27	3.25	3.20	3.24	3.23	3.05	3.10	3.12	2.81	2.89	2.93

[IMAGE OMITTED. SEE PDF]

Analysis

Atmospheric and Oceanic Circulation Metrics

The Greenland blocking index (GBI) uses the 500-hPa geopotential height ($${\mathrm{Z}}_{500}$$) to estimate blocking over the Greenland region (Fang, 2004). Strong and persistent blocking can result in extreme summer melt at the ice sheet surface (Hanna et al., 2014). The revised GBI from Hanna et al. (2018) is used, calculated by subtracting the area-averaged $${\mathrm{Z}}_{500}$$ over the Arctic region (6$$0{}^{\circ}$$N to 8$$0{}^{\circ}$$N) from the area-averaged $${\mathrm{Z}}_{500}$$ over the Greenland region (6$$0{}^{\circ}$$N to 8$$0{}^{\circ}$$N, and 8$$0{}^{\circ}$$W to 2$$0{}^{\circ}$$W). The resulting time series is then standardized relative to the last 80 years of the pre-industrial period. Here only the JJA mean GBI is considered.

The North Atlantic Meridional Overturning Circulation (NAMOC) index measures the strength of the NAMOC, which is predicted to weaken with the addition of GrIS meltwater to the ocean (Muntjewerf, Sellevold, et al., 2020; Vizcaíno et al., 2010). The NAMOC index is defined as the maximum of the overturning stream function north of 2$$8{}^{\circ}$$N and below 500-m depth.

Melt/Albedo Feedback

We use the melt/albedo feedback (or albedo feedback) calculations following Box et al. (2012). The albedo feedback $$\left({\alpha }_{\mathit{feedback,a}}\right)$$ is quantified by a regression of 20 annual samples of detrended anomalies of summer (JJA) average net shortwave radiation $$\left(S{W}_{\mathit{net}}\right)$$ and near-surface air temperature $$\left({T}_{\mathit{air}}\right)$$, in units of W $${\mathrm{m}}^{-2}$$ $${\mathrm{K}}^{-1}$$, with anomalies indicated by the ^′ character in $${{\Delta }}^{\prime }$$: 1 $${\alpha }_{\mathit{feedback,a}}={{\Delta }}^{\prime }S{W}_{\mathit{net}}/{{\Delta }}^{\prime }{T}_{\mathit{air}}$$

The regression uses annual pairs of anomalies instead of successive values in the time series. This method of pairing is illustrated in Figure S1 in Supporting Information S1. Since this definition of albedo feedback excludes lags that may influence albedo changes (e.g., a low-albedo year preconditioning the next year for low albedo), an alternative formulation of albedo feedback, referred to as the bulk albedo feedback $$\left({\alpha }_{\mathit{feedback,b}}\right)$$, is considered, which is the change in $$S{W}_{\mathit{net}}$$ over the change of $${T}_{\mathit{air}}$$: 2 $${\alpha }_{\mathit{feedback,b}}={\Delta }S{W}_{\mathit{net}}/{\Delta }{T}_{\mathit{air}}$$

Equilibrium Line Altitude

The equilibrium line is characterized by a zero annual mean SMB, which separates the ice sheet accumulation zone from the ablation zone. To calculate the average Equilibrium Line Altitude (ELA) of the GrIS in our model, we use the following algorithm: Iterate over the grid cells in the ablation zone (ablation grid) using the annual mean SMB field. If any of the neighboring grid cells have positive SMB, linear-interpolate the elevation where SMB equals zero based on the SMB and elevation values of the positive SMB grid (with grid area $${a}_{p}$$) and the ablation grid (with grid area $${a}_{n}$$). Record this computed elevation as one ELA value, with the approximate length of the shared edge of the two grid cells as its weight. The approximate edge length is calculated as $$\left(\sqrt{{a}_{n}}+\sqrt{{a}_{p}}\right)/2$$. The final average ELA is the length-weighted average of all the recorded ELA values.

Results

The results have been structured as follows. Section 3.1 describes the response of the GrIS in ARCTIC, with the driving factors analyzed in Section 3.2. Then Section 3.3 compares ARCTIC with the two $$1{}^{\circ}$$ runs. To illustrate spatial change patterns over time and analyze the differences between F09M, F09, and ARCTIC, this analysis focuses on three time periods: years 131–150, 231–250, and 331–350. Years 131–150 represent the $${\text{CO}}_{2}$$ stabilization period. Years 231–250 and 331–350 occur one and two centuries later, respectively, with the latter marking the end of the simulation.

GrIS Response in the VR Run

In general, our findings are similar to those in Muntjewerf, Sellevold, et al. (2020). The differences are only in magnitude and timing, which will be analyzed in detail in Section 3.3. A brief overview of the response of some key variables is given below.

The 1% $${\text{yr}}^{-1}$$ $${\text{CO}}_{2}$$ increase causes a nearly linear rise in global average near-surface temperature (0.3 K per decade) with 80% greater warming for the Arctic and 10% greater warming for the GrIS (Figure 2b). By $${\text{CO}}_{2}$$ stabilization, the GrIS has warmed by 5.0 K, with an additional 3.8 K of warming by the end of the simulation.

Consistent with Muntjewerf, Sellevold, et al. (2020), our simulation reveals a strong increase in mass loss acceleration after $${\sim} $$110 years, rising from 2.4 Gt $${\text{yr}}^{-2}$$ before year 110 to 13.0 Gt $${\text{yr}}^{-2}$$ by year 150, and then decreasing gradually. Over 350 years, the ice sheet area shrinks by 17% (Figure 2c). The cumulative contribution to global mean sea level rise reaches 53 mm by year 150 and 831 mm by year 350 (Table 1). The MB trend is dominated by the decreasing SMB, which becomes negative around year 130 (Figure 2d). Ice discharge gradually decreases as marine-terminating outlet glaciers thin, decelerate, and even transition to land-terminating glaciers.

The decline in SMB, dominated by enhanced surface melt (Figure 2e), aligns with earlier findings (Muntjewerf, Sellevold, et al., 2020). During years 131–150 and 331–350, surface melt exceeds four and eight times the pre-industrial value, respectively (Table 2). While total precipitation initially increases due to increased rainfall, the trend reverses during the last century as decreasing snowfall dominates. By the end of the simulation, total precipitation is $${\sim} $$25% higher than the pre-industrial value (Table 2). Refreezing increases before $${\text{CO}}_{2}$$ stabilization due to more available liquid water from increased rainfall and surface melt, and then stabilizes as snow cover becomes saturated. Sublimation is relatively small throughout the simulation and is therefore not discussed further.

Table 2 Annual Ice Sheet-Integrated Surface Mass Balance and Components Mean (Gt $${\text{yr}}^{\mathit{-}\mathit{1}}$$)

	Last 20 years of CTRL	Years 131–150	Years 231–250	Years 331–350
$$\mathtt{F}\mathtt{09}\mathtt{M}$$	$$\mathtt{F}\mathtt{09}$$	$$\mathtt{A}\mathtt{R}\mathtt{C}\mathtt{T}\mathtt{I}\mathtt{C}$$	$$\mathtt{F}\mathtt{09}\mathtt{M}$$	$$\mathtt{F}\mathtt{09}$$	$$\mathtt{A}\mathtt{R}\mathtt{C}\mathtt{T}\mathtt{I}\mathtt{C}$$	$$\mathtt{F}\mathtt{09}\mathtt{M}$$	$$\mathtt{F}\mathtt{09}$$	$$\mathtt{A}\mathtt{R}\mathtt{C}\mathtt{T}\mathtt{I}\mathtt{C}$$	$$\mathtt{F}\mathtt{09}\mathtt{M}$$	$$\mathtt{F}\mathtt{09}$$	$$\mathtt{A}\mathtt{R}\mathtt{C}\mathtt{T}\mathtt{I}\mathtt{C}$$
SMB	–	701	651	−745	−620	−369	−2,213	−1,752	−1,398	−2,552	−2,254	−2,159
Precipitation	942	1,026	955	1,047	1,273	1,200	1,106	1,374	1,277	1,156	1,308	1,265
Snowfall	850	934	869	782	961	930	699	932	913	695	819	849
Rain	92	92	86	265	312	270	406	443	364	461	488	416
Refreezing	142	295	307	680	872	758	784	956	824	781	1,019	830
Melt	485	476	468	2,147	2,400	1,986	3,662	3,606	3,083	4,009	4,073	3,806
Sublimation	–	52	57	60	53	70	34	34	52	19	19	32

Changes in surface melt and ice dynamics reshape the ice sheet. SMB is greatly reduced over the ice sheet periphery, exceeding 5 m $${\text{yr}}^{-1}$$ along the southeast and west margins by the end of the simulation (Figures 3f–3h). In the ice sheet interior, SMB increases due to the locally enhanced snowfall. This drives significant thinning toward the margins and slight thickening over the high interior (Figures 3j–3l). Ice flow redistributes mass, with surface ice velocity decreasing along the ice sheet margins due to the thinning of outlet glaciers and increasing over the steeper slopes caused by the thickness changes (Figures 3n–3p). For instance, surface ice velocity along the 2,500 m contour line, initially 35 m $${\text{yr}}^{-1}$$ before warming, increases by 4%, 22%, and 52% over the three periods, respectively.

[IMAGE OMITTED. SEE PDF]

Drivers for Melt Change

The drivers of surface melt change have been extensively analyzed by Muntjewerf, Sellevold, et al. (2020). To confirm the dominant role of ablation zone expansion, we expand upon their analysis by subdividing the total integrated surface melt into three components: melt within the original ablation zone, identified using a mask averaged over the last 20 years of the pre-industrial period; melt within the extended ablation zone (the ablation zone of each year minus the pre-industrial ablation zone); and melt within the accumulation zone. As evident from Figure 2e, the accelerated melt increase around year 110 is primarily due to the extended ablation zone. At this time, there is a faster increase of the ablation area (Figure 2c). In contrast, melt within the original ablation zone and accumulation zone remains relatively small and constant over time. Melt in the original ablation zone increases during the first 160 years due to extended melt periods and decreased surface albedo. Then it decreases as the ice sheet shrinks (Figure 2c) and formerly glaciated grid cells gradually become ice-free. Melt in the accumulation zone follows a similar pattern, influenced by expanding melt extent and shrinking accumulation area.

Expansion of the ablation zone initiates a strong melt/albedo feedback. Darker surfaces with lower albedo are exposed (Figures 3b–3d), absorbing more shortwave radiation and further enhances surface melting. As proposed in Muntjewerf, Sellevold, et al. (2020), the quasi-parabolic shape of the ice sheet also accelerates the ablation zone's expansion as it approaches the interior plateau. The rapidly expanded ablation zone also increases turbulent heat fluxes. By the end of the simulation, net shortwave radiation provides 37% of the total additional melt energy (Table S1 in Supporting Information S1), followed by sensible heat flux (23%), net longwave radiation (22%), latent heat flux (10%), and ground heat flux (9%).

To examine the role of large-scale atmospheric circulation on summer GrIS surface melt, we calculated the GBI for ARCTIC (Figure 4a). The GBI shows a negative trend before $${\text{CO}}_{2}$$ stabilization, indicating weakened summer blocking over the Greenland region. This finding aligns with Sellevold and Vizcaíno (2020), who used an AMIP style configuration of CESM2.1 without a dynamic GrIS under the same 1% $${\text{yr}}^{-1}$$ $${\text{CO}}_{2}$$ warming scenario. A more negative GBI is associated with a reduction in surface melt. After $${\text{CO}}_{2}$$ stabilization, the GBI shows no significant trend. Figure 4b shows the linear regression between GrIS-integrated JJA melt and the JJA GBI. Both variables were filtered using a 10-year high-pass filter to focus on sub-decadal timescales. On sub-decadal timescales, the GBI explains 40% of the annual variability of summer surface melt. This suggests that the GBI, and more broadly the atmospheric circulation pattern, is not a driver of melt acceleration but instead counteracts some effects of global warming on surface melt before $${\text{CO}}_{2}$$ stabilization and explains certain interannual melt variability.

[IMAGE OMITTED. SEE PDF]

Impacts of Enhanced Resolution

Large-Scale Climate

We first examine the large-scale climate representations between the regionally refined ARCTIC run and the $$1{}^{\circ}$$ runs. Figures 5a–5c show the summer mean 500 hPa geopotential height for ARCTIC during the three time periods - years 131–150, 231–250, and 331–350. As the atmosphere warms, the 500 hPa geopotential height increases over the northern high latitudes. The differences between ARCTIC and the $$1{}^{\circ}$$ runs show a consistent spatial pattern across all periods (Figures 5d–5i). Compared with F09M and F09, ARCTIC exhibits significantly lower 500 hPa geopotential height over Greenland except during years 331–350 relative to F09M. This lower geopotential indicates weaker blocking and possibly increased cyclonic flows over Greenland in ARCTIC. ARCTIC also shows significantly higher geopotential over subpolar regions compared to the $$1{}^{\circ}$$ runs, with this anomaly extending further north over time.

[IMAGE OMITTED. SEE PDF]

Figure 6 compares summer temperatures over the northern hemisphere between ARCTIC and the $$1{}^{\circ}$$ runs. The lower-troposphere summer virtual temperature, computed by equating a layer mean virtual temperature with the 500–1,000 hPa geopotential thickness, is higher in ARCTIC across much of the northern hemisphere (Figures 6a–6f). This is consistent with the typical response to increasing horizontal resolution (also reducing physics time step) in general circulations models (Pope & Stratton, 2002; Roeckner et al., 2006) including CAM (Herrington & Reed, 2020): increasing the horizontal resolution warms the climate, since higher resolved vertical velocities generate more condensational heating. However, F09M and F09 exhibit a warmer lower troposphere than ARCTIC, centered over the GrIS and extending to the Canadian Archipelago and Alaska during years 131–150 (Figures 6a and 6d) and to East Siberia during years 231–250 (Figures 6b and 6e). Near-surface temperature differences exhibit a distinct spatial pattern compared to the lower troposphere (Figures 6g–6l). Except in some terrestrial regions such as parts of Siberia and Eurasia, ARCTIC is significantly cooler than the $$1{}^{\circ}$$ runs at the near-surface level. This results from the cooler pre-industrial climate in ARCTIC (Table S1 in Supporting Information S1), likely due to tuning in which the albedo of snow over sea ice was increased. Regions with a cooler lower troposphere in ARCTIC during years 131–150 and 231–250 also exhibit much lower near-surface temperature. However, near-surface temperature differences over the GrIS are smaller due to its perennial ice and snow cover. Some regions at the ice sheet periphery are warmer in ARCTIC, possibly due to differences in cloud conditions, which will be discussed in Section 3.3.3. We note that ARCTIC has a cooler baseline climate over the GrIS, and its impact will be discussed in Sections 3.3.3 and 3.3.4.

[IMAGE OMITTED. SEE PDF]

The evolution of NAMOC appears to be insensitive to the increased atmospheric resolution. The three simulations show very similar NAMOC indexes, with even their initial differences diminishing over time (Figure 7a). NAMOC weakens significantly during the period of $${\text{CO}}_{2}$$ increase, with NAMOC index decreasing by $${\sim} $$0.13 Sv $${\text{yr}}^{-1}$$, and then gradually stabilizes, with NAMOC index remaining around 5 Sv for nearly two centuries. The NAMOC index begins to decline at the onset of $${\text{CO}}_{2}$$ increase, well before the rapid rise in freshwater flux around year 110, and then remains stable even as freshwater flux continues to increase (Figure 7b). A similar relationship is observed in simulations under SSP5-8.5 forcing (Muntjewerf, Petrini, et al., 2020), suggesting a relatively limited role of additional freshwater input from the GrIS in NAMOC weakening compared to global warming and non-GrIS freshwater flux sources in CESM2.

[IMAGE OMITTED. SEE PDF]

ARCTIC shows a slower decline in Northern Hemisphere sea ice extent than the $$1{}^{\circ}$$ runs (Figure 7c). In the $$1{}^{\circ}$$ runs, the Arctic becomes ice-free (sea ice extent < one million $${\text{km}}^{2}$$) in September after $${\sim} $$60 years, whereas this occurs roughly five decades later in ARCTIC. All the simulations show an accelerated annual sea ice decline around year 90 (Figure 7c). This acceleration may be due to the ice-free Arctic Ocean in summer absorbing more radiation, which increases ocean heat storage and slows sea ice formation during colder seasons. However, the exact timing of this phenomenon requires further investigation. Before the Arctic becomes ice-free year-round in the last century, ARCTIC maintains a larger sea ice extent than the $$1{}^{\circ}$$ runs due to its slower rate of sea ice decline. This difference is primarily driven by the larger initial sea ice thickness in ARCTIC (Figure 7d), which likely results from increasing the albedo of snow over sea ice in ARCTIC during tuning. The impact of atmosphere resolution on Arctic sea ice conditions is indirect and model dependent (Selivanova et al., 2024), thus further research is needed to investigate the sensitivity of Arctic sea ice simulations to atmospheric model resolution in CESM2.

SMB Evolution and Ice Sheet Changes

The ice sheet evolution in the three simulations follows a similar overall pattern but differs in magnitude and timing. Figure 8b compares the evolution of the GrIS-integrated SMB of the three simulations. In ARCTIC, SMB decreases more slowly than in the $$1{}^{\circ}$$ runs, with the annual mean SMB drop by the end of the simulation being 306 Gt $${\text{yr}}^{-1}$$ and 226 Gt $${\text{yr}}^{-1}$$ ($${\sim} $$10%) smaller than in F09M and F09, respectively (Table 1). The SMB difference is primarily due to reduced surface melt in ARCTIC compared to the $$1{}^{\circ}$$ runs, particularly during years 180–260 (Figure 8a). During years 231–250, GrIS-integrated melt in F09M and F09 exceeds that in ARCTIC by 579 Gt $${\text{yr}}^{-1}$$ and 523 Gt $${\text{yr}}^{-1}$$ ($${\sim} $$18%), respectively (Table 2). Differences in precipitation and refreezing between the simulations are relatively small (Figure 8a). F09 has higher precipitation and refreezing than ARCTIC, partially offsetting the greater melting and reducing SMB differences. In F09M, the adjustments that direct cold-rain to surface runoff and increase sub-grid roughness over Greenland (see Section 2.5) cause a substantial reduction in precipitation, averaging 204 Gt $${\text{yr}}^{-1}$$ less than F09 over 350 years. Refreezing is also reduced in F09M due to the removal of snowfall that could be converted from cold-rain, which reduces the porous space necessary for effective meltwater retention and refreezing. It may also weaken the thermal insulation of snow layers, making it harder to maintain the colder temperatures that favor refreezing. Sublimation differences between the simulations are small (Table 2) and thus not discussed further.

[IMAGE OMITTED. SEE PDF]

The slower SMB decrease in ARCTIC results in a slower MB decrease (Figure 8c), leading to a slower ice volume decline (Figure 8d) and a smaller contribution to global mean sea level rise (Figure 8e). Over the 350 years, F09M and F09 contribute 267 and 145 mm ($${\sim} $$20%) more sea level rise than ARCTIC, respectively (Table 1). We acknowledge that the three simulations start with different initial ice volumes and mass balances, which will be discussed in Section 4.

Averaged over the 350-year period, the smaller melt in ARCTIC is most pronounced over the western and northern gently sloping ice sheet surfaces, where F09M and F09 melt over 300 mm more ice annually (Figures 9e and 9f). In contrast, the southernmost GrIS exhibits more melt in ARCTIC, reversing the larger spatial pattern. The spatial pattern of melt differences between ARCTIC and the $$1{}^{\circ}$$ runs remains consistent despite their differing initial conditions (Figure 8d). The drivers of this melt difference pattern will be analyzed in Section 3.3.3.

[IMAGE OMITTED. SEE PDF]

For total precipitation, distinct patterns are detected, with F09M estimating lower precipitation and F09 estimating higher precipitation than ARCTIC, particularly along the south and southeast coasts (Figures 9a and 9b). The smoother topography in the $$1{}^{\circ}$$ grid allows more moisture to penetrate further into the ice sheet from the southeast. As mentioned above, directing cold-rain to surface runoff and increasing sub-grid roughness in F09M significantly reduces total precipitation. Over the higher plateau where meltwater almost all refreezes, the $$1{}^{\circ}$$ runs exhibit stronger refreezing given their larger surface melt (Figures 9c and 9d). In low-elevation regions, differences in refreezing between ARCTIC and F09 is relatively small. While for F09, the reduced snow accumulation significantly limits refreezing, particularly along the south and southeast coasts.

After 350 years, ARCTIC exhibits weaker ice sheet thinning particularly over the western and northern gently sloping ice sheet surfaces (differences $${ >} $$100 m) compared to the $$1{}^{\circ}$$ runs (Figures 9g and 9h). In the southernmost GrIS, higher precipitation in F09 reinforces lower melt, leading to reduced ice sheet thinning. While for F09M, reduced snowfall and refreezing exceed lower melt values, leading to enhanced thinning in that region.

Surface Energy Balance

To explain the differences in meltwater production between ARCTIC and the $$1{}^{\circ}$$ runs, we examine the energy fluxes that reach the ice sheet surface.

As a key impact factor of surface incident radiation, clouds over the GrIS are significantly altered by the enhanced resolution. Along the coast, the steeper topography in ARCTIC moves clouds more offshore (Figures S2a and S2b in Supporting Information S1), consistent with the findings of Herrington et al. (2022). By contrast, in the interior, likely due to lower pressure over Greenland (Figure 5) that leads to weaker large-scale subsidence and cloud dissipation, higher cloud fractions occur in ARCTIC, especially during years 231–250. In the aggregate, the enhanced cloud cover due to weaker blocking in the interior outweighs any reduction from decreased moisture intrusion. To further evaluate the impact of clouds, we examine cloud liquid and ice water paths, which are more directly related to cloud optical thickness and radiative scattering properties than the simpler cloud fraction metric. We find that the spatial pattern of total cloud liquid water path differences over Greenland largely mirrors that of cloud fraction differences, except for larger differences in southern Greenland (Figures 10a and 10b). The differences in cloud ice water path is relatively small, thus is not shown here.

[IMAGE OMITTED. SEE PDF]

Surface incident radiation is affected by cloud conditions and atmospheric temperature. ARCTIC exhibits higher incident shortwave radiation along the margins, especially in the north (Figures 10c and 10d), which is primarily driven by thinner clouds or reduced cloud cover in these regions (Figures S2c and S2d in Supporting Information S1). Conversely, ARCTIC exhibits lower incident longwave radiation across the ice sheet (Figures 10g and 10h). These lower values are most prominent along the northern and southeastern margins, driven primarily by lower atmospheric temperature (Figure 6; Figures S2g and S2h), with reduced incident longwave radiation due to clouds playing a secondary role (Figures 10e and 10f).

The major differences in surface energy balance (SEB) are driven strongly by differences in albedo. Despite the higher incident shortwave radiation, especially along the ice sheet margins, ARCTIC exhibits reduced net shortwave radiation inland relative to the $$1{}^{\circ}$$ runs, with the largest reduction over the gently sloping ice surface (or transitional area) in the western and northern basins (Figures 11c and 11d). The albedo in these regions is higher in ARCTIC than in the $$1{}^{\circ}$$ runs (Figures S2i and S2j in Supporting Information S1), leading to less shortwave radiation absorption at surface, which becomes pronounced in the last two centuries. In contrast, the weaker net longwave radiation in ARCTIC aligns with the weaker incident longwave radiation (Figures 10g and 10h) and remains consistent through time (Figures 11e and 11f).

[IMAGE OMITTED. SEE PDF]

The simulations also exhibit non-negligible differences in turbulent heat fluxes (Figures 11g–11j). Larger sensible heat fluxes are observed in ARCTIC along the ice sheet margins due to higher near-surface temperatures (Figures 6g–6l), and the differences increase through time. The gently sloping ice surfaces in the western and northern basins exhibit smaller sensible and latent heat fluxes in ARCTIC. Smaller sensible heat fluxes in these regions result from lower near-surface temperatures in ARCTIC (Figures 6g–6l; Figure S3a in Supporting Information S1); reduced latent heat fluxes result from lower atmospheric moisture (Figures S3b and S3c in Supporting Information S1). The spatial distribution of summer evaporation differences reveals weaker evaporation in ARCTIC over the North Atlantic and North America (Figures S4e and S4f in Supporting Information S1), which have been recognized as major moisture sources for the GrIS (Nusbaumer et al., 2019). Despite significant differences in sea ice extent (Figure 7c and Figures S4b and S4c in Supporting Information S1), evaporation differences over sea ice regions remain relatively small (Figures S4e and S4f in Supporting Information S1).

To evaluate the aggregate contribution of each surface energy flux to the summer melt energy differences, the area-weighted mean fluxes over the GrIS were calculated. Figures 12a and 12b display the mean energy flux differences between ARCTIC and the $$1{}^{\circ}$$ runs. Before warming starts, summer melt energy differences between ARCTIC and the $$1{}^{\circ}$$ runs are near zero. However, the contributions of SEB components are substantially different. In ARCTIC, reduced net longwave radiation due to the cooler baseline climate is offset by stronger net shortwave radiation and sensible heat fluxes. Therefore, to identify the dominant energy driver for melt energy discrepancies during warming, initial energy flux differences must be removed.

[IMAGE OMITTED. SEE PDF]

With initial differences removed, it is evident that the reduced net shortwave radiation in ARCTIC dominates in lowering melt energy (Figures 12c and 12d; Figure S6 in Supporting Information S1). In comparison to the net shortwave radiation values reported by ARCTIC (blue solid line; Figures 12c and 12d), we calculated net shortwave radiation using incident shortwave radiation from ARCTIC, but surface albedo fields from the $$1{}^{\circ}$$ runs (blue dashed lines; Figures 12c and 12d). The comparison shows that absorbed shortwave radiation in ARCTIC would rise substantially with $$1{}^{\circ}$$ albedo fields. In comparison, negative net longwave radiation anomalies in ARCTIC are smaller with the initial differences removed, given their relatively consistent pattern over time. Stronger negative latent heat and ground heat flux anomalies further reduce melt energy in ARCTIC. The negative ground heat flux anomalies are primarily due to reduced refreezing over the plateau as a result of less melting in ARCTIC. Within the last century, increasingly positive sensible heat flux anomalies, together with reduced negative latent heat flux anomalies, narrow melt energy differences between ARCTIC and the $$1{}^{\circ}$$ runs. Figures 12e and 12f show the cumulative energy flux anomalies. This analysis confirms that reduced melt energy in ARCTIC is primarily driven by higher albedo values weakening net shortwave radiation.

Melt/Albedo Feedback and the Impact of Ice Sheet Hypsometry

The higher albedo in ARCTIC is a result of weaker albedo feedback. We first examine the bulk albedo feedback (Equation 2), which captures the accumulated change in net shortwave radiation relative to near-surface temperature, thereby accounting for time-lagged effects such as melt preconditioning. Over the 350 years, the bulk albedo feedback in ARCTIC (Figure 13a) and the $$1{}^{\circ}$$ runs (not shown) is positive across the GrIS, indicating increased shortwave radiation absorption as near-surface temperature rises. The positive bulk albedo feedback is most pronounced in the ablation zones, with peaks at lower elevations in the western and northern basins ($${ >} $$15 W $${\mathrm{m}}^{-2}$$ $${\mathrm{K}}^{-1}$$). Compared to F09M and F09, ARCTIC exhibits weaker bulk albedo feedback (Figures 13b and 13c), indicating that less shortwave radiation is absorbed for the same near-surface temperature increase. The spatial patterns of bulk albedo feedback differences align with those of net shortwave radiation differences (Figures 11c and 11d), suggesting that this pattern is primarily driven by differences in net shortwave radiation rather than near-surface temperature. Although the bulk albedo feedback peaks at lower elevations, the differences between ARCTIC and the $$1{}^{\circ}$$ runs are larger over the higher ablation zones near the equilibrium line. We also examine the albedo feedback (Equation 1) using annual paired detrended anomalies of net shortwave radiation and near-surface temperature over the three 20-year periods. Compared to bulk albedo feedback, albedo feedback reflects the interplay of physical mechanisms but is more sensitive to interannual variability, such as fluctuations in snowfall, making it harder to interpret. During years 231–250 and 331–350, when albedo feedback is strong (Figures S7b and S7c in Supporting Information S1), its differences exhibit a similar but more variable spatial pattern compared to the bulk albedo feedback (Figures S7e, S7f, S7h, and S7i in Supporting Information S1). The differences between bulk albedo feedback and albedo feedback highlight the role of melt preconditioning subsequent melt seasons.

[IMAGE OMITTED. SEE PDF]

The albedo feedback parameter presumes linearity over the range of temperature changes it is applied to, however albedo is a non-linear function of absolute temperature. Therefore it is possible that the $$1{}^{\circ}$$ runs have a larger albedo feedback parameter because their baseline temperatures are warmer, and the albedo feedback begins earlier on in the simulations and remains engaged over a longer duration of the experiment than the ARCTIC run. To isolate the impact of different baseline climate, we normalized the summer mean GrIS surface albedo (Figure 14a) by both near-surface air temperature (Figure 14b) and lower troposphere virtual temperature (Figure 14c). The differences in albedo are small between the simulations for cooler temperatures, but ARCTIC exhibits higher albedo values for warmer temperatures, with this increase more pronounced for albedo conditioned on lower troposphere virtual temperature. Temperature of the lower troposphere better represents the baseline climate since the main differences in SEB comes from less clear-sky longwave radiation in ARCTIC, which primarily arises from troposphere temperatures within a few kilometers above the surface (Ohmura, 2001). In comparison, air temperature at 2 m is largely impacted by surface temperature at or close to the melting point, thus might not be a good representation of baseline climate. Nevertheless, the albedo feedback parameter is defined using the 2 m air temperature, which can be visualized by a straight line connecting the starting and ending points in Figure 14b. The slopes of the straight lines in Figure 14b are steeper for the $$1{}^{\circ}$$ runs (F09 $${=}-0.020$$ $${\mathrm{K}}^{-1}$$; F09M $${=}-0.022$$ $${\mathrm{K}}^{-1}$$; ARCTIC $${=}-0.017$$ $${\mathrm{K}}^{-1}$$) consistent with their larger albedo feedback compared to the ARCTIC run.

[IMAGE OMITTED. SEE PDF]

The impact of the cooler baseline in ARCTIC on the albedo feedback can be understood by adding the initial differences in temperature from the $$1{}^{\circ}$$ runs to the starting and ending temperatures in the ARCTIC run, and using the polynomial fit of albedo to temperature to find the new albedos from these new temperatures. These alternative ARCTIC scenarios are shown as dotted lines in Figure 14b. The slopes of these alternative scenarios (ARCTIC+dT2m $${=}-0.018$$ $${\mathrm{K}}^{-1}$$ in both cases) are more similar to the ARCTIC run than the $$1{}^{\circ}$$ runs, indicating that the cooler baseline in ARCTIC is not responsible for its reduced albedo feedback.

The difference in albedo feedback between ARCTIC and the $$1{}^{\circ}$$ runs can largely be attributed to how surface topography is represented at varying grid resolutions. Figures 15a and 15b show the cumulative area-surface elevation relationships of the GrIS in CISM and CAM, respectively. A steeper slope in the cumulative area-elevation relationship indicates less area increase per meter of elevation rise, reflecting steeper topography. The similar slopes in Figure 15a indicate that the ice sheet topography is similar in CISM among the simulations. When the CAM ice sheet topography is updated based on the CISM topography using the CESM topography software (Lauritzen et al., 2015), a smoothing and flattening process of the raw topography is applied, which is required to prevent grid-scale numerical instabilities in the atmospheric dynamical core (Herrington et al., 2022; Lauritzen et al., 2015). This results in a flatter CAM ice sheet topography relative to the original ice sheet topography in CISM (Figures 15a and 15b), with elevated ice sheet margins and a lower interior (Figure S8 in Supporting Information S1), and also extends the ice sheet margins in CAM beyond the true ice sheet margins (Figure S9 in Supporting Information S1). Since the smoothing and flattening is stronger for lower-resolution grids, F09M and F09 always have flatter CAM ice sheet topography than ARCTIC within the elevation range of ELA variation (Figure 15b), indicating a greater increase in area for the same elevation rise. The slope differences between ARCTIC and the $$1{}^{\circ}$$ runs become even larger at the end of the simulation.

[IMAGE OMITTED. SEE PDF]

The topography smoothing and flattening of the GrIS in CAM is nontrivial, especially in lower-resolution grids. To quantify the ice sheet surface elevation differences between CAM and CISM, the root mean square error (RMSE) between CAM and CISM (remapped to CAM grid for calculation) surface elevation fields over the GrIS was calculated. As shown in Figure 15c, CAM-CISM elevation differences are much smaller in ARCTIC, with an RMSE less than half of those in the $$1{}^{\circ}$$ runs, indicating that ARCTIC better captures the ice sheet topography in CAM and requires fewer temperature and moisture corrections during EC downscaling. To assess how these discrepancies compare with elevation changes due to warming, we also computed the RMSE between annual CISM elevation fields and their initial conditions (light-colored lines in Figure 15c). In the $$1{}^{\circ}$$ runs, warming-induced elevation reduction reaches a similar magnitude to the CAM-CISM RMSE only by the end of the simulation, whereas in ARCTIC, this occurs a century earlier. These findings highlight the significant impact of CAM's topography smoothing and flattening on GrIS representation, particularly in lower-resolution grids.

EC downscaling in CESM2 with the default fixed lapse rates could cause biases in albedo correction, especially for large CAM-CISM topography discrepancies. Compared to RACMO2.3, the default setting of EC downscaling results in a melt energy gradient (negative) smaller in magnitude due to biases in radiation downscaling (Sellevold et al., 2019). Sellevold et al. (2019) suggested the radiation downscaling biases is because CESM2 underestimates the albedo gradient, since it fails to capture very low bare ice albedo due to a fixed bare ice albedo value (0.5). In our study, all the simulations use the same bare ice albedo. It will not further decrease when downscaled to lower elevations along the margins (where CAM has positive elevation biases compared to CISM). However, over the interior, where CAM has negative elevation biases, elevation correction during downscaling underestimates albedo increase. The larger elevation correction in the $$1{}^{\circ}$$ runs results in greater underestimation of snow albedo compared to ARCTIC (Figure 14a; Figures S2i and S2j in Supporting Information S1), leading to stronger albedo feedback and larger melt. It indicates that the current EC downscaling with fixed lapse rates is not sufficient to correct the biases resulting from topography smoothing and flattening in CAM, especially for lower-resolution grids.

The higher snow albedo in ARCTIC results in slower expansion of the ablation zone. Ryan et al. (2019) demonstrates the dominant role of Greenland's seasonally fluctuating snowline in reducing ice sheet albedo compared to bare ice albedo reduction caused by melt processes. Here, instead of the end-of-summer snowline elevation and bare ice area, we use the ELA and ablation area to avoid potential challenges in snowline classification. Although they differ due to processes such as superimposed ice formation (Cogley et al., 2011), a significant correlation exists between the ELA and end-of-summer snowline elevation (Fausto & the PROMICE team*, 2018). Over the 350 years, the average CISM ELA rises from $${\sim} $$750 m to over 2,000 m, with F09M having the largest ELA increase and ARCTIC having the smallest (Figure 16a). The smoothed and flattened topography and the baseline climate affect the average CAM ELA. ARCTIC exhibits slower ablation zone expansion in both CISM and CAM (Figure 16b). Figure 16c shows the relationships between annual ablation area and mean ELA. In CAM, the steeper regression slope indicates a smaller ablation area increase with the same amount of ELA increase. Since CISM receives the interpolated SMB calculated by CLM, a steeper regression slope is also detected using CISM SMB in ARCTIC. This results in more melt in the $$1{}^{\circ}$$ runs, making their ice sheet topography in CISM slightly steeper than ARCTIC (Figure 15a), opposite to CAM topography. We also tested that calculating the annual mean ablation area and ELA in CAM using ARCTIC outputs remapped to the f09 grid (not shown) via the Earth System Modeling Framework first-order conservative remapping algorithm (ESMF Joint Specification Team, 2021) results in a similar steeper regression slope in ARCTIC, confirming that the result is independent of the grid used to calculate ELA.

[IMAGE OMITTED. SEE PDF]

In addition to topography, other factors can also affect albedo feedback. Along the southeast coast, the larger precipitation (as snowfall) in F09 (Figure 9b) may more effectively slow albedo reduction, resulting in weaker bulk albedo feedback compared to ARCTIC (Figure 13c). Moreover, the impact of clouds on net shortwave radiation cannot be eliminated from the bulk albedo feedback calculation. However, since the cloud pattern differences that increase downward shortwave radiation in ARCTIC (Figures 10c and 10d) would lead to an opposite pattern from what is observed, they are not considered a contributor to the albedo feedback differences.

Given the regional dependence of ELA and surface topography, we also analyzed the cumulative area-elevation relationship in CAM for individual GrIS drainage basins (Figure S10 in Supporting Information S1) as defined by Rignot and Mouginot (2012). The slope difference between ARCTIC and the $$1{}^{\circ}$$ runs is greatest in the steepest Central East and South East Basins (Figures S10c and S10d in Supporting Information S1), but due to their narrow ablation zones and large precipitation, it does not result in substantial albedo feedback differences. In basins like the South West, where albedo feedback is weaker in ARCTIC, steeper cumulative area-elevation slopes are consistently observed compared to the $$1{}^{\circ}$$ runs (Figure S10e in Supporting Information S1).

Summary and Discussion

In this study, we examine the impact of enhanced atmospheric resolution on the sensitivity of the GrIS to a changing climate, using a coupled Earth system/ice sheet model (CESM2.2-CISM2.1) with refined resolution over the Arctic. The variable-resolution grid has a horizontal resolution of 1/$$4{}^{\circ}$$ over the broader Arctic region and $$1{}^{\circ}$$ elsewhere. The simulation was conducted under a multicentury idealized warming scenario (1% $${\text{CO}}_{2}$$ increase to 4$${\times} {\text{CO}}_{2}$$) and compared with two reference simulations using the $$1{}^{\circ}$$ lat-lon grid. One of these reference simulations is from the baseline study by Muntjewerf, Sellevold, et al. (2020), enabling an exploration of the impact of enhanced horizontal resolution.

Despite differences in the magnitudes and timing of responses, the behavior of the GrIS under the warming climate in the variable-resolution run is similar to the findings of Muntjewerf, Sellevold, et al. (2020). A “break point” in the mass loss rate trends occurs around year 110, driven by the activation of a strong melt/albedo feedback and enhanced turbulent heat fluxes that accelerates surface melt. Our study confirms the findings of the 1° model with regard to the mechanisms of melt increase in response to warming.

The enhanced resolution generally reduces the mass loss of the GrIS. In the variable-resolution run, the cumulative contribution from the GrIS to global mean sea level rise is 53 mm by year 150 and 831 mm by year 350, about 40% and 20% smaller than the $$1{}^{\circ}$$ runs, respectively. The sea level rise contribution from our variable-resolution run by year 150 is also smaller compared to other CESM simulations by 2100 under the CMIP RCP8.5 and SSP5-8.5 scenarios, when $${\text{CO}}_{2}$$ concentration approaches 4$${\times} {\text{CO}}_{2}$$. In comparison, Lipscomb et al. (2013) reports 76 mm and Muntjewerf, Petrini, et al. (2020) reports 109 mm sea level rise. When compared to projections using stand-alone GrIS models forced by CMIP5 GCM outputs, our estimate approaches the lower bound (90$$\pm $$50 mm during 2015–2100; Goelzer et al., 2020). The annual MB differences between the variable-resolution run and the $$1{}^{\circ}$$ runs (Table 1) exceed the standard deviation of annual MB from the last 100 years of the pre-industrial simulations ($$\pm $$0.30 mm $${\text{yr}}^{-1}$$ for ARCTIC, $$\pm $$0.27 mm $${\text{yr}}^{-1}$$ for F09), and the final fully coupled segment of the long spun-up pre-industrial Earth system/ice sheet state ($$\pm $$0.23 mm $${\text{yr}}^{-1}$$; Lofverstrom et al., 2020). Therefore, this difference can be attributed primarily to differences in grid resolution rather than model internal variability.

The reduced mass loss in the variable-resolution run primarily results from reduced surface melt during summer, concentrating in the western and northern transitional zones. Compared to the variable-resolution run, the larger surface melt in the $$1{}^{\circ}$$ runs is primarily driven by greater underestimation of snow albedo, characterized by stronger albedo feedback and faster expansion of the ablation zone. The smoothing and flattening of the ice sheet topography in CAM in the $$1{}^{\circ}$$ runs leads to larger CAM-CISM topography discrepancies over the GrIS. Since EC downscaling with the default fixed lapse rates exhibits a biased low albedo gradient compared to RACMO2.3 (Sellevold et al., 2019), larger topography correction in the $$1{}^{\circ}$$ runs during EC downscaling leads to greater underestimation of snow albedo. In addition to albedo, the steeper topography in ARCTIC also moves cloud further offshore. If without the differences in albedo, net shortwave radiation would have increased in ARCTIC due to enhanced incident shortwave radiation from cloud changes. Therefore, considering the crucial impacts of topography, future sea level projections based on coupled models with coarse resolution may be overestimated due to their inability to resolve the topography of Greenland adequately. Stand-alone ice sheet model simulations may benefit from using SMB and temperature forcing from regionally refined climate models, which better resolve steep topographic gradients, thus bypassing the use of RCMs for downscaling. Since the current EC downscaling in CESM2 with fixed lapse rates is insufficient at correcting the biases resulting from topography smoothing and flattening in CAM, EC downscaling with varied lapse rates or different downscaling methods are worth investigating.

Comparisons between these simulations are complicated by differences in grid resolution, physics time step, and dynamical core. Like increasing resolution, reducing the physics time step can also increase resolved vertical velocities and thus condensational heating. By comparing two AMIP-style CESM2.2 simulations using the same quasi-uniform $$1{}^{\circ}$$ SE grid but different time steps, Herrington et al. (2022) showed that the simulation with a reduced time step produced a warmer troposphere over nearly all latitudes. In a comparison of two simulations using the same time step but different grids - the ARCTIC grid and the quasi-uniform $$1{}^{\circ}$$ grid, Herrington et al. (2022) revealed that the temperature increase caused by enhanced resolution was confined to the refined Arctic region and had a larger magnitude. Therefore, the differences in the simulated climate within the refined Arctic region in this study are more likely a result of changing the horizontal resolution, although the impact of changing the physics time step cannot be ruled out. Over Antarctica, simulations using the FV and SE dynamical cores on $$1{}^{\circ}$$ grids show no significant differences in major cloud properties and circulation (Datta et al., 2023). Herrington et al. (2022) compared a set of AMIP style simulations with six different grids from the FV and SE dynamical cores and found larger historical melt and precipitation biases (excess melt and precipitation) in the $$1{}^{\circ}$$–$$2{}^{\circ}$$ runs using the SE dynamical core than those using the FV dynamical core, and these melt biases are greatly reduced with regional refinement (e.g., ARCTIC grid). This adds confidence to our attribution of the smaller melt in the variable-resolution run to finer resolution rather than differences in the dynamical core. However, it is impossible to fully disentangle the impacts of changing resolution and dynamical core in the current simulations, as these factors vary simultaneously.

Another uncertainty in this study is the initial conditions. The pre-industrial simulation of ARCTIC was branched off from an initial state similar to that of the F09M 1% $${\text{CO}}_{2}$$ simulation. This aims to achieve a near-equilibrium state of the GrIS after the grid change. As a more direct comparison, F09 underwent a similar spin-up process as ARCTIC. The resulting initial ice sheet conditions before the idealized warming period differ across the three simulations: F09 has a larger initial ice volume than ARCTIC, while F09M has a smaller initial ice volume than ARCTIC (Figure 8d). We find that ARCTIC exhibits a slower SMB decrease and lower melt compared to the $$1{}^{\circ}$$ runs over multicentury timescales, regardless of whether their initial ice volume is larger or smaller. Additionally, ARCTIC has a cooler pre-industrial climate than F09M and F09. The impact of this cooler climate persists throughout the idealized warming simulation, as shown in Figures 6g–6l. Moreover, we note that all three simulations exhibit small positive drifts in GrIS MB before the warming scenario begins (Table 1). The cooler initial climate in ARCTIC leads to a longer adjusting period for the GrIS to transition to mass loss, which degrades the precision of specific estimates for GrIS sea level rise contribution. However, our results suggest that atmospheric temperature differences are not the main driver of the different GrIS responses among the simulations. Taking these factors into account, we conclude that our findings are robust. A possible future application is a full spun-up simulation using the ARCTIC grid, which has the potential to provide a more realistic pre-industrial GrIS topography.

One limitation of the current model configuration lies in the ice-ocean interface. The direct impact of oceanic thermal forcing on ocean-terminating ice fronts is not included, and the floating criterion used in the calving parameterization is highly idealized. The limited understanding and implementation of processes such as calving and submarine melting in ice sheet models have been identified as a major source of uncertainty for future projections of the GrIS (Goelzer et al., 2020). Oceanic forcing can enhance solid ID (Holland et al., 2008; Wood et al., 2018). Using another coupled Earth system/ice sheet model, EC-Earth-PISM, under the same 1% $${\text{yr}}^{-1}$$ $${\text{CO}}_{2}$$ warming scenario, Madsen et al. (2022) showed that even by embedding a constant oceanic thermal forcing and a simple geometric calving criterion, ID would first increase and then decrease, resulting in a much smaller reduction after 350 years. Most modeling studies that include oceanic forcing estimate a secondary future sea level rise contribution from ice dynamics compared to SMB for the entire ice sheet (Aschwanden et al., 2019; Fürst et al., 2015; Price et al., 2011). However, with improved bathymetry and bed topography mapping, Choi et al. (2021) showed that ice dynamics could contribute as much as, or even more than, SMB to GrIS mass loss over this century. These studies illustrate the importance of properly including ice-ocean interactions for future GrIS projections on century to multicentury timescales. Currently, the functionality for ice sheet-ocean interactions is being developed in CESM and will be included in future model versions. The ARCTIC grid, or grids with even higher-resolution refinement, will also be beneficial in providing more accurately resolved atmospheric forcing for modeling ice dynamics in narrow fjords.

Overall, our study demonstrates the value of employing variable-resolution grids for coupled climate-GrIS modeling, offering critical insights into ice sheet-climate interactions. It highlights the importance of grid resolution in modeling the evolution of the GrIS on multicentury timescales, particularly in capturing topography-related processes and feedbacks, and thus advances projections of the GrIS's future sea level rise contribution.

Acknowledgments

AH and AG contributions to this material are supported by the National Center for Atmospheric Research (NCAR), a major facility sponsored by the NSF under Cooperative Agreement no. 1852977. Computing and data storage resources, including the Cheyenne supercomputer (Computational and Information Systems Laboratory, 2017), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR, as part of the System for Integrated Modeling of the Atmosphere project. ZY, ACS, JT, and RT are supported by the NSF (Grant OPP, 1952199) and HDR iHARP institute (Grant OAC, 2118285). We thank the editor Stephen Griffies, associate editor Nick Golledge, reviewer Miren Vizcaino, and the other two anonymous reviewers for their very useful comments that helped improve the manuscript.

Data Availability Statement

CESM2 is an open-source model, available via the CESM GitHub repository (CESM Project, 2024). The F09 experiment was done with tag “cesm2.1.3”; the ARCTIC experiment was done with tag “cesm2_2_alpha06d.” The data presented in this manuscript are stored in an open repository Zenodo (Herrington & Yin, 2024), and the code for generating the plots are available at Yin (2024).

References

Aschwanden, A., Fahnestock, M. A., & Truffer, M. (2016). Complex Greenland outlet glacier flow captured. Nature Communications, 7(1), [eLocator: 10524]. [DOI: https://dx.doi.org/10.1038/ncomms10524]