North Atlantic Cooling is Slowing Down Mass Loss of Icelandic Glaciers

Icelandic glaciers have been losing mass since the Little Ice Age in the mid‐to‐late 1800s, with higher mass loss rates in the early 21st century, followed by a slowdown since 2011. As of yet, it remains unclear whether this mass loss slowdown will persist in the future. By reconstructing the contemporary (1958–2019) surface mass balance of Icelandic glaciers, we show that the post‐2011 mass loss slowdown coincides with the development of the Blue Blob, an area of regional cooling in the North Atlantic Ocean to the south of Greenland. This regional cooling signal mitigates atmospheric warming in Iceland since 2011, in turn decreasing glacier mass loss through reduced meltwater runoff. In a future high‐end warming scenario, North Atlantic cooling is projected to mitigate mass loss of Icelandic glaciers until the mid‐2050s. High mass loss rates resume thereafter as the regional cooling signal weakens.


10.1029/2021GL095697
2 of 11 eruption of Eyjafjallajökull that covered surrounding ice caps with dark volcanic tephra, reducing surface albedo (Gudmundsson et al., 2012;Möller et al., 2019). In sharp contrast, mass loss has slowed down since 2011 (Aðalgeirsdóttir et al., 2020;Ciraci et al., 2020;Hock et al., 2019;Wouters et al., 2019), but as of yet little is known about whether this slowdown is temporary, if it will persist or further strengthen in the future.
To address this, we use the Regional Atmospheric Climate Model RACMO2.3 (see Section 2.1). For the recent past , the model is forced by climate reanalyses including the latest product of the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 (Hersbach et al., 2020;RACMO2.3-ERA henceforth;Noël et al., 2015). For projections throughout the 21st century , RACMO2.3 is forced by the output of the Community Earth System Model CESM2 (see Section 2.2) under a high-end Shared Socioeconomic Pathway SSP5-8.5 emission scenario (RACMO2.3-CESM2 henceforth; Noël et al., 2021). Using a high-end emission scenario enables us to explore the response of Icelandic glaciers to unmitigated atmospheric warming until the end of the century. Daily SMB and its components, including snowfall and rainfall, surface melt, refreezing, retention and runoff, are further statistically downscaled from the output of RACMO2.3 at 11 km to a gridded glacier mask (RGI Consortium, 2017) and digital elevation model (Porter et al., 2018) re-sampled to 500 m horizontal resolution (Figure S1a in Supporting Information S1). Statistical downscaling (Noël et al., 2016) corrects daily melt and runoff for elevation differences on the relatively coarse 11 km grid using modelled elevation gradients, and for underestimated bare ice and volcanic tephra albedo using remotely sensed measurements from the Moderate Resolution Imaging Spectroradiometer (MODIS; Figure S1b in Supporting Information S1; see Section 2.3). We use this new, daily data set, thoroughly evaluated using a combination of in situ (see Section 2.4) and remote sensing measurements (see Section 2.5), that reconstructs the recent  and projects the future (1958-2099) SMB of Icelandic glaciers at 500 m spatial resolution to explore whether the post-2011 mass loss slowdown will further persist under extreme warming in the course of the 21st century.

