Spatially heterogeneous effect of climate warming on the Arctic land ice

Global warming has already substantially altered the Arctic cryosphere. Due to the Arctic warming amplification, the temperature is increasing more strongly, leading to pervasive changes in this area. Recent years were notably marked by melt records over the Greenland Ice Sheet, while other regions such as Svalbard seem to remain less influenced. This raises the question of the current state of the Greenland Ice Sheet and the various ice caps in the Arctic for which few studies are available. Here, we run the regional climate model (RCM) Modèle Atmosphérique Régional (MAR) at a resolution of 6 km over four different domains covering all Arctic land ice to produce a unified surface mass balance product from 1950 to the present day. We also compare our results to large-scale indices to better understand the heterogeneity of the evolutions across the Arctic and their links to recent climate change. We find a sharp decrease of surface mass balance (SMB) over the western Arctic (Canada and Greenland) in relationship with the atmospheric blocking situations that have become more frequent in summer, resulting in a 41 % increase of the melt rate since 1950. This increase is not seen over the Russian Arctic permanent ice areas, where melt rates have increased by only 3 % on average, illustrating a heterogeneity in the Arctic SMB response to global warming.

1 Introduction

The warming amplification of the Arctic has led to an average temperature rise of +3.8 ^∘C poleward of 66.5^∘N since 1979, 4 times larger than the global average (Rantanen et al., 2022). While this warming contributes to a higher melting rate of glaciers and ice caps (e.g., Fettweis et al., 2017; Noël et al., 2018), it has also raised the atmospheric humidity, leading to more solid precipitation in winter (Przybylak, 2002; Førland et al., 2002). In combination with large-scale atmospheric circulation variations, changes in average melting and precipitation rates modify the surface mass balance (SMB) of the Arctic land ice, i.e., the Greenland Ice Sheet, Arctic ice caps and major peripheral glaciers.

The SMB is the difference between the total amount of precipitation (solid and liquid) plus condensation or riming and the ablation by meltwater runoff and evaporation or sublimation. It is a component of the total ice mass budget of permanent ice areas together with the ice discharge driven by the ice dynamics. Note that the definition we use for SMB was formerly used for climatic mass balance in Cogley et al. (2011). However, the Arctic SMB is more sensitive to seasonal climate variations, and its importance in the total mass budget is expected to increase relative to the ice discharge, at least over the Greenland Ice Sheet (Fürst et al., 2015). Furthermore, the combined Arctic permanent ice areas (excluding the Greenland Ice Sheet) are the major contributors to sea level rise after the ice sheets (Gardner et al., 2013; Moon et al., 2020).

While the warming trend is global, the different studies carried out over the Arctic indicate a regional heterogeneity in the response of SMB to the climate of the last decade. The higher frequency of blocking anticyclonic events has increased the summer melt rate over the Greenland Ice Sheet and the Canadian ice caps (Fettweis et al., 2013; Lenaerts et al., 2013; Noël et al., 2018; Fettweis et al., 2017; Topál et al., 2022; Rajewicz and Marshall, 2014). On the contrary, recent North Atlantic cooling has decreased glacier mass loss rates in Iceland (Noël et al., 2022). In Svalbard, most studies indicate a significant negative mass balance trend in recent decades (e.g., Schuler et al., 2020; van Pelt et al., 2019), although Lang et al. (2015) found a stable mass balance instead.

Figure 1 MAR domains used over the Arctic (solid white boxes) and integration sub-regions for analysis (dashed boxes). AWS locations are shown with a dot colored as a function of the years of data they provide.

High-resolution dynamical downscaling has enhanced the estimations of SMB across the Arctic by providing continuous results in space and time compared to in situ observations (and satellite data). However, a unified estimate is still lacking over all the permanent land ice areas of the Arctic using the same method. Moreover, SMB estimates over the Russian High Arctic remain very scarce. Here, we present the results from a series of dynamical downscaling simulations at high resolution (6 km), covering all the Arctic regions with permanent Arctic land ice (Baffin, Devon, Ellesmere, Iceland, Svalbard, Greenland, Franz Joseph Land, Novaya Zemlya and the Russian High Arctic islands) using the Modèle Atmosphérique Régional (MAR) regional climate model (RCM). The aims of the study are (1) to present a unified SMB product derived from the same method over all the Arctic and (2) to highlight the links between SMB changes over different regions and general climate patterns.

2 Methods

2.1 MAR

