## 1. Introduction

With about two-thirds of global mean sea level rise being due to ice mass loss in the period 1900–2018 and the rest being mainly attributed to thermal expansion (Frederikse and others, Reference Frederikse2020; Fox-Kemper and others, Reference Fox-Kemper2021), the knowledge of future glacier and ice-sheet mass loss is of utmost importance. The Greenland Ice Sheet holds an ice volume of 7.4 m of global mean sea level equivalent (SLE) (Morlighem and others, Reference Morlighem2017) and has contributed 10.8 ± 0.9 mm to the mean sea level rise between 1992 and 2018, with 49.7% of the mass loss caused by dynamic discharge into the ocean and the rest by melt and runoff (The IMBIE Team, 2020; Fox-Kemper and others, Reference Fox-Kemper2021). Polar amplification has caused the Arctic to warm nearly four times as fast as the rest of the Earth (Rantanen and others, Reference Rantanen2022) and with the global temperatures expected to rise further due to increased anthropogenic forcing, the sea level contribution from the Greenland Ice Sheet will continue to increase in the coming centuries. Global warming also induces increased variability of climate, including more extreme events of both warm and cold temperatures, which are expected to be more frequent in the future (Seneviratne and others, Reference Seneviratne2021). Increased climate variability is expected to cause the ice sheet to melt even faster than today (Mikkelsen and others, Reference Mikkelsen, Grinsted and Ditlevsen2018; Beckmann and Winkelmann, Reference Beckmann and Winkelmann2022). Atmospheric warming causes the surface melting to increase when the temperature is above the freezing point, while a similar cooling would not affect the melting for temperatures below the freezing point. Consequently, a surface temperature varying around the freezing point will tend to cause the melting to increase and lead to net mass loss. However, the response of the ice-sheet mass balance to increased temperature variability is complex and highly non-linear, and if the snow accumulation rate increases in warmer climates, higher variability could also result in a larger ice sheet (Albrecht and others, Reference Albrecht, Winkelmann and Levermann2020). An earlier study investigated this asymmetric response to temperature fluctuations using a simplified, perfectly plastic ice-sheet model and concluded that the Greenland Ice Sheet would become 0.5–1 m SLE smaller when including the inter-annual temperature variability than without this variability (Mikkelsen and others, Reference Mikkelsen, Grinsted and Ditlevsen2018). It is unclear what the influence of including temperature variability would be in a model of the Greenland Ice Sheet that includes ice dynamics and discharge. Most models in the Ice Sheet Modelling Intercomparison Project for the IPCC's sixth Assessment Report (ISMIP6) had a mass loss lower than the observed in the historical period below (Aschwanden and others, Reference Aschwanden, Bartholomaus, Brinkerhoff and Truffer2021), and in some cases, this was traced back to a lack of variability in the climate model forcing over the historical period (Goelzer and others, Reference Goelzer2020). Another recent study found that the Greenland Ice Sheet volume decreased when extremely warm summers were included in the temperature forcing (Beckmann and Winkelmann, Reference Beckmann and Winkelmann2022), but the effect of cold spells was not included.

In this paper, we examine the effect of inter-annual variability in atmospheric temperatures on the present-day Greenland Ice Sheet geometry and volume by using a thermodynamically coupled three-dimensional ice-sheet model to account for ice-dynamical feedbacks. Our goal is to revise the previous results for a simplified ice-sheet model (Mikkelsen and others, Reference Mikkelsen, Grinsted and Ditlevsen2018) by using a dynamic ice flow model with realistic geometry and mass-balance forcing. This is needed to assess the importance of including climate variability in ice-sheet simulations over the historical period and can inform efforts to develop a new protocol for projections of ice-sheet mass loss in ISMIP7 (Aschwanden and others, Reference Aschwanden, Bartholomaus, Brinkerhoff and Truffer2021).