Regional Atmospheric Climate Model
We use output of the Regional Atmospheric Climate Model (RACMO2.3) over a domain including Greenland, the Canadian Arctic, Svalbard and Iceland at 11 km spatial resolution for the period 1958-2019 (Noël et al., 2015). The model incorporates the dynamical core of the High Resolution Limited Area Model (HIRLAM; Undèn et al., 2002) and the physics package cycle CY33r1 of the European Centre for Medium-Range Weather Forecasts-Integrated Forecast System (ECMWF-IFS, 2008). The model is specifically adapted to represent surface processes in glaciated regions, and includes a 40-layer snow module simulating melt, percolation, retention and refreezing into snow layers, and subsequent surface runoff (Ettema et al., 2010). RACMO2.3 represents drysnow densification (Ligtenberg et al., 2011), drifting snow erosion , and snow albedo based on grain size, cloud optical thickness, solar zenith angle, and impurity content (soot; Van Angelen et al., 2012). Here we estimate the SMB of Icelandic glaciers as the local (kg m −2 ) and the spatially integrated (Gt yr −1 ) sum of, where PR is total precipitation, that is, the sum of snowfall (SF) and rainfall (RA) precipitation, RU is the surface runoff, SU and ER are the total sublimation and the erosion of drifting snow, which are small and do not significantly contribute to SMB. Surface runoff is estimated as the fraction of rain and melt (ME) that is not retained or refrozen (RF) in firn as, The above definition of SMB including 'internal accumulation' from refreezing and retention (RF), is referred to as 'climatic mass balance' in (Cogley et al., 2011).
For the contemporary climate, RACMO2.3 is forced by ERA-40 (1958ERA-40 ( -1978Uppala et al., 2005), ERA-Interim (1979Dee et al., 2011) andERA5 (2019;Hersbach et al., 2020) reanalyses within a 24-grid-cell-wide relaxation zone at the model lateral boundaries. Forcing consists of temperature, pressure, specific humidity, wind speed and direction being prescribed at the 40 model atmospheric levels every 6 hr. Upper atmospheric relaxation is active (Van de Berg & Medley, 2016). Sea ice extent and sea surface temperature are also prescribed from the reanalyses on a 6-hourly basis. The snowpack is initialised in September 1957 using snow temperature and density profiles from the offline Firn Densification Model (FDM; Ligtenberg et al., 2018). Ice albedo is prescribed from the 500 m Moderate-resolution Imaging Spectroradiometer (MODIS) 16-day surface albedo product (MCD43A3) as the lowest 5% annual values averaged for the period 2000-2015, and clipped between 0.30 for dark bare ice and 0.55 for clean bright ice under perennial snow. For future climate, RACMO2.3 is forced by CESM2.

Community Earth System Model
We   . Observed mass change integrated over the three largest Icelandic ice caps using annual in situ SMB measurements (blue) (Aðalgeirsdóttir et al., 2020), and monthly GRACE(-FO) (red) and CryoSat-2 (green) are shown for comparison. CryoSat-2 mass change is offset by 50 Gt for clarity. Uncertainties in the GRACE and CryoSat-2 products are shown as coloured bands (see Section 2.5). Mass loss is converted into mm sea level rise equivalent (right axis  1958-1994, 1995-2010 and 2011-2019. Note that modelled cumulative mass change in (b) does not consider non-surface mass balance processes such as basal melt and calving . except for ice dynamics, that is, excluding calving processes and subsequent glacier thinning and retreat. The model runs at a spatial resolution of 0.9 × 1.25° (∼111 km) and exclusively prescribes land use changes, atmospheric greenhouse gas and aerosol emissions (Eyring et al., 2016). A climate projection from CESM2 under a high-end emission scenario (SSP5-8.5) is used to force RACMO2.3 at 11 km every 6 hr throughout the 21st century . Sea ice extent and sea surface temperature are also prescribed from CESM2 on a 6-hourly basis. The product is further statistically downscaled to 500 m spatial resolution.

Statistical Downscaling
To accurately resolve SMB patterns over the relatively small glaciers of Iceland, the daily output of RACMO2.3 forced by ERA reanalyses (1958( ) and CESM2 (1958( -2099 are statistically downscaled to a 500 m ice mask derived from the Randolph Glacier Inventory (RGI Consortium, 2017) version 6 and the 100 m spatial resolution ArcticDEM (Porter et al., 2018), re-sampled onto a 500 m grid ( Figure S1a in Supporting Information S1). The downscaling procedure corrects individual SMB components (except for total precipitation), that is, primarily surface melt and runoff, for elevation and albedo differences on the relatively coarse model grid at 11 km resolution. These corrections reconstruct individual SMB components on the 500 m topography using daily-specific elevational gradients estimated at 11 km. Runoff is further adjusted using a re-sampled 500 m  Figure S1b in Supporting Information S1). In line with Gascoin et al. (2017) and Gunnarsson et al. (2021), the minimum albedo in MODIS imageries for tephra-covered snow in the accumulation zone is set to 0.30, while for tephra-covered bare ice in the ablation zone, the minimum is set to 0.15. Total precipitation is bilinearly interpolated from the 11 km onto the 500 m grid without additional corrections. The statistical downscaling technique is further described in (Noël et al., 2016).

In Situ Records
We use 1196 in situ surface mass balance measurements covering the period 1991-2019 and collected at ∼65 sites ( Figure 1c) across the main Vatnajökull ice cap . Annual mass balance is estimated from surface height change between consecutive stake measurements in September, while summer and winter balances are estimated from surface height change between consecutive stake measurements in May-September and September-May, respectively, considering the density of the ice/snow that was lost/gained. For a meaningful comparison, modelled SMB is summed over the exact overlapping period in days. Since 'internal accumulation' is not measured by stakes, while it is included in our model SMB definition (Equation 2), we subtract modelled refreezing and retention from the downscaled SMB in the accumulation zone, that is, where SMB > 0 on average for 1958-2019. For local comparison with stake measurements, we selected the downscaled grid cell with the smallest elevation difference among the closest pixel and its eight adjacent neighbours. We rejected one site (HOSP) situated outside the Randolph Glacier Inventory ice mask, and three other sites (BR0-BR1-BR2, see inset map in Figure S2c in Supporting Information S1) situated outside the RACMO2.3 ice mask at 11 km in a sector of Vatnajökull that experiences the highest mass loss rates . The latter three sites are shown for illustration but were not considered in the regression statistics (red dots in Figures S2a, S2b and S2d in Supporting Information S1). To evaluate modelled 2 m air temperature (T2m) in RACMO2.3-ERA, we use 51,833 daily in situ measurements from 27 automatic weather stations (AWS) installed on the four main Icelandic ice caps ( Figure S2e in Supporting Information S1). Measurements cover the period 1994-2019 and span the summer season, but sometimes further extend from mid-April to end of December. We find good agreement between modelled and measured 10.1029/2021GL095697 5 of 11 2 m air temperature (R 2 = 0.84) with a relatively small cold bias (−0.48°C, i.e., mean difference between model and observations) and small RMSE (2.15°C), highlighting the ability of RACMO2.3 to represent the near-surface temperature above Icelandic glaciers ( Figure S2f in Supporting Information S1).

Remote Sensing Records
We  (Figure 1a and Figure S2b in Supporting Information S1), respectively.

Model Uncertainty
The A similar uncertainty is derived from model comparison with seasonal in situ measurements (see Figures S2c and S2d in Supporting Information S1). Compared to stake measurements, the SMB product shows a mean bias of 310 mm w.e. in winter (210 mm w.e. in summer), that is, representative of snowfall accumulation (runoff ablation). Integrated over the long-term accumulation (ablation) area of Icelandic glaciers, this results in an annual SMB uncertainty of 2.2 Gt yr −1 (Figures S2c and S2d in Supporting Information S1). Compared to remotely sensed estimates from GRACE (R 2 = 0.98) and CryoSat-2 (R 2 = 0.93), we find a mean modelled SMB bias of 0.61 and 0.06 Gt per month ( Figure S2b in Supporting Information S1).

Results
The present-day SMB product has been thoroughly evaluated using a combination of in situ and remote sensing measurements, showing excellent agreement (0.81 < R 2 < 0.98; Figures S2a and S2b in Supporting Information S1). In situ measurements consist of 1,196 records of annual, winter and summer balance from ∼65 sites located on the main Vatnajökull ice cap (1991-2019; Figure 1b and Figures S2a, S2c and S2d in Supporting Information S1; Aðalgeirsdóttir et al., 2020;Pálsson et al., 2020). The evaluation of summer and winter balances suggest a small underestimation of winter accumulation (310 mm w.e. yr −1 on average, shown in Figure S2c in Supporting Information S1), that is partly compensated by an underestimation of summer ablation (210 mm w.e. yr −1 , shown in Figure S2d in Supporting Information S1). The ablation underestimation can be partly ascribed to the fact that RACMO2.3 at 11 km resolution does not well capture high melt rates of small, low-lying outlet glaciers, especially in the southeast of Vatnajökull where they almost reach sea level. In addition, the SMB product does not account for non-surface mass balance processes (e.g., basal melt, calving) that are currently estimated to contribute ∼2 Gt yr −1 . To further evaluate glacier-integrated SMB, the downscaled product is also compared to remotely sensed mass change estimates from the Gravity Recovery and Climate Experiment (GRACE) satellite, its Follow On mission (GRACE-FO), as well as CryoSat-2 since 2002 ( Figure 1a and Figure S2b in Supporting Information S1). There is a good agreement between the downscaled SMB product and both satellite estimates, and the resulting total bias is of similar order as the non-surface mass balance contribution . In spite of the good agreement with in situ and satellite gravimetry and altimetry products, it is important to note that (large) local uncertainties in individual SMB components in space and time remain. These uncertainties must be addressed by enhancing the amount of available measurements through continued survey of Icelandic glaciers. Nonetheless, our new present-day SMB product accurately captures both the long-term mass loss trend and the seasonal cycle of winter mass gains and summer mass losses (Figure 1a).
Linear regression shows that variations in the SMB of Icelandic glaciers are mainly driven by fluctuations in meltwater runoff (R = 0.94 in Figure 2a). In turn, runoff varies linearly as a function of spatially averaged glacier 2 m air temperature (T2m; R = 0.92 in Figure 2b). As shown for the Greenland ice sheet (Fausto et al., 2016), the latter relationship reflects that although absorbed shortwave radiation is the main energy source for melt in Iceland, fluctuations in downward longwave radiation and sensible heat flux following the advection of cold/warm air masses explain most of the interannual variability in 2 m air temperature above glaciers, and thus mass change. As a third connection, we find that modelled Icelandic glacier 2 m air temperature strongly correlates with sea surface temperature (SST) measured by ERA reanalyses in the North Atlantic Ocean (R = 0.88 for the period 1958-2010 in Figure 2c), notably in a region south of Greenland called the Blue Blob (Josey et al., 2017).
Here, we focus on the northern part of this cooling area, referred to as Northern Blue Blob (NBB), fringing the southern tip of Greenland and Iceland (blue contour in Figure 3a). As a result, SMB and runoff from Icelandic glaciers also strongly correlate with measured sea surface temperature in the NBB area ( Figure S4a in Supporting Information S1). Surface temperature in the NBB remained almost unchanged until the mid-1990s (grey line in Figure 4b) and thereafter started to increase simultaneously and at the same rate as Iceland glacier 2 m air temperature (+0.07°C yr −1 ; orange line in Figure 4b). After 2011, cooling in the NBB (locally −0.15°C yr −1 in Figure 3a) strongly reduced atmospheric warming above Icelandic glaciers (Figure 4b), hence decreasing meltwater runoff (Figure 1c). We conclude that the post-2011 development of the Blue Blob resulted in the slowdown of glacier mass loss in Iceland (Figure 1a).
To predict the influence of the North Atlantic Ocean temperature evolution on the glacier mass balance throughout the 21st century, we use downscaled output of RACMO2.3-CESM2 under a high-end warming scenario (Figure 4a). Using such a scenario enables us to explore the glacial response to the most extreme future atmospheric and oceanic temperature forcing. First we ascertained that RACMO2.3-CESM2 reproduces (1) Icelandic glaciers to be in approximate mass balance before 1995, (2) the onset of mass loss in the period 1995-2010, and (3) the post-2011 slowdown (Figure 4a). In addition, we confirmed that RACMO2.3-CESM2 realistically captures the spatial correlations between sea surface temperature in the NBB and Iceland glacier 2 m air temperature (Figures S3a and S3b in Supporting Information S1). It is important to note that CESM2 does not assimilate observations, as opposed to ERA reanalyses, and interannual variability will consequently differ between the two forcings. Therefore, it is to be expected that both spatial correlations ( Figures S3a, S3b and S4a, S4b in Supporting Information S1) and modelled SMB from RACMO2.3-CESM2 and RACMO2.3-ERA are slightly different (∼4 Gt on average for 1958-2019; Figure 4a). Nevertheless, both products show similar correlation patterns, SMB variability and trend (−0.2 Gt yr −2 for 1958-2019). Under the SSP5-8.5 scenario, RACMO2.3-CESM2 projects sustained cooling in the NBB area from 2011 to 2052 (Figure 3b). By the mid-2050s, the NBB area is predicted to cool by ∼3°C compared to 1958-2010 (cyan dots in Figure 2c and black line in Figure 4b), resulting in a ∼2°C decrease in glacier 2 m air temperature (R = 0.89; cyan dots in Figure 2c and red line in Figure 4b).

Discussion and Conclusions
We interpret the asymmetric cooling between ocean and atmospheric temperature as the superposition of a background warming trend with a superimposed regional perturbation (Figure 2c). This mechanism is illustrated by idealised linear temperature regressions in Figure S5 in Supporting Information S1. Since the mid-1990s, Icelandic glaciers and the North Atlantic Ocean have experienced a similar background warming (red line in Figure  S5a in Supporting Information S1). Superimposed on this, the North Atlantic is currently undergoing and will undergo a pronounced regional cooling in the Blue Blob area from 2011 until the mid-2050s (blue line in Figure  S5a in Supporting Information S1). As the regional temperature decline exceeds the background warming, this results in a net surface cooling in the NBB area (black line in Figure S5a in Supporting Information S1). This local oceanic cooling is currently observed in the North Atlantic (Figure 3a), and has been linked to the development of the Warming Hole, that is, a long-term (multi-decadal) cooling of the North Atlantic Ocean that has been ascribed to a potential weakening of the Atlantic Meridional Overturning Circulation (AMOC) (Caesar et al., 2018;Keil et al., 2020;Rahmstorf et al., 2015;Zhu & Liu, 2020). While the significance of the long-term Warming Hole cooling and the role played by the AMOC remain uncertain (Fox-Kemper et al., 2021), the strong cooling recently observed in the Blue Blob is strengthened by regional atmospheric cooling resulting from interannual variability of the North Atlantic Oscillation (Josey et al., 2017). Applying a linear fit to the data suggests that Icelandic glaciers experience only 60% of the NBB cooling signal. We interpret this as atmospheric cold content being partly dissipated by the advection of air masses between the NBB and Iceland (orange line in Figure S5a in Supporting Information S1). Consequently, the decline in glacier 2 m air temperature is approximately half that of sea surface temperature in the NBB (Figure 4b), explaining the asymmetrical (z-shaped) correlation (cyan dots in Figure 2c and Figure S5b in Supporting Information S1). After the mid-2050s, North Atlantic cooling is projected to halt in CESM2, so that the background warming signal prevails again in the NBB (Figure 3c). This results in a parallel increase of sea surface and Iceland glacier 2 m air temperature by ∼3.5°C until 2100 (R = 0.96; red dots in Figures 2c and 4b). In line with our results, Gervais et al. (2018) projected a persistent cooling in the NBB until the 2050s, followed by a temperature increase under a high-end warming scenario (see their Figure 4a). However, they do not identify the driving processes of the latter warming. We hypothesise that potential mechanisms include reduced cold water and sea ice transport from the Labrador Sea and Fram Strait into the NBB area following atmospheric warming in the second half of the century.
Based on past and projected sea surface temperature fluctuations in the NBB (black line in Figure 4b), we identify three different stages for the mass balance of Icelandic glaciers under a high-end emission scenario. In 1958-2010 (grey shade in Figure 4), the 2 m air temperature over Icelandic glaciers increases according to the background warming (Figure 3b), enhancing runoff ( Figure 4a) and bringing Icelandic glaciers from approximate mass balance  to mass loss (1995-2010; grey line in Figure 4c). In 2011-2052 (cyan shade in Figure 4), Icelandic glacier mass loss slows down, the result of pronounced cooling in the NBB that reduces runoff and even returns glaciers to a state of approximate mass balance in the late 2040s (cyan line in Figure 4c). From the mid-2050s onward, warming is projected to prevail again in the NBB, enhancing runoff and resuming mass loss (red line in Figure 4c). As a result, Icelandic glaciers could lose about 30% of their total volume by 2100 ( Figure 4c). If, as we project, the slowdown of mass loss persists until the 2050s, Icelandic glaciers could preserve ∼200 Gt of ice (or 5% of the total volume) by the mid-2050s relative to a trajectory of sustained mass loss at the 1995-2010 rate (dashed line in Figure 4c). Under a high-end warming scenario, mass loss accelerates again in the mid-2050s (−14.0 ± 2.5 Gt yr −1 for 2053-2099). At the latter mass loss rate, and without additional feedbacks, Icelandic glaciers could completely vanish by the end of the 23rd century.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.

Data Availability Statement
Data required to reproduce the figures and tables presented in the study are available on Zenodo (https://doi. org/10.5281/zenodo.5836476). Data include (1)