MAR is a 3D atmosphere–snowpack RCM initially designed for polar regions (Gallée and Schayes, 1994). It has been used in multiple studies and has proven to be reliable in reconstructing the recent SMB changes over the Greenland (Fettweis et al., 2017, 2020) and Antarctic (Agosta et al., 2019) ice sheets or smaller ice caps (Svalbard, Lang et al., 2015).

MAR resolves the primitive equations using the hydrostatic approximation and has a vertical sigma coordinate system. MAR also includes the 1D surface scheme SISVAT (Soil Ice Snow Vegetation Atmosphere Transfer; De Ridder and Gallée, 1998; Ridder and Schayes, 1997; Gallée and Duynkerke, 1997; Gallée et al., 2001; Lefebre et al., 2003) which describes the surface properties and their evolution through their interactions with the atmosphere. The snow and ice module of SISVAT describes the snowpack metamorphism and properties (such as temperature, liquid water content and density) of the first 20 m of permanent ice areas divided into 30 layers of snow, firn or ice. Since MAR is not coupled here with an ice sheet model, the topography and ice extent are fixed in the model throughout the entirety of the simulations. Pixels are considered to be ice covered only if they have at least 50 % of their area covered by ice.

In this study, MAR version 3.11.5 (hereafter MARv3.11.5) is used to reconstruct SMB changes over the Arctic ice caps and ice sheet. The improvements of this version are described in Kittel et al. (2021). A general summary of the modules and schemes used in MAR can also be found in Fettweis et al. (2017). Figure 1 presents the four integration domains (without the relaxation zone) used to run MAR over all the permanent Arctic ice areas at a 6 km horizontal resolution using 24 vertical layers in the atmosphere, with the first level at 2 m above surface. We used four different integration domains in order to reduce the computational cost of a 70-year-long simulation, enabling such a high horizontal resolution. The model parameters and setup are kept the same over all the domains. While a 6 km resolution might be too low to fully resolve the elevation of the smaller ice areas, the average hypsometry of the model grid ice pixels remains close to observations in every sub-region (see Fig. S1 in the Supplement). The biggest discrepancy can be seen in the Franz Joseph sub-region, where our grid, on average, underestimates the real ice elevation, though it remains a relatively small bias and most likely does not have a major influence on the results. Even though the hypsometries do agree, the model resolution can still affect the surface mass balance in strong-topographic-variations area, affecting shading, wind drift and precipitation.

2.2 Reanalysis

The ERA-5 reanalysis (Hersbach et al., 2020) is used as a forcing field to prescribe MAR boundary conditions every 6 h at each vertical level and over the ocean (temperature, u and v components of the wind, humidity, surface pressure, sea ice concentration and sea surface temperature). We chose ERA-5 because it has the advantage of being continuous from 1950 to the present day, and it is performing well over Greenland (Delhasse et al., 2020).

2.3 Data and evaluation methods

MAR has been used often over Greenland (e.g., Fettweis et al., 2020; Lambin et al., 2022) and Svalbard (Lang et al., 2015) but less frequently over Canada, Iceland or the Russian Arctic. As our study deals with three new MAR domains, more attention is given to the evaluation of the results against field observations. First, the simulations are evaluated for their performance in reproducing real atmospheric conditions (in particular the 2 m temperature – T2m; 2 m pressure – P2m; and wind speed – WS). Then, the reconstructed SMB is compared to the few observations available, namely observations from satellite altimetry (from Hugonnet et al., 2021) to compare the regionally integrated SMB from 2000 to 2020 over land-terminating glaciers, along with the SMB dataset from Machguth et al. (2016) over the Greenland Ice Sheet, available on the PROMICE website.

2.3.1 Evaluation of the atmosphere

Over the different domains, 102 automatic weather stations (AWSs) were used to evaluate the MAR simulations. The localization of the AWSs is shown in Fig. 1. Daily average values were used to compare observations to MAR simulations. For the modeled values, daily means were extracted as a distance-weighted mean between the four nearest MAR pixels. To avoid bias coming from ocean pixels (where the SST is prescribed to MAR from ERA5), only land MAR pixels were considered for the evaluation. Finally, the AWSs with an elevation difference of more than 200 m with the four-nearest-pixels average were excluded to avoid artificial biases driven by the elevation difference (16 excluded in total). Mean bias, root mean squared error (RMSE), centered root mean squared error (CRMSE) and correlation (r) between observed and modeled values were computed.

2.3.2 SMB evaluation