We construct an ensemble of temperature-forcing time series that mimics the mean annual temperature variability found in the NOAA-CIRES 20th Century Reanalysis (20CR), version 2c (Compo and others, Reference Compo2011) following the method of Mikkelsen and others (Reference Mikkelsen, Grinsted and Ditlevsen2018). We add this variability to the seasonal mean temperatures of the period 1960–1989 from the RACMO model (Noël and others, Reference Noël, van de Berg, Lhermitte and van den Broeke2019). We then perform an ensemble of ice-sheet simulations, where each ensemble member is forced with a corresponding temperature-forcing time series ensemble member. We compare this simulated ice-sheet ensemble to a control run forced with constant climate without inter-annual temperature variability to assess the impact on the ice-sheet volume, and we identify the most sensitive ice-sheet areas. The spatial differences in the response to the inter-annual temperature forcing are then analyzed in six regional sectors of the Greenland Ice Sheet. Finally, we investigate how temperature variability influences the response of the ice sheet to an abrupt warning.

## 2. Methods

### 2.1. The ice-sheet flow model

To determine the influence of inter-annual temperature variability on the modeled present-day Greenland Ice Sheet we use the three-dimensional Parallel Ice Sheet Model (PISM) (www.pism.io). PISM is an open-source thermodynamically coupled ice-sheet model that applies two shallow approximations to the Stokes equation, both assuming a small ice-sheet aspect ratio. The Shallow Ice Approximation (SIA) assumes that the horizontal shear stresses and the longitudinal deviatoric stresses are negligible compared to the vertical shear stresses. This assumption holds well for most of the interior of the Greenland Ice Sheet where there is limited basal sliding. The Shallow Shelf Approximation (SSA), on the other hand, assumes that the shear stresses are negligible compared to the longitudinal or membrane stresses of the ice flow (Schoof, Reference Schoof2006). In the hybrid model incorporated in PISM, the horizontal velocity of the ice is given by a weighted average of the SIA velocity and the SSA velocity using the SSA as a sliding law (Bueler and Brown, Reference Bueler and Brown2009). The constitutive relation for the rheology of ice used in PISM is based on Glen's flow law (Glen, Reference Glen1955), modified by Nye (Nye, Reference Nye1957) and Lliboutry and Duval (Lliboutry and Duval, Reference Lliboutry and Duval1985), or formally the Glen–Paterson–Budd–Lliboutry–Duval law:

where *n* is the flow law exponent, which we take to be 3 for both shallow approximations, *τ* _{e} is the effective deviatoric stress and *τ* _{ij} are the components of the deviatoric stress tensor. *A* is the ice softness which is given piecewise for cold and temperate ice by two Arrhenius functions of the pressure, *P*, temperature, *T* and liquid water fraction, *w*, (Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012). The enhancement factor, *E*, is set to 3 for the SIA and 1 for the SSA. The study by Aschwanden and others (Reference Aschwanden, Fahnestock and Truffer2016) recommended an SIA enhancement factor of *E* = 1.25 for the Greenland Ice Sheet when tuned to match the ice flow velocity of the outlet glaciers, but we use their recommended set of parameters for *E* = 3, which overall provides the best fit to the Greenland Ice Sheet ice volume, see discussion. The basal sliding velocity **u**_{b} from the SSA is related to the basal shear stress **τ**_{b} and the yield stress through a power friction law, where we use a pseudo-plasticity exponent *q* of 0.6 following Aschwanden and others (Reference Aschwanden, Fahnestock and Truffer2016). The yield stress is given by the Mohr–Columb model and is modified by a till friction angle to make the till softer in lower-lying regions. We parameterize the till friction angle as a continuous function of the bed topography that increases linearly from $5^\circ$ to $40^\circ$ between 700 m below and 700 m above sea level following Aschwanden and others (Reference Aschwanden, Fahnestock and Truffer2016). For bed topography, we use IceBridge BedMachine Greenland, version 4 (Bamber and others, Reference Bamber2013), while the geothermal heat flux distribution is taken from Shapiro and Ritzwoller (Reference Shapiro and Ritzwoller2004). Our model setup does not account for calving through a representation of the physical processes, instead, we remove all ice that exceeds the present-day boundary using the prescribed front retreat option in PISM, similar to the ‘retreat implementation’ used in ISMIP6 (Nowicki and others, Reference Nowicki2020).

### 2.2. Climate forcing

For the surface boundary conditions, we use the 2 m air temperature field and precipitation field from RACMO2.3p2 at 5.5 km spatial resolution (Noël and others, Reference Noël, van de Berg, Lhermitte and van den Broeke2019). We use a 12-month climatology based on the multi-year monthly averages of temperature and precipitation for the period 1960–1989, in which the Greenland Ice Sheet volume is considered to be close to equilibrium (The IMBIE Team, 2020). This is used as the climate forcing in the model initialization. The temperature and precipitation fields are used to calculate the surface mass balance (SMB). All precipitation is taken to be snow accumulation for temperatures below $0^\circ$C, while the fraction of snow is taken to linearly decrease to zero between $0^\circ$C and $2^\circ$C. Surface melt is calculated with a positive degree-day scheme which assumes that the amount of melt is proportional to how much the temperature exceeds the freezing temperature. The positive degree day factors used are 3.3 mm K^{−1} d^{−1} for snow and 8.8 mm K^{−1} d^{−1} for ice. In addition to the inter-annual variability and the seasonal cycle, the temperatures also fluctuate daily, so while the monthly mean temperature might be below freezing, some days it will be above, resulting in melting. Therefore, it is assumed that the daily temperature is given by a normal distribution around $\bar {T}$ with variance over the month of $\sigma _{\rm pdd}^2$. The number of positive degree days, *D*, for a time interval Δ*t* = *t* _{2} − *t* _{1} is given by

We use a constant uniform positive degree day standard deviation of *σ* _{pdd} = 5 K (Payne and others, Reference Payne2000). To adjust the temperature for changes in surface altitude throughout the simulation, we use a lapse rate of 6.5 K km.

### 2.3. Including inter-annual temperature variability in the forcing field

To assess the effect of inter-annual temperature variability on the Greenland Ice Sheet volume evolution, an ensemble of 50 temperature time series with a length of 10 ka, denoted surrogates, are constructed based on the statistics of the 20CR over the period 1851–2014 (Compo and others, Reference Compo2011). These are then added to the 12-month climatology from RACMO with mean monthly temperatures averaged over 1960–1989, which does not include inter-annual variability. The resulting climate-forcing fields represent the 1851–2014 inter-annual temperature variability over Greenland, and we use these climate-forcing fields to investigate the effect of realistic climate variability in a model of the Greenland Ice Sheet. By using an ensemble of idealized, detrended forcing fields, we can assess the modeled ice-sheet response and take into account the stochastic nature of the climate variability.

The inter-annual temperature variations of the temperature surrogates, *T*, are modeled as a function of time, *t*, according to the first-order auto-regressive model (Hasselmann, Reference Hasselmann1976; Frankignoul and Hasselmann, Reference Frankignoul and Hasselmann1977)

where *ϕ* is the auto-regressive parameter and $\varepsilon _t$ is white noise drawn from a Gaussian distribution with zero mean and variance $\sigma ^2_\varepsilon$. We use the statistics of the full period of the 20CR to generate the 10 ka time series. The variability found in the 20CR dataset is non-uniform over Greenland as shown in Figure 1 (a), so we follow the procedure by Mikkelsen and others (Reference Mikkelsen, Grinsted and Ditlevsen2018) and average over a box spanning from 68 to 80°N and from 25 to 60°W. We then de-trend this time series by subtracting its least squares linear fit and set *c* = 0, such that the temperature surrogates have zero mean. From the resulting time series, we estimate the parameters in (3) to be *ϕ* = 0.67 and $\sigma _\epsilon = 0.92$ K. The 50 temperature surrogates constructed using these parameters all have a standard deviation of 1.24 K. The de-trended time series and a subset of the constructed temperature surrogates used to force the positive degree day model are shown in Figure 1 (b). To account for the variability being non-uniform, we scale the time series to match the temporal standard deviation in each grid point, i.e. we compute the outer product of the standard deviation map in Figure 1 (a) and the generated time series divided by their standard deviation. In addition to the 50 surrogates based on the 20CR data, we also construct an ensemble of 50 surrogates with twice the standard deviation representing another background climate state with more pronounced extreme events denoted 20CRx2.