As precipitation and snow surface processes are the most challenging variables to represent in climate models, large biases can arise between models and observations when simulating the SMB. It is crucial to evaluate the modeled SMB over the different regions, although direct observations remain scarce, especially over the Russian Arctic.

Hugonnet et al. (2021) developed a global product of glacier elevation change from 2000 to 2019 using NASA's Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER). Their glacier mass change product consists of monthly mass loss estimates integrated over sub-regions of the Randolph Glacier Inventory (RGI). It contains mass balance (MB) estimates for all the land ice in Canada, Iceland, Svalbard, the Russian archipelagoes and the Greenland periphery. Because the MB is the difference between the SMB and the dynamical iceberg discharge (in the case of marine-terminating glaciers), we selected data for only land-terminating glaciers using their classification in the RGI 6.0. Annual modeled SMB was then integrated over all the glaciers to be compared with MB estimates.

This altimetry dataset is useful for evaluating the SMB over large remote regions of our study, where very sparse in situ observations are available. We use the in situ SMB dataset from Machguth et al. (2016) to evaluate MAR, as done in Fettweis et al. (2020) over Greenland, as the mass loss by iceberg discharge over the Greenland Ice Sheet is significant compared to over smaller Arctic ice caps. This dataset contains historical SMB measurements from more than 3000 stakes over the ice sheet. It is quite different from the evaluation using the Hugonnet et al. (2021) MB product (annual spatially integrated data) so the results will not be intercomparable but will give another estimate of the performance of MAR.

Moreover, we evaluated the annual modeled SMB on given glaciers using the World Glacier Monitoring Service (WGMS) dataset. The same method was applied as above, using RGI 6.0 glacier geometries intersected with MAR pixels to spatially integrate the specific SMB over a given glacier. Whilst the spatial coverage of the dataset is low, it has the advantage of comparing the SMB measurements directly as opposed to the altimetry product.

3 Evaluation

3.1 Climate evaluation

Table 1 presents the resulting mean bias, CRMSE and r over all regions for the near-surface temperature, pressure and wind speed. The results are computed annually and seasonally (JJA for summer and DJF for winter) and for each MAR domain. As the datasets for Svalbard and the Russian archipelagoes are inside the same domain, we separated them for the evaluation as the observational datasets are different. The height measurement for wind speed was not always available. We then use the 2 m wind speed from MAR for the comparison, potentially leading to inherent biases.

Table 1 Evaluation results (bias, centered root mean squared error (CRMSE) and correlation coefficient (r)) over all regions, annually and seasonally. Standard deviation is provided as a ± value. There are 19 AWSs in Canada, 6 in Iceland, 69 in Greenland, 4 in Svalbard and 7 in the Russian Arctic.

The correlation coefficient between MAR and the observed 2 m pressure (P2m) is mostly larger than 0.9 over all regions. The high negative bias over Svalbard and Greenland is imputable to an often high altitude difference between AWSs and MAR pixels (see Table S1 for a comparison of altitude and temperature correction over Svalbard). This difference does not influence the correlation, which is the only relevant statistical value concerning pressure. The 2 m temperature is also reproduced very well in each domain (the annual correlation coefficient is artificially driven up because of the seasonal cycle). There is, however, a general negative bias compared to observations. Moreover, the temperature is better reproduced in winter than in summer (r=0.91 vs. 0.81), because temperature variability is lower as the near-surface temperature is close to 0 ^∘C most of the time and less driven by the general circulation dynamics than in winter. Finally, the modeled 2 m wind speed has a bias lower than ±1.5 m s⁻¹, and the performances of the MAR reconstruction are homogeneous over all domains. However, wind speed observations are particularly sensitive to local site effects, which are not resolved at a resolution of 6 km (as seen by the correlation values lower than 0.7; Lambin et al., 2022). Moreover, we do not have the information regarding the height of the measurements, which can also influence the comparison. The quality of the reanalysis products (ERA-5) depends largely on the number of observations that were assimilated. Because our study goes up to 1950, it is worth evaluating the precision as a function of time: the further back we go, the fewer the observations that will be available. For example, ERA-5 has been shown to perform less prior to 1979 above the Antarctic because of the scarcity of satellite observations over the continent (Marshall et al., 2022).

Figure 2 Time evolution of the near-surface pressure (P2m) and temperature (T2m) 10-year correlation coefficient distribution among AWSs between daily observed and modeled values. Boxes show the median and quartiles of the correlation distribution amongst stations, and the whisker extent shows the rest of the distribution outside of outliers (diamonds).