### 2.4. Model initialization

The ice-sheet model is initialized prior to the experiments with the inter-annual temperature time series. As a first step in the model initialization, the ice-sheet model is run for 50 ka forced with the 12-month climatology at 9.6 km resolution followed by a 20 ka run at 4.8 km resolution to obtain a state–state Greenland Ice Sheet close to equilibrium. The steady-state volume fluctuates around equilibrium with a standard deviation of 2.95 mm SLE for the 9.6 km resolution and 0.33 mm mm SLE for the 4.8 km resolution. As mentioned above, we restrict the ice sheet to the domain of the present-day ice sheet and remove all ice beyond the present-day margin. This boundary condition leads to a larger amount of ice discharge compared to observations but is required to keep the geometry of the ice sheet close to the observations. The resulting steady-state Greenland Ice Sheet has decreased in volume to 7.07 m SLE relative to the observed present-day Greenland Ice Sheet. It is thinner in the interior and thicker close to the margin compared to the observed Greenland Ice Sheet as shown in Figure 2 with an overall RMSE of the elevation difference between modeled and observed thickness of 247 m.

The overall shape of the ice sheet depends on the SIA enhancement factor, *E*. As mentioned above, we use *E* = 3 in our study, which overall provides the best fit to the Greenland Ice Sheet, but we also made an additional set of runs using *E* = 1.25 to test the effect of the enhancement on our results. As shown in Figure 3, decreasing the softness to *E* = 1.25 in the steady-state run results in an ice sheet that is thicker than both the observed ice sheet and the one modeled with *E* = 3. Consequently, the total SMB over the Greenland Ice Sheet is 21.2 Gt a^{−1} larger with *E* = 1.25 compared to *E* = 3 due to the surface elevation feedback; this results in an equally larger discharge at the front, since we are investigating steady-states.

### 2.5. Experiments

After the model initialization, the simulation is continued with an ensemble of runs where the temperatures are uniformly offset from the 12-month climatology by the two ensembles of temperature surrogates. Like the ensembles of temperature surrogates, we denote the two ensembles of Greenland Ice Sheet simulations, 20CR and 20CRx2. We run the ensembles for 10 kyrs with a resolution of 4.8 km until the ensemble means have reached new steady states. In addition, an unperturbed control run is also considered to assess the effect of the inter-annual temperature variability.

The second experiment starts from the same initialized ice sheet as the above experiment but with an instantaneous and uniform change in temperature forcing of [0, 0.5, 1, 1.5, 2] K in order to assess the response of the Greenland Ice Sheet to changes in temperature with and without inter-annual temperature variability. These simulations are all carried out over 500 years after the instantaneous temperature change.

## 3. Results

### 3.1. Steady state

The evolution of ice volume of the two ensembles is shown in Figure 4. After 10 ka, the ice volumes are 1.9 ± 0.4 and 11.5 ± 1.4 cm SLE smaller relative to the control for the 20CR and 20CRx2 simulations, respectively, where the first number denotes the ensemble mean and the second number the ensemble standard deviation. The volume changes correspond to mean elevation changes of 1.8 ± 0.3 and 10.8 ± 1.3 m for the 20CR and 20CRx2, respectively. This shows that the modeled ice-sheet response to variability is strongly non-linear. A doubling of the variability results in an approximately six times larger response.

As shown in Figure 5, the temperature-variability-induced mass losses from the SMB and basal mass balance are fully compensated by reduced discharge at the end of the simulation, resulting in the steady states. This result shows that the increased melt in ensemble runs is being compensated for by reduced discharge, i.e. a negative feedback caused by calving dynamics.

The thickness difference between the control run and the ensemble mean is shown in Figures 6(a,b) and the ensemble standard deviations in Figures 6 (c,d). The ensemble mean simulated ice sheet is thinner than the control run and the largest difference is in the north region, but for the 20CRx2 the difference is also considerable on the west coast. From this, we see a clear thinning induced by the temperature variability, especially in the North for both ensembles. We divide the Greenland Ice Sheet into the six basins (Mouginot and Rignot, Reference Mouginot and Rignot2019) shown in Figure 6, with their respective mean ice thickness changes listed in Table 1. The Northern (NO) basin holds 11.6% of the initialized ice volume but contributes to about 40% of the mass loss during the simulations with inter-annual variability in the temperature forcing. In contrast, the Southeast (SE) basin holds 12.8% of the initialized ice volume but contributes only about 1% of the total mass loss during the simulations. Figure 7 shows the SMB for the control run and the ensemble mean for the first 10 years. The distribution of mass loss negatively correlates with the SMB, i.e. areas with high mass loss correspond to areas with a low SMB.

The effect of including climate variability on the Greenland Ice Sheet volume is similar for the *E* = 1.25 runs compared to the *E* = 3 runs; the mean ensemble ice volume is 2.0 ± 0.4 cm smaller for the 20CR run and 8.8 ± 0.9 cm smaller for the 20CRx2 run compared to the control simulation with *E* = 1.25.

We also tested a uniform variability over Greenland, following Mikkelsen and others (Reference Mikkelsen, Grinsted and Ditlevsen2018), which resulted in a slightly larger effect.

### 3.2. Response to instantaneous changes in temperature forcing

The transient behavior of the simulations when forced with an instantaneous change in the temperature forcing, with and without the inter-annual variability is shown in Figure 8 where the average rate of change in mass for the first 100, 300 and 500 years are shown. The mass loss rate increases with increasing temperature step. For the first 100 years, the average rate of mass loss per degree warming is ~0.15 mm SLE a^{−1} K^{−1} for the control run, and 0.014–0.17 mm SLE a^{−1} higher than the control run when adding the inter-annual variability in the temperature forcing depending on the size of the temperature step. For the 20CRx2, the rate of change in mass is almost twice as large as the rate of the control run. For each experiment with an instantaneous change in temperature, the control run lies within the range of the ensemble runs for the first 100 years, close to the lower bound of mass loss rate of the 20CR runs for the 300 and 500 years, and below the mass loss rates of the 20CRx2 runs for the 300 and 500 years.

## 4. Discussion

In this study, we investigate the effect of climate variability on the steady-state volume of the Greenland Ice Sheet. Most models in the ISMIP6 study had a mass loss lower than observed over the historical period (Goelzer and others, Reference Goelzer2020). In some cases, this difference could be explained by the fact that climate forcing over the historical period did not include climate variability on interannual and decadal timescales, but the effect from the historical mass loss being inaccurate was effectively removed by considering only the difference from a control run (Goelzer and others, Reference Goelzer2020). In order to improve projections of the ice-sheet mass loss, one step could be to include the observed variability in the forcings over the historical period in order to improve the model initialization. Our results provide a quantified assessment of this effect.

The main goal of our study was to revise the previous results for a simplified Greenland Ice Sheet (Mikkelsen and others, Reference Mikkelsen, Grinsted and Ditlevsen2018) by considering a more appropriate state-of-the-art ice flow model for the Greenland Ice Sheet. We found that the total volume of the Greenland Ice Sheet is reduced by 1.9 ± 0.4 cm SLE after 10 ka when including the inter-annual variability of the 20CR in temperature forcing. This is much less than the result of 0.5–1 m in Mikkelsen and others (Reference Mikkelsen, Grinsted and Ditlevsen2018), showing that a more complex model in our study, including ice-flow dynamical feedback, reduces the effect of temperature variability on the ice-sheet volume. Our study differs from Mikkelsen and others (Reference Mikkelsen, Grinsted and Ditlevsen2018) by the complexity of the ice-sheet model. They use a simplified model (Oerlemans, Reference Oerlemans2003), assuming that ice flows as perfectly plastic material and that the ice sheet is axially symmetric with a bedrock that slopes linearly downward. They coupled their model to atmospheric variations only by changing the equilibrium line altitude linearly according to variations in global temperature and then calculating the corresponding change in the steady-state volume. The simplified ice-sheet model used by Mikkelsen and others (Reference Mikkelsen, Grinsted and Ditlevsen2018) does not have ice–ocean interaction or calving included, a process that plays a significant role for the Greenland Ice Sheet, where only half of the mass loss is due to ablation (The IMBIE Team, 2020).