While we could expect a better agreement after 1979 in our ERA5 forced MAR simulation due to the assimilation of satellite data in the reanalysis, there is not a significant evolution of the correlation coefficient as a function of time for the P2m (Fig. 2). The latter is constant at approximately r=0.99 from 1950 to the present day. The analysis nevertheless reveals a slight increase in the correlation coefficient of the 2 m temperature (from 0.96 in 1955 to 0.97 in 1985) as a result of satellite datasets assimilation (in particular, sea ice cover (SIC) and sea surface temperature (SST)) after 1979, which strongly influences the reconstructions elsewhere, mostly where observational data were scarce (Marshall et al., 2022). However, the good number of older observations in the Arctic (compared to in the Southern Hemisphere) explains the good performance of MAR forced by ERA-5 before 1979 (Hersbach et al., 2020) and even before 1957, the International Geophysical Year (see Table S2).

Figure 3 Statistical distribution of average yearly SMB values over land-terminating glaciers according to the RGI 6.0 between 2000 and 2020 for modeled values (in orange) and observed values (MB estimates from Hugonnet et al., 2021, in blue). Note that Greenland (grd) only includes peripheral glaciers. Boxes show the median and quartiles of the distribution amongst years, and the whisker extent shows the rest of the distribution outside of outliers (diamonds).

3.2 SMB

Figure 3 shows the statistical distribution of annual modeled SMB values (mod) and MB satellite observation (obs) estimates over land-terminating glaciers in all sub-regions for the period 2000–2020. There are some biases in the annual mean values, which are positive over Svalbard, the Greenland periphery and Ellesmere Island and negative over Baffin, Devon and Iceland. The main difference is in the variability of the values, where modeled interannual variability is systematically higher than the observed one. This could be related to the lower interannual variability of altimetry products because of snowpack densification and ice dynamics (Li et al., 2023).

Land-terminating glaciers represent only a small fraction (10 % when accounting for the Greenland Ice Sheet, 43 % when not) of all the ice areas studied here. The main bias of this evaluation comes with the integration of the 6 km MAR pixels over small glaciers (especially with small ice tongues) with a strong spatial SMB gradient or a notable sensitivity to site effects. Finally, while the RGI is generally precise in the classification of land- and marine-terminating glaciers, it is sometimes less accurate (as in northern Svalbard, for example), which could explain the slight positive bias of the simulated SMB values as some ice discharge would be included in the observation dataset.

Over the Greenland Ice Sheet, the evaluation using the PROMICE dataset yields a correlation of 0.93 between the model and the observations. The average bias is +0.03 m w.e yr⁻¹ (for an average observational value of −0.86 m w.e yr⁻¹), and the RMSE is 0.43 m w.e yr⁻¹.

Figure 4 Modeled and observed annual mean specific mass balance over different glaciers of a given region using the WGMS dataset, with n being the number of observations (Greenland periphery and Franz Joseph Land are not included due to having no observations).

Figure 5 Annual (thin line) and 20-year running mean (thick line) of the annual integrated SMB (blue), runoff (RU, orange) and snowfall (SF, green) over (a) Greenland, (b) Iceland, (c) Svalbard, (d) Franz Joseph Land, (e) Novaya Zemlya, (f) the Russian High Arctic islands, (g) Baffin, (h) Devon and (i) Ellesmere.

Figure 4 shows the statistical distribution of modeled and WGMS SMB observations amongst measurements for every sub-region where data were available. The variabilities of the observations are closer to modeled variabilities than in the case of the altimetry product, which is in line with their lower interannual variability. The main bias of this comparison (added to the 6 km MAR resolution as mentioned above) comes from the low spatial extent of the WGMS dataset for some sub-regions. As such, a detailed list of all measurements can be found in the Supplement (Fig. S2). Finally, a more refined comparison between point-stake measurements and altitudinally downscaled modeled SMB over Svalbard is also available in the Supplement (Fig. S3). The values of RMSE in this study for this point-stake evaluation (0.73 m w.e. yr⁻¹) are a bit larger than those of Østby et al. (2017) (0.59 m w.e. yr⁻1) and of van Pelt et al. (2019) (0.43 m w.e. yr⁻¹), though the latter two studies calibrated their models to reduce discrepancies with the stake data.

4 Results

Our simulations show that the Arctic experiences an overall yearly SMB anomaly of −96.4 Gt yr⁻¹ over 2000–2020 compared to the reference period of 1950–1970. This value becomes even more negative when considering the recent past evolution, with an anomaly of −154 Gt yr⁻¹ between 1975–1995 and 2000–2020. This total SMB decrease is mainly driven by Greenland (due to it being by far the largest ice area). However, Greenland runoff has increased by 35 % between 1975–1995 and 2000–2020 but has, on average, increased by 45 % over the other regions. This difference implies that there is a clear interest in analyzing the different Arctic sub-regions independently to better identify the driving processes involved.

Table 2 Regional averages and variability of the SMB, runoff (RU) and snowfall (SF) integrated over permanent ice areas for different time periods (in Gt yr⁻¹).

4.1 Integrated SMB changes

Figure 5 shows the evolution of the integrated SMB, SF and RU over all subregions. The Baffin Island, Devon and Ellesmere Island ice caps and glaciers have been losing mass since 1950. Over Baffin Island, this has been accelerating in recent years, with the SMB going from −22.1 Gt yr⁻¹ between 1950 and 1970 to −33 Gt yr⁻¹ from 2000 to 2020 (comparable to results from Noël et al., 2022). The snowfall has remained stable across the whole period, while the runoff has increased significantly (from 39.6 to 52.8 Gt yr⁻¹). Further north, the Devon ice cap has seen roughly the same evolution as that of Baffin Island. The SMB decreased from −4 Gt yr⁻¹ over 1950–1970 to −6.5 Gt yr⁻¹ over 2000–2020 as a consequence of higher runoff (+2.4 Gt yr⁻¹) but stable snowfall. The same evolution also occurred over Ellesmere Island, where the 30 % increase in runoff lead to a decrease in SMB (from -9.3 to −16.8 Gt yr⁻¹) over 2000–2020 compared to 1950–1970.

Similarly, the SMB decreased over Greenland, Iceland and Svalbard over 2000–2020 compared to the period before 1970. The strong increase in runoff (anomaly of +119.1 Gt yr⁻¹) over the Greenland Ice Sheet despite higher snowfall (+42.5 Gt yr⁻¹) has resulted in a lower SMB (from 343.7 to 267 Gt yr⁻¹.). Over Iceland, the increase in runoff is not compensated for at all by snowfall that remained stable, leading to an SMB decrease of 1.4 Gt yr⁻¹. Over Svalbard, the net SMB was, on average, positive (1.7 Gt yr⁻¹) before 1970 but negative (−0.8 Gt yr⁻¹) after 2000 as a result of an increase in runoff (+8.3 Gt yr⁻¹).

On the other side of the Arctic, the SMB increased over the Franz Joseph Land archipelago, Novaya Zemlya and the Russian High Arctic islands over 2000–2020 compared to 1950–1970. As a result of higher snowfall (+0.9 Gt yr⁻¹) and a stable RU, the SMB is now higher over the Franz Joseph Land archipelago. It is also higher over Novaya Zemlya (−5.2 to −4.2 Gt yr⁻¹) for the same reasons. Finally, over the Russian High Arctic islands, the SMB increased steadily from −1.8 to 0.6 Gt yr⁻¹ because of both an increase in SF (5.1 to 6.5 Gt yr⁻¹) and a decrease in RU (8.5 to 8 Gt yr⁻¹). It is the only region where the RU decreased overall in the simulation period.

The previous paragraphs suggest a similar temporal evolution for different SMB components and/or regions. We now present normalized values of the 20-year running mean SMB, snowfall and runoff over all regions for a better intercomparison regardless of the size of the different regions (Fig. 6.).

Figure 6 The 20-year running mean of the normalized time series (mean-subtracted and divided by standard deviation) of (a) SMB, (b) RU and (c) SF over all regions.

Snowfall increased everywhere from 1950. The Russian archipelagos (Franz Joseph, Novaya Zemlya and the Russian High Arctic archipelagos) have seen the largest relative increase of snowfall, with an acceleration since 1995. To a lesser extent, this can also be observed for Svalbard and Greenland. However, our results suggest that this reached a peak around 1985 over Iceland and around 1995 for the Canadian regions (Ellesmere, Devon and Baffin islands).

Svalbard excepted, all the regions experienced a large decrease in runoff before a significant increase. Runoff decreased until 1975 over Greenland and the Canadian Arctic and until 1985 for the Russian archipelagos. On the contrary, the runoff is steadily increasing throughout the whole period over Svalbard. While Iceland, Greenland and the Canadian Arctic experienced an increase in the runoff since 1975, it is clear that the climatological average runoff increase accelerated unequivocally in all regions from 2000.