In our model setup, the influence of inter-annual temperature variability on the total mass balance of the Greenland Ice Sheet is reduced by ice discharge, since an increase in the melt due to temperature variability will result in thinning resulting in reduced discharge. The feedback from discharge to the mass loss induced by including inter-annual temperature variability is thus negative. This negative feedback is not accounted for in the study by Mikkelsen and others (Reference Mikkelsen, Grinsted and Ditlevsen2018), which may explain why they find a larger sensitivity to inter-annual temperature variability than we find here. However, the fixed margin calving scheme implemented in our model is not based on a physical calving criterion but only ensures that the Greenland Ice Sheet does calve at the present-day margin.

We speculate that the response to temperature variability is dampened by the negative feedback between SMB change and ice discharge. Mikkelsen and others (Reference Mikkelsen, Grinsted and Ditlevsen2018) lack the important negative feedback from discharge and their result should therefore be considered an upper bound on the sensitivity. Here, we use mask-based calving which results in particularly strong negative feedback as it can accommodate any change in horizontal flux. We, therefore, argue that our results should be taken as a lower bound.

The modeled steady-state ice sheet fluctuates about the equilibrium volume with a standard deviation of 3.05 mm SLE for the 9.6 km resolution run, which reduces to 0.38 mm SLE for the 4.8 km resolution run. Meanwhile, the volume of the ensemble members fluctuates with a mean standard deviation of 3.70 and 10.2 mm SLE for the 20CR and 20CRx2 runs, respectively. It is therefore important to consider simulations at a resolution of at least 4.8 km in order for the mass-loss signal not to be drowned in the fluctuations around the steady state.

The response to the inter-annual temperature variability in ice volume is greatest in the NO basin (Table 1, Fig. 6) where the volume is reduced by 7.8 ± 0.9 and 45 ± 5 mm SLE for the 20CR and 20CRx2 run, respectively, corresponding to a mean elevation reduction of 13 ± 2 and 76 ± 8 m compared to the control run. This is in stark contrast to the SE basin where the volume reduces by 0.2 ± 0.4 and 1.5 ± 1.0 mm SLE for the 20CR and 20CRx2 simulations, respectively, corresponding to a mean elevation reduction of 0.3 ± 0.5 and 2.1 ± 1.4 m compared to the control run. The spatial distribution of mass loss correlates with the spatial distribution of SMB, which is positive for all basins but balanced by an equally large discharge at the ice front.

In Figure 7 the spatial distribution of the SMB is plotted for the control run and the ensemble means, both for July and the annual mean. For July, the ablation zone reaches far inland in the NW, CW, NO and NE basins, while it is more narrowly confined closer to the ice-sheet margin in the SE basin. When including the inter-annual variability in the temperature forcing, the change in SMB is greatest near the equilibrium line altitude where the asymmetry of the positive degree day scheme is most pronounced. The SE basin accounts for 49.6% of the total SMB, while the NO basin only accounts for 2.8% of the total SMB for the control run. Including the inter-annual variability in temperature forcing decreases the mass balance of the ensemble mean by 0.4% for the SE basin and by 6% for the NO basin for the 20CR ensemble relative to the control run. Those numbers increase to 1.3 and 23% for the 20CRx2 ensembles.

The 20CRx2 ensemble simulations are constructed to have twice the inter-annual temperature variability found in the 20CR temperatures for the period 1851–2014 (Fig. 1). In a future warmer climate, the variability of climate may increase, and these simulations are included in order to test the influence of a higher variability than in the recent past. The temperature variability is not uniform over Greenland but is greatest in the Northern areas, which we also found to be the most sensitive. The northern basins have less precipitation than the southern basins, and are thus more sensitive to climate changes that lead to increased surface melting. The SW basin has contributed most to the recent observed mass loss (The IMBIE Team, 2020), but this could be due to a combination of several different factors, not considered in our study that focuses solely on including the temperature variability.