Finally, the SMB evolution can be divided into three main periods over all regions. A first period where it increased as runoff decreased, then a second period with a stabilization (increase in both runoff and snowfall) and then a third where the strong increase in runoff led to a large decrease in SMB. It is important to mention that, Svalbard excepted, all regions had a higher SMB over 1975–1995 than over 1950–1970 or 2000–2020 (Table 2).

Table 3 Correlation coefficient between annual values of SMB, RU and SF over all sub-regions.

The SMB evolution is relatively associated with the runoff evolution, only of the opposite sign (see Table 3). This indicates that, although snowfall has increased, melt and runoff variations are the main drivers of the recent SMB over the whole Arctic. Greenland, Iceland and the Canadian Arctic follow the same pattern, with a slight increase from 1960 to 1975, followed by a decrease from 1975 to 2000 that accelerates afterward. The Russian archipelagos, on the other hand, experienced a large increase from 1960 up to 1980, followed by a stabilization between 1980 and 2000 and a slight decrease afterward. Only Svalbard stands out as having a relatively stable SMB (increase in both runoff and snowfall compensating for each other) throughout the whole simulation period, though we still find a SMB linear trend of −0.04 Gt yr⁻² (this is in line with some previous results (e.g., Lang et al., 2015), although there remain significant discrepancies between studies over Svalbard (see Table 4)).

Noël et al. (2018)Noël et al. (2018)Fettweis et al. (2020)Lenaerts et al. (2013)Noël et al. (2022)Noël et al. (2022)Noël et al. (2020)Noël et al. (2020)Radić and Hock (2011)van Pelt et al. (2019)Aas et al. (2016)Lang et al. (2015)

Table 4 Regionally integrated SMB comparisons between this study and past studies. No past estimates were found for the Russian Arctic.

Comparing our results to previous studies in Table 4, we found integrated annual SMB values that were close to what Noël et al. (2018) and Lenaerts et al. (2013) found over the Canadian Arctic. Our results are also comparable to what Fettweis et al. (2020) and Noël et al. (2022) found over Greenland and Iceland, respectively. One main discrepancy concerns Svalbard, where our results suggest a significantly higher integrated SMB compared to Radić and Hock (2011) and Aas et al. (2016). Our results are also higher than those of Noël et al. (2020) for the 2013–2018 period but are lower for the 1958–1985 period. As said above, uncertainties are still large over Svalbard, and not all studies agree, though some progress has been made towards identifying a clear tendency (Schuler et al., 2020).

4.2 Spatial tendencies

Generally, glaciers, ice caps and ice sheets tend to see their equilibrium lines (annual SMB equal to zero) rise because of global warming. This tendency is often driven by the increase in surface melt at lower altitudes. This phenomenon can be seen in Greenland where the ablation zone experienced an SMB decrease of up to −350 kg m⁻² yr⁻¹ on average between 1960 and 2000 (Fig. 7a). At the same time, the northeastern interior of the Greenland Ice Sheet has experienced an SMB increase of +50 kg m⁻² yr⁻¹ as a result of more snowfall (see Fig. S4).