In our study, we did not account for the effect of variation in precipitation. However, we conducted an additional test where we ran an ensemble forced by precipitation anomalies from the RACMO data at 1 km spatial resolution (Noël and others, Reference Noël, van de Berg, Lhermitte and van den Broeke2019) constructed in the same manner as the temperature surrogates. We found the effect to be negligible compared to the effect of inter-annual temperature variability. We also found temporal anomalies in temperature and precipitation to be approximately uncorrelated when averaging over Greenland.

Our study invites further work to investigate the effect of spatial variability in temperature and precipitation forcing fields on the ice sheet and its response to climate changes. Future studies could include the feedback of the atmosphere as the Greenland Ice Sheet evolves (i.e thinning and retreat) in response to temperature and precipitation changes and investigate regional differences in the variability.

## 5. Conclusion

We ran an ensemble of ice-sheet simulations to determine the influence of inter-annual temperature variability on the Greenland Ice Sheet volume. Each ensemble member was forced by a realization of temperature anomalies given by an auto-regressive model based on the inter-annual temperature variability found in the 20CR dataset over the period of 1851–2014 and added to a 12-month climatology from RACMO. Furthermore, we tested the influence of temperature variability on projections by considering a series of ensemble runs where we abruptly increased the temperature by a step change of [0.5, 1.0, 1.5, 2.0] K. We repeated the experiments with twice the inter-annual temperature variability found in the 20CR dataset; these simulations are denoted 20CRx2.

Accounting for inter-annual temperature variability in simulations of the Greenland Ice Sheet leads to a smaller steady-state ice-sheet volume in our study compared to the simulation forced with a constant climatology (all else being equal), as well as higher mass loss rates with a spread among ensemble members when forced with an abrupt change in temperature. We find that the steady-state volume decreases by 1.9 ± 0.4 and 11.5 ± 1.4 cm SLE for the 20CR and 20CRx2 simulations, respectively, when forced with a variable temperature forcing compared to a constant temperature forcing. For the first 100 years, the average rate of mass loss per degree warming is ~0.15 mm SLE a^{−1} K^{−1} for the control run, and 0.014–0.17 mm SLE a^{−1} higher than the control run for the first 100 years when adding the inter-annual variability in the temperature forcing depending on the size of the temperature step.

Our results are based on a complex ice flow model with ice dynamical feedbacks. Compared to previous results with a simpler ice flow model Mikkelsen and others (Reference Mikkelsen, Grinsted and Ditlevsen2018), we find a smaller sensitivity of the simulated ice sheet to variability in the climate forcing. We argue that our result should be interpreted as a lower bound. In a future warmer climate, the climate variability is projected to increase (Seneviratne and others, Reference Seneviratne2021), potentially leading to further increase of the ice-sheet mass loss.

The sensitivity of the Greenland Ice Sheet to inter-annual temperature variability is found to be spatially dependent. The NO basin, experiencing less precipitation and lower SMB than other Greenland Ice Sheet basins, loses 0.9–1.1% of its total volume, corresponding to 40% of the total mass reduction found between the control run and the ensemble mean.

Our results show that temperature variability in the forcing of the ice-sheet model affects both the Greenland Ice Sheet equilibrium volume and its transient response to instantaneous changes in temperature. The present-day Greenland Ice Sheet is the result of past climate changes and variability, and thus projections must take into account the effect on both the initial, present-day ice sheet, as well as in projections of its future evolution. In this regard, including climate variability is important, since it contributes to both the evolution and the uncertainty of future projections of mass loss.

## Acknowledgements

We acknowledge funding from the Independent Research Fund Denmark through the project GreenPlanning (grant no. 0217-00244B). We thank the editor Professor Douglas MacAyeal and an anonymous referee for their comments, which helped improve the manuscript.