Figure 7 Annual SMB anomalies between the 1990–2020 and 1950–1979 periods over (a) the whole Arctic, (b) Novaya Zemlya, (c) Svalbard, (d) Devon, (e) Iceland, (f) the Russian High Arctic islands, (g) Ellesmere, (h) Baffin and (i) Franz Joseph Land. Hashed areas denote where the anomaly has a low significance value with regard to its variance (using Student's t test with 90 % p value). The equilibrium line between ablation and accumulation areas for the 1990–2020 period is shown with a dashed purple line.

This tendency is not present over the southern Canadian ice caps (Devon, Baffin), where the SMB decreased nearly everywhere by at least 100 kg m⁻² yr⁻¹. In Iceland, the Vatnajökull ice cap has seen an overall decrease in SMB, except in its southern part. Looking at Svalbard, there is a difference between the southern part of the region where SMB decreased significantly in the ablation zone and the northern, higher, and colder part of the region where SMB increased by more than 200 kg m⁻² yr⁻¹ due to larger snowfall (see Fig. S4). Finally, the Russian archipelagos (Franz Joseph, Novaya Zemlya and the Russian High Arctic) experienced a general increase of SMB over nearly their whole surface, being an ablation or accumulation area.

Overall, there is a difference in the SMB evolution between the western part of the Arctic (Canada and Greenland) where the SMB decreased and the eastern part of the Arctic (Svalbard and the Russian archipelagoes) where the SMB increased after 2000 compared to the period before 1970. This difference in tendency is very clear between 1950 and 1979 and remains during the recent period, though it is less pronounced.

5 Discussion

5.1 Correlation to large-scale indices

Between the dry center of Greenland and the marine Russian archipelagos, the wide variety of climates across the Arctic cryosphere may mean that its response to climate change is not homogeneous spatially. This can be already shown by comparing the climate of the recent past in Greenland (Fettweis et al., 2017) and, for example, Iceland (Noël et al., 2022). In the latter region, it has been shown that the North Atlantic cooling has contributed to stabilizing the SMB of Iceland since 2010, while over Greenland, melt rates were increased by the recurring atmospheric blocking situation gauged by negative North Atlantic Oscillation (NAO) conditions.

Figure 8 Correlation over 1950–2020 of integrated summer (JJA) ME (a) and annual SF (b) over all sub-regions and annual large-scale atmospheric and oceanic indices. AMO: Atlantic Multi-decadal Oscillation, GBI: Greenland Blocking Index, NAO: North Atlantic Oscillation, POL: Polar–Eurasian pattern, EAWR: East Atlantic–western Russia index, SCA: Scandinavian pattern, WP: West Pacific pattern, AO: Arctic Oscillation, PNA: Pacific North American index, PDO: Pacific Decadal Oscillation, NP: North Pacific index, SST: annual average sea surface temperature over 70^∘ N, SIC: annual average sea ice concentration over 70^∘ N.

With the same idea of linking SMB variations to large-scale changes, we selected a wide variety of atmospheric indices, averaged over the whole year, to compare with the annual time series of SMB variables for every Arctic region. Figure 7 shows the correlation of annually averaged atmospheric indices to summer (JJA) melt (a) and snowfall (b), the main drivers of SMB over the different regions. Two more oceanic indices were also added, namely the annual average sea surface temperature over 70^∘ N (SST) and the annual average sea ice concentration (SIC) over 70^∘ N. Overall, a lot of indices do not correlate with the melt rates or snowfall rates.

We see, however, a strong correlation between the melt rates in the western part of the Arctic (Greenland and Canada) and the GBI and AMO indices. This has already been observed in the recent past (e.g., Fettweis et al., 2013) for Greenland. It implies that the blocking situation, which increases melt over Greenland, also strongly impacts the Canadian Arctic. We can also observe the anticorrelation between the GBI index and snowfall in Iceland. This might be related to the northerly flow induced by the anticyclonic conditions over Greenland when the GBI is high (with an often low NAO), as has already been suggested by past studies (Matthews et al., 2015; Fettweis et al., 2013; Rajewicz and Marshall, 2014).

On the other side of the Arctic, no such significant correlation to atmospheric indices is found. We see, however, a common correlation (respectively, anticorrelation) between Svalbard's and the Russian archipelagos' snowfall and average Arctic SST (respectively, Arctic SIC). This suggests that a warmer ocean and less ice-covered ocean has likely resulted in higher evaporation and then more snowfall. We found a strong correlation between snowfall and the temperature of the atmosphere around this region (not shown), which also implies higher-saturation water vapor pressure and, further, more precipitation. However, our results do not enable a statement of whether the additional humidity mainly comes from the neighboring ocean, from more remote areas or from a combination of both sources.

Figure 9 Inter-regional correlation of 1950–2020 annual SMB (a), ME (b) and SF (c).

Correlating annual values of SMB between all sub-regions (see Fig. 9a) confirms the existence of two distinct sub-region groups of similar evolutions: Greenland along with the Canadian Arctic on one side, and all of the eastern Arctic from Svalbard to the Russian High Arctic islands on the other. We see again that the SMB correlation is mainly driven by ME. We also see that SF correlates more between the regions of the eastern Arctic than of the western Arctic, noticeably between Novaya Zemlya, Franz Joseph Land and the Russian High Arctic islands.

Figure 10 Cumulative modeled SMB and observed MB from 1990 to 2020 for (a) Greenland (using IMBIE dataset) and from 2000 to 2020 for (b) Iceland, (c) Svalbard, (d) Franz Joseph land, (e) Novaya Zemlya, (f) the Russian High Arctic islands, (g) Baffin, (h) Devon and (i) Ellesmere.

5.2 Mass balance comparison and calving rates

As we studied processes taking place at the surface of permanent ice areas, integrating the SMB spatially does not reflect the total ice mass balance (MB). More specifically, the increase in SMB does not imply an increase in ice mass as altimetry and gravimetry measurements demonstrated that all the regions studied here are still losing mass. Though the scale is much smaller than in the Antarctic, ice calving can make up large proportions of the total ice loss (sometimes called dynamic ice loss) in some Arctic subregions. For example, ice discharge was roughly equal to melting in Greenland between 2008 and 2012 (Enderlin et al., 2014). To assess the importance of calving against SMB over the different regions, we compare our SMB estimates to integrated MB products. Over the Greenland Ice Sheet, we used the IMBIE mass balance dataset (The IMBIE Team, 2020). It consists of mass change measurements from satellite gravimetry and satellite altimetry from 1997 up to 2012. For the other regions, we used the altimetry dataset of (Hugonnet et al., 2021) mentioned previously.

By comparing the cumulative MB and the cumulative SMB (Fig. 10), we can estimate the calved volume over all the sub-regions of the study. Over the Canadian Arctic (Devon, Ellesmere and Baffin), the dynamic ice loss is relatively low (even close to zero in the case of Baffin); thus, the ice mass loss can be considered to be mainly driven by the surface mass balance and then the atmospheric conditions. This low dynamical ice loss can be explained because only a few glaciers are marine terminating in Baffin Island. This is, however, not the case over Ellesmere and Devon islands, where the surface ratio of marine-terminating glaciers is close to 50 %. There, the low dynamical ice loss could be explained by the SST of the waters surrounding the North Arctic Canada, which has not significantly warmed yet compared to atmospheric temperature over the glaciers. Contrarily, in the eastern Arctic, while the SMB continues to be positive over Franz Joseph or the Russian High Arctic islands and has overall increased since 1950 (see Fig. 7), the ice mass is still decreasing rapidly (up to an MB of −5 Gt yr⁻¹ over Franz Joseph Land). It is also the case over Svalbard and Novaya Zemlya, where a relatively constant SMB since the beginning of the 21st century goes along with a steady decrease of the total ice mass. This can be explained by the rapid Arctic Ocean warming that increases the calving rates rapidly, particularly near Svalbard and Novaya Zemlya, where its warming is the most pronounced, with more than 0.8 ^∘C per decade (Li et al., 2022). Note that, while the Greenland Ice Sheet's SMB is positive, lower recent values have resulted in stronger mass loss, as highlighted by Fig. 10a. Greenland aside, the average Arctic ice MB has been −111 Gt yr⁻¹ since 2020, while the average SMB has been −62 Gt yr⁻¹.

6 Conclusions

Considering all the land ice over the Arctic, our simulations reveal that the annual surface mass balance decreased by 120 Gt yr⁻¹ between the period of 1950–1979 and 2000–2020. This overall mass loss has been accelerating by −4 Gt yr⁻² from 2000. It is mainly driven by melt, which has, on average, increased by 21 % since 1950. This melt increase is, however, heterogeneous spatially, with an increase of 41 % for Greenland but only 9 %, on average, over the Russian sub-regions where snowfall accumulation has increased by 28 %. Along with Svalbard, those regions have experienced a general increase of their SMB when looking over the whole simulation period. However, record-low SMBs have been observed everywhere during the last decade, such as in 2020 for all of the eastern Arctic, Devon and Ellesmere or in 2019 for Greenland and Devon.

We have also identified two distinctive sub-region groups (Baffin, Devon, Ellesmere and Greenland and Svalbard, Franz Joseph, Novaya Zemlya and the Russian High Arctic) that seem to have the same links to climatological drivers and that went under a comparable SMB evolution. We have shown that melt is correlated to GBI over Greenland and northern Canada. Snowfall over the latter group seems to be correlated to the average Arctic ocean temperatures, while it is not the case elsewhere. No atmospheric large-scale indices seem to be correlated to its evolution. While these links have been established for the annual mean SF and ME time series, more work remains to be done to understand what is driving the surface mass balance over those two groups. This is especially the case over the Russian Arctic, where only a few studies have been carried out.

Finally, our results suggest rapid changes in the Arctic land ice. While some regions in the Arctic have gained mass at their surface (but are still losing mass taking into account the ice dynamics), these conclusions could be totally different in the years to come. For instance, most recent years were marked by several negative records over the Russian sub-regions. A repeat year of such extreme melting could quickly reverse the trend in these regions and lead to a general loss of surface mass throughout the Arctic. It will therefore be important to continue to study the Arctic land ice and to update these results regularly.