Arctic glaciers record wavier circumpolar winds

Glaciers in the Arctic respond sensitively to climate change, recording the polar amplification of global warming with increasing mass loss. Here, we use glacier mass balances in Svalbard and northern Arctic Canada to categorize tropospheric variability and the associated summer circulation over the Arctic. We establish a link between annual glacier mass balances and their respective atmospheric forcings since 1950 using GRACE/GRACE-FO satellite data (2002–2021), as well as regional climate models and reanalysis data (1950–2019). We find that asynchronous behaviour of mass balance between the regions has become very likely since the early 2000s, exceeding the range of previous decadal variability. Related tropospheric circulation exhibits more meridional patterns, a greater influence of meridional heat advection and a wavier summer circulation. The traceable impact on glacier mass balances emphasizes the importance of dynamic next to thermodynamic climate changes for the future of glacier mass loss, Arctic ecology and societal impacts. Glacial mass is dependent on the balance between melt and snow accumulation, which is impacted by rising Arctic temperatures. Glacier mass balances in Svalbard and northern Canada were asynchronous since the 1990s, related to changes in patterns of atmospheric heat advection.

G lobal warming has disproportionately affected the northern polar latitudes (Arctic amplification) [1][2][3] , and record temperatures have triggered sensitive responses in the mass storage of Arctic glaciers 4 . Glacier mass depends on the balance between accumulation and ablation processes: snowfall, refreeze and retention of meltwater add mass to the glacier system or redistribute it, whereas snowmelt, evaporation and sublimation subtract mass from it. The difference between accumulation and ablation processes is summarized in the surface mass balance (SMB). For marine-terminating glaciers, such as 163 tidewater glaciers in Svalbard (SVA) 5,6 and over 300 in the Canadian Arctic Archipelago 7 , mass is additionally removed by submarine melt and calving (terminal discharge, D). Regional climate models specialized in snow and ice processes, such as MAR3.11 8 or RACMO2.3 9,10 , have been shown to successfully reproduce satellite 11 and in situ observations of glacier and ice sheet mass balance when forced by global reanalysis data 12 such as the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v.5 (ERA5) 13 .
Here, we focus on Arctic Canada North (ACN) and SVA as part of the circumpolar ring of glacier-covered islands and landmasses ( Fig. 1) affected by disproportionate Arctic warming, with their substantial contributions to sea level rise 4 . These glacier regions are separated by roughly 100° of longitude with ice surfaces covering an area of about 34,000 km 2 (SVA) and 105,000 km 2 (ACN) 4 , respectively. Located at about the same latitude (74°-83° N), climatic conditions over both archipelagos are influenced by the same circulation regime characterized by year-round jet-stream dynamics and the winter polar vortex 14 . As ACN is situated east and SVA west of the Greenland Ice Sheet, both regions have been influenced by more frequent blocking conditions in recent years 15,16 .
Temperatures in Svalbard are projected to rise by 6-7 °C until 2100 under the IPCC scenario representative concentration pathway (RCP) 4.5 and up to 10 °C under business-as-usual conditions (RCP 8.5) compared with 1971-2000 17 . Precipitation is said to increase by 45-65% and to contain a progressively lower fraction of snowfall, causing unchanged or more snow storage under RCP 4.5 and less snow for RCP 8.5 over most of Svalbard (50-100% in many places) 17 . Combined with a shorter accumulation season 18 , the overall SMB will become increasingly negative. Dynamic loss of ice through calving accounts for about 21% of the ablation of tidewater glaciers in Svalbard 6,19 .
Temperatures in the Canadian Arctic Archipelago are projected to rise similarly, by about 6.5 °C until 2100 under the moderate RCP 4.5 scenario, causing the loss of about 18% of its glacier volume, mainly through runoff 20 . Between 1996 and 2015, temperatures have increased more in northern than southern Arctic Canada and, in combination with reduced precipitation, drove an acceleration in mass loss of 2.1 ± 0.8 Gt yr -2 in ACN 21 with a marked transition to increased mass loss in 2004 11 .

Glacier mass balances of ACN and SVA
Time-variable gravity field measurements of the GRACE/ GRACE-FO satellite missions 22,23 prove continuous mass loss, with rates amounting to -30.9 ± 8.3 Gt yr -1 (±2σ) in ACN and -13.2 ± 2.3 Gt yr -1 in SVA for 2002-2021 ( Fig. 1a,b). The GRACE/ GRACE-FO twin orbiters record the inter-satellite distance variations caused by the differential gravitational attraction as they pass over Earth features of varying density at slightly different times. Repeated measurements provide a unique monthly time series of gravity field variations, from which mass changes, such as ice loss to the ocean, can be deduced with a typical spatial resolution of 300 km and an accuracy of 2 cm water equivalent 23,24 . The GRACE/ GRACE-FO uncertainties for mass changes in ACN and SVA are shown in Extended Data Fig. 1.
Apart from the long-term mass loss, GRACE/GRACE-FO data reveal that SVA's annual mass-balance anomalies consistently deviate in the opposite direction from those of ACN (Fig. 1c). The annual balances represent the mass change from March to March of the following year (Methods), selected to encompass the full melt season and to centre on the minimum Arctic ice conditions Arctic glaciers record wavier circumpolar winds Ingo Sasgen 1 ✉ , Annette Salles 1,2 , Martin Wegmann 3,4 , Bert Wouters 5,6 , Xavier Fettweis 7 , Brice P. Y. Noël 5

and Christoph Beck 8
Glaciers in the Arctic respond sensitively to climate change, recording the polar amplification of global warming with increasing mass loss. Here, we use glacier mass balances in Svalbard and northern Arctic Canada to categorize tropospheric variability and the associated summer circulation over the Arctic. We establish a link between annual glacier mass balances and their respective atmospheric forcings since 1950 using GRACE/GRACE-FO satellite data (2002-2021), as well as regional climate models and reanalysis data (1950-2019). We find that asynchronous behaviour of mass balance between the regions has become very likely since the early 2000s, exceeding the range of previous decadal variability. Related tropospheric circulation exhibits more meridional patterns, a greater influence of meridional heat advection and a wavier summer circulation. The traceable impact on glacier mass balances emphasizes the importance of dynamic next to thermodynamic climate changes for the future of glacier mass loss, Arctic ecology and societal impacts.
in September. Out of 16 annual mass balances (2003-2016 and 2019-2020) and one biennial mass balance (2017 and 2018; bridging GRACE/GRACE-FO gap, displayed as white dotted line in Fig.  1a,b), only two years are exceptions from this counter-directional mass-balance behaviour.
To examine whether the same characteristics hold true for the time before the GRACE/GRACE-FO era, we reconstruct the SMB of the two regions using the regional climate model MAR3. 11 (1950-2019), six-hourly forced by ERA5 13 , run at a resolution of 15 km with its coupled multilayer snow and ice energy balance model SISVAT 8 .
Cumulative monthly SMB from MAR3.11 typically lies within 2σ error margins of the GRACE/GRACE-FO time series for their common period as shown on Fig. 1a,b (for alternative SMB models, see Extended Data Fig. 2). Note that we account for long-term ice-dynamic changes and the retreat of marine-terminating glaciers on the basis of the difference between the GRACE/GRACE-FO trend and the mean MAR3.11 SMB for 2002-2019, while additionally considering the mass loss of 4.4 ± 0.8 Gt related to the episodic surge loss of Austfonna Glacier in 2013 in Svalbard 25 (Fig. 1b). This estimated discharge D* is then subtracted from the SMB time series   to approximate the complete mass balance (SMB-D*) for ACN and SVA back to 1950, which is the first available year of our MAR3.11 dataset. At this point, we do not aim to fully reconstruct the time evolution of the dynamic discharge for both regions due to our focus on interannual changes, but note the widespread acceleration of D* in ACN 7 and evidence for increased discharge of selected glaciers in SVA 26 .
To isolate the interannual mass variations related to atmospheric dynamics and to remove long-term climate change signals, including the gradual ice-dynamic changes in D* 7,26 , we detrend the monthly time series for the entire period 1950-2019 using a polynomial fit (order n = 3). With the polynomial fit, we remove mass loss trends and accelerations from the time series (orders n = 1 and n = 2, respectively) and improve the recovery of interannual variations, particularly at the beginning and at the end of the period (order n = 3). As before, we then estimate annual mass balances on the basis of the detrended monthly SMB-D* estimates from MAR3.11.

Synchronicity of mass balances in the regions
The annual mass-balance anomalies in SVA and ACN reveal interesting transitions in temporal relationship between the regions (Fig. 2a). While their annual mass-balance anomalies tend to have the same sign from 1960 to the early 1970s and in the mid-1990s, they show asynchronous behaviour between 1975 and 1985 and a distinct transition beginning in the late 1990s, continuing until today (Fig. 2a). Thus, the period up to the year 2002 is characterized by moderately co-varying or opposing mass-balance behaviour in SVA and ACN, with 15 yr running correlations ranging typically between -0.4 and 0.4 (Fig. 2b). Beginning in the year 2002, however, mass-balance anomalies in the two glacier regions show distinct alternations: in years with positive anomalies in ACN (reduced melt), SVA experiences negative anomalies (enhanced melt) and vice versa, pointing to a new state of atmospheric forcing conditions. Correlations become significant with values <-0.7 starting with the 15 yr period 1999-2013 (significance threshold r = ± 0.55 at the 95% confidence level, accounting for autocorrelation 27 ).
To analyse the relationship with atmospheric conditions, we categorize the yearly mass-balance anomalies of ACN and SVA into four types on the basis of their sign and synchronicity: concurrent mass-balance anomalies in both regions, either positive (+|+) or negative (-|-), and contrary mass-balance anomalies, either positive in ACN and negative in SVA ( + | -) or the other way around (-|+). We count the years belonging to each of the four categories within running 15 yr bins, approximately the time span of continuous data from the first GRACE mission (2002-2017). The MAR3.11 reconstruction reveals that bins centred between 1960 and 1973 include a majority of years with same-sign anomalies (9 or more out of 15 years; Fig. 2b). Starting in the early 2000s, opposite mass-balance anomalies become dominant in both glacier regions (for example, 10 or more out of 15 years starting 2001). Considering bins between the late 1970s and the late 1990s, asynchronous mass-balance years prevail in the first decade, while no single category dominates in the following decade (Fig. 2b).

Atmospheric drivers of mass change in the past and today
We analyse the atmospheric conditions in years with dominantly synchronous (1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977) and asynchronous (2004-2018) mass-balance behaviour, corresponding to the periods with the   highest positive and negative mass-balance correlations (Fig. 2b). We focus on the months June, July and August (JJA) as fluctuations in the summer melt rate cause the greatest year-to-year variability in the annual mass balances of both regions (75% for ACN and 68% for SVA) 11,17,18 (Extended Data Fig. 3). Consequently, we investigate atmospheric variables characterizing the summer circulation and associated melt potential in the Arctic, namely the geopotential height of the 500 hPa pressure surface (z500), the vertically integrated northward heat flux and the 2 m air temperature (t2m; Extended Data Fig. 4) from NCEP/NCAR, ERA5 and ERA20C monthly reanalysis data (Methods).
The correlation of z500 and northward heat flux with the annual mass-balance anomalies for both regions and time windows 1963-1977 and 2004-2018 are shown in Fig. 3. The early period, characterized by synchronous mass-balance behaviour, shows correlations with the largest values over the respective glacier region, suggesting a regional atmospheric forcing of the mass balance (Fig. 3a). Negative correlations indicate that reduced z500s in JJA lead to higher mass balances (less melting) and vice versa. However, mass balances also correlate significantly with z500 in other areas of the Arctic (P < 0.05). Between 2004 and 2018, by contrast, the correlations for ACN and SVA reveal nearly complementary patterns for both regions (Fig. 3b). Interestingly and dissimilar to the early period, the (anti-) correlation with the mass balance does not peak over the glacier regions themselves. Instead, extensive areas of statistically significant, alternating correlation flank the regions, leaving areas with coefficients around zero west of ACN and over SVA. Employing GRACE/GRACE-FO mass-balance anomalies instead of SMB anomalies for 2004-2018 produces very consistent correlation patterns with z500 and northward heat flux (Extended Data Fig. 5). This shift in correlation patterns points to an increased importance of longitudinal pressure gradients and heat transport for the mass balances of ACN and SVA.
To further analyse this, we utilize the JJA vertically integrated northward heat flux as an indicator of the additional energy available for enhanced melt as a result of the heat advection across latitudes. The correlation patterns confirm two distinct channels of heat transport to govern the recent mass-balance years (Fig. 3d), not pronounced in the early years for ACN and to a lesser extent for SVA (Fig. 3c). Reduced northward heat flux slightly south of ACN prevails during times of positive mass-balance anomalies (reduced melt conditions), while SVA experiences enhanced northward heat flux (enhanced melt conditions), and the other way around. This alternating relationship with its spatial similarity has manifested itself only in recent years.
Next, we investigate the effect of uncertainties in the SMB reconstruction (1950-2019) on the categorization mass-balance years. For this, we randomize the categorization on the basis of random sampling (N = 5,000) of the annual SMB anomalies in the presence of SMB uncertainties (Methods). From the resulting ensemble of annual categorizations, we infer the probability of any 15 yr bin containing a majority of either synchronous or asynchronous mass-balance years, with the latter shown in Fig. 4 (symmetry for synchronous years inferred). The SMB reconstruction suggests that the synchronicity underwent decadal variations, swinging back and forth between synchronous and asynchronous states. Before 2000, there is no clear evidence of a dominant synchronicity as at least 10% of the random samples fall into the respective opposite category. Starting in the early 2000s, asynchronous behaviour has become very likely (following the IPCC definitions 28 ), with probabilities exceeding the significance threshold of 5% for the 15 yr bin centred on 2004 (Fig. 4).
This underpins the strong anti-correlations shown in Fig. 2b, promoted by an increasing magnitude of mass-balance variability seen in the recent two decades (Figs. 1c and 2a). The probability distributions of the number of years for each individual category are shown in Extended Data Fig. 6. The distinctly higher probability of asynchronous years dominating the recent 15 yr bins, compared with their likelihood in the years before 2000, is further depicted in Extended Data Fig. 7.

Changing summer circulation in the Arctic
On the basis of the four categories of mass-balance years shown in Fig. 2b, together with their likelihood shown in Fig. 4, we examine the corresponding atmospheric circulation in summer (JJA), focusing on the geopotential height of 500 hPa and the wind velocity at 250 hPa, reflecting the polar front jet-stream level. Figure 5a,b shows two composites, (1) the JJA z500 mean of positive (+|+) minus the JJA mean of negative (-|-) years and (2) the JJA z500 mean of opposite-sign mass-balance years, positive in ACN (+|-) minus positive in SVA (-|+). The same-sign composite difference (Fig. 5a) exhibits a large circumpolar area of significant negative z500 values over the Arctic Ocean, amplifying the zonal pattern of pressure differences and mid-tropospheric winds (Fig. 5a). By contrast, the composite difference of opposite-sign years depicts a wave pattern of z500 anomalies across the Arctic domain (Fig. 5b). The difference between (+|-) and (-|+) years, therefore, corresponds to intensified meridional flow and more northward heat flux over the Arctic North Atlantic and towards SVA with its inherent impact on the near-surface air temperature and related glacier melt. For same-sign mass-balance anomalies, the wind composites at 250 hPa (jet-stream level) show dominantly circumpolar flow of air masses affecting both regions in years with same-sign mass balances (Fig.  5c). A consistent zonal westerly flow over Northern Europe indicates a canonical maritime summer climate state, similar to a negative summer North Atlantic Oscillation state. The perturbation in the other years differs greatly between ACN and SVA (Fig. 5d), showing a pronounced planetary wave structure passing south of Greenland. During years with more positive mass balances in ACN, air masses tend to move southwards over ACN and northwards over SVA. The latter circulation patterns occur more frequently in recent  28 . Bins with less than 5% synchronous years prevail since 2004 (highlighted in yellow). As one bin consists of 15 consecutive years, a minimum of 8 (a)synchronous mass-balance years determine the category. Since all years are sorted into one of the two categories, the probability of synchronous mass-balance behaviour is P s = 100% - P a (not shown). years (Fig. 5e). A high-pressure anomaly is located above Northern Europe, deflecting maritime air masses and promoting stationary weather patterns. In addition, the increase in frequency of the (-|+) pattern shown in Fig. 2b and Extended Data Fig. 6 is in line with enhanced blocking conditions over Greenland. Individual composites for each category underlying the composite differences are shown in Extended Data Figs. 8 and 9. The uncertainties in the z500 composites arising from uncertainties in the SMB reconstruction are presented in Extended Data Fig. 10.

Conclusions
We show that years with predominantly meridional circulation have become more frequent in the past two decades, superseding zonal circulation patterns since the late 1990s (Figs. 4 and 5e). The occurrence of years associated with zonal circulation drops from roughly 73% in the time frame 1963-1977 to only 20% during 2004-2018 (Supplementary Table 1). Today, years with primarily meridional circulation occur about three to four times more often than years with zonal circulation, pointing towards the fact that melt and thus glacier mass balances in the North Atlantic and northeast Canadian Arctic are increasingly controlled by meridional airflow during summer.
The interannual and inter-regional temperature variability over the investigated Arctic glaciers has increased with the meridionalization of atmospheric flow and is superimposed on the strong thermodynamic warming of the Arctic. We have shown that the intensification of meridional heat advection from 1950 to 2019 has directly impacted the glacier mass balances of ACN and SVA, promoting more frequent exceptional mass loss and greater changes from one year to the next. For example, in 2011 and 2012, GRACE measured substantial losses in ACN, accompanied by weak losses in SVA, while recording the reverse for the following year, 2013 (Fig. 1c).
We regard this asynchrony in the mass-balance behaviour to be a characteristic of the 'new Arctic' 3 , in line with a statistically significant increase of quasi-stationary planetary wave patterns 29 and enhanced blocking situations over Greenland 15,16 , affecting mid-latitude weather patterns over the Northern Hemisphere 2,30-35 . Other factors include trends or non-stationarity in the summer North Atlantic Oscillation 36 or trends in synoptic variance 37 . In general, a combination of natural decadal variability and anthropogenic effects is expected.
The exact impact of increasing temperature variability on Arctic systems and their potential ensuing tipping points is yet unknown,  as are the detailed mechanisms of the polar to mid-latitude linkages.
To date, global models fail to reproduce the observed atmospheric characteristics of the new Arctic 31,38 . Therefore, projections may grossly underestimate the climate impacts on the Arctic cryosphere, including 39 ecosystems and indigenous communities [40][41][42] .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41558-021-01275-4.  43 , we combine the SDS datasets by weighting each monthly solution coefficient-wise with the inverse of the uncertainty associated with each GRACE/GRACE-FO coefficient. For this, we first estimate the overall noise level for each solution on the basis of the degree power in the noise-dominated spectral range of spherical harmonic degrees 40 to 60, after subtracting bias, trend, annual, and semi-annual and temporal variations longer than four months (using a moving-average filter) from the GRACE/GRACE-FO coefficients' time series. The resulting empirical estimate of the noise level is used to calibrate the GRACE/ GRACE-FO coefficient uncertainties, provided by the SDS centres along with the Level 2 data, which are then used to determine the optimal weight of each solution for each GRACE/GRACE-FO coefficient. The combination is done for the detrended coefficient time series; the trend is later restored by applying equal weights to trends obtained from all SDS centres. We carry out the combination on detrended data as the differences in the trends appear to be systematic and arising from different processing choices of the SDS centres (Extended Data Fig. 1). Therefore, the monthly weights derived from coefficients beyond degree and order 40 are not representative for the relative uncertainty of the trends. Applying these weights to the different temporal linear signals would introduce artificial monthly temporal variability. Therefore, we remove the trends before combination and restore them, assuming equal weights for each SDS centre, in the combined solution. Extended Data Fig.  1 shows time series of GRACE/GRACE-FO mass changes and their uncertainties, respectively, for individual SDS solutions, as well as their combination.

GRACE/GRACE-FO corrections.
We use Stokes potential coefficients up to degree and order 60 and apply the following common corrections to the Level 2 data: (1) insertion of degree-1 coefficients-which are not recovered by GRACE/ GRACE-FO-provided by the SDS on the basis of ref. 44 , an improvement to the original method proposed by ref. 45 (for each SDS dataset, these coefficients are available as GRACE Technical Note 13 46 ) and (2) replacement of highly uncertain C 20 coefficients from GRACE/GRACE-FO by more accurate estimates from satellite-laser ranging 47 (GRACE Technical Note 11 48 ). Following the recommendation of ref. 49 adopted by the GRACE/GRACE-FO SDS centres, we replace the C 30 coefficients from GRACE/GRACE-FO with the satellite-laser-ranging estimate provided in the GRACE Technical Note 14 50 , starting August 2016. We apply the ICE6G Model 51 to correct for the linear trends in the gravity field caused by glacial isostatic adjustment. The linear trend caused by hydrological variations is subtracted on the basis of the average of the two models presented in ref. 4 . No correction is applied for glacial isostatic adjustment caused by the retreat after the Little Ice Age as it approximates zero for both regions 4 .

GRACE/GRACE-FO gravimetric inversion.
We estimate mass changes from the GRACE/GRACE-FO coefficients of degree j and order m representing the gravity field, C GRACE jm (ti) at the time t i of the ith measurement epoch using forward-modelling-based inversion approach in the spatial domain, an adaption of the approach previously applied to Greenland 52,53 and Antarctica 54 . The regions of ACN and SVA and other glaciated areas, n = 1…N, are geographically defined using a global subdivision of glacier regions 55,56 and the Randolph Glacier Inventory Version 4.0 57 . First, we model the spatial representation of the gravity field signal arising from a known uniform mass distribution of 1 Gt in magnitude within each region, C n (Ω), where Ω stands for the spherical colatitude ϑ and longitude φ, and hence Ω = (ϑ,φ). Then, the mass change is obtained by estimating an adjustment factor, m n (t i ), to the known mass change of 1 Gt, minimizing  Fig. 1).
Signal leakage out of the region of interest constitutes an additional source of uncertainty 58 . While the forward-modelling approach mitigates this leakage effect, it does not eliminate it completely. The fundamental resolution of the GRACE/ GRACE-FO data used here (maximum degree and order 60) cannot fully capture the mass-change information at the spatial scales of the Randolph Glacier Inventory Version 4.0 57 . Therefore, mass-change estimates differ due to the processing/ post-processing strategy applied. To quantify uncertainties associated with processing/post-processing choices, we analyse the differences between our estimates and those obtained by ref. 4 (updated to November 2019 by the authors; 179 monthly solutions). We estimate the uncertainty related to the processing choices in the gravimetric inversion to be ±15 Gt and ±5 Gt (1σ) for ACN and SVA, respectively.
The total median uncertainty of monthly mass changes for ACN and SVA obtained by quadratic summation of both uncertainty sources is estimated to be ±18 Gt and ±8 Gt (1σ), respectively. Propagation of these uncertainties to the annual mass-balance rates yields a mean uncertainty of ±12 Gt yr -1 and ±5 Gt yr -1 (1σ). For comparison, using the gridded Level 3 data product CSR RL06 mascons 59 (April 2002 to April 2021) yields annual rates that differ from our estimates within ±4 Gt yr -1 for ACN and ±5 Gt yr -1 for SVA.
Uncertainties of the rates of mass change for ACN and SVA for 2002-2021 (Fig. 1a,b) amount to ±4.2 Gt yr -1 and ±1.1 Gt yr -1 (1σ), respectively. The uncertainties are the quadratic sum of propagated total monthly uncertainties and the maximum difference between rates of mass change obtained here and those using the processing of ref. 4 . For this estimate, we randomly select 163 of the 179 monthly solutions following a bootstrap approach (n = 200) and estimate a maximum difference between the rates of ±1.3 (ACN) and ±0.7 (SVA), using the same selection of months. Using the gridded Level 3 CSR RL06 mascons yields a rate of mass change of -32.7 Gt yr -1 and -13.2 Gt yr -1 for ACN and SVA, respectively, well within uncertainties of the values presented here (-30.9 ± 8.3 Gt yr -1 for ACN and -13.2 ± 2.3 Gt yr -1 for SVA).
MAR3.11 SMB. The SMB and its components are obtained from the output of the regional climate model MAR3.11 8 . MAR3.11 output is obtained at 15 km horizontal resolution for the period 1950-2019, six-hourly forced at its lateral boundaries with ERA5 reanalysis 13 . Our analysis shows that alternative MAR3.11 simulations, using National Centers for Environmental Prediction (NCEP)/ National Center for Atmospheric Research (NCAR) forcing at a spatial resolution of 20 km, provide similar results, albeit with a slightly larger difference to the GRACE/GRACE-FO data, making the ERA5 forced time series the preferential model output. Simulations based on the regional climate model RACMO2p1 21 for 1958-2015 are consistent with MAR3.11 output using ERA5, in line with a high level of agreement between MAR3.9/MAR3.10 and RACMO2p1 over the Greenland Ice Sheet 12,52 .
Grid cells of non-or partially glaciated land surfaces in ACN and SVA are excluded using the land surface mask of MAR3.11. From the monthly fluxes, cumulative storage changes are calculated by summation over time. The MAR3.11 monthly time series is interpolated onto the GRACE/GRACE-FO time epochs, t i , when both data are directly compared. To account for ice-dynamic discharge of marine-terminating glaciers, linear trends of GRACE/GRACE-FO and MAR3.11 SMB from 2002 to 2019 are estimated and their difference obtained as discharge D*, reconciling trends in both datasets for the common observation period. For SVA, a mass loss of 4.4 ± 0.8 Gt related to the episodic surge event of Austfonna Glacier in 2013 25 is subtracted from the SMB output before estimating D*. To complete the mass balance, the same value of D* is assumed in the MAR3.11 time series starting 1950, neglecting a possible gradual increase in ice-dynamic loss 26 . However, as we remove a polynomial fit of order n = 3 from the mass balances to focus on interannual changes, any long-term changes of D* are removed, thus not altering the results of the study.
Uncertainty of SMB estimate. We empirically estimate the uncertainty of the cumulative monthly SMB on the basis of the comparison with GRACE/ GRACE-FO data. For this, we assume that the mean monthly difference between the reconciled GRACE/GRACE-FO and SMB time series arises from the uncertainties both in GRACE/GRACE-FO and in the SMB data, σ 2 Total = σ 2 G + σ 2 SMB . The GRACE/GRACE-FO uncertainty, σ 2 G , is estimated on the basis of the detrended difference between two independent time series for each glacier system, one of which is derived in this paper while the other one is an update to ref. 4 . Therefore, the GRACE/GRACE-FO uncertainty estimate includes stochastic noise induced by measurement errors as well as difference caused by processing choices made during the inversion of the gravity field changes for mass changes. The resulting monthly SMB uncertainties are ±12 Gt for ACN and ±9 Gt for SVA (±1σ), leading to an uncertainty of the annual balance anomalies of ±7 Gt yr -1 and ±5 Gt yr -1 , respectively. These uncertainties are comparable to the independent estimates of ±12 Gt yr -1 for ACN 21 and ±7 Gt yr -1 for SVA 8 , which, however, represent full SMB temporal components in the annual mass balance, not mass-balance anomalies.
NCEP/NCAR, ERA5 and ERA20C data. We analyse NCEP/NCAR and ERA5 monthly reanalysis atmospheric data, with a spatial resolution of 2.5° × 2.5° and 0.25° × 0.25°, respectively, to study the climatic drivers of mass balance in ACN and SVA. In particular, we investigate the following variables: geopotential height at 500 hPa (z500; NCEP), vertically integrated northward heat flux (ERA5/ERA20C), wind velocity at 250 hPa (NCEP) and t2m (NCEP). As ERA5 monthly reanalysis datasets reach back only to 1979 (at the time of this analysis), monthly data for northward heat flux is extended to 1950 using ERA20C data to obtain the full time series from 1950 to 2020. For consistency, ERA20C and ERA5 output is resampled to the same equiangular grid of 1° × 1° using geographic bi-linear interpolation. Correlations with the mass-balance anomalies as well as the composites are calculated for each grid point using the JJA mean of each selected variable for the years in question. To assess the consistency between datasets, we calculate the spatial correlations between the annual mass balances of ACN/SVA and the northward heat flux on the basis of ERA20C and ERA5 within the overlapping period 1979-1995. The resulting standard deviation of the difference of the correlation coefficients is <0.05 in both cases, indicating acceptable consistency between datasets. The results are plotted with the North-Polar equal-area azimuthal projection centred on the longitude of 0° (EPSG: 3897).

Correlation of atmospheric variables and mass balances.
We calculate the Pearson correlation coefficient (lag zero) of the summer (JJA) mean of the atmospheric variables z500 and northward heat flux (Fig. 3) and t2m (Extended Data Fig. 4) for the periods 1963-1977 and 2004-2018 with the respective annual mass balance in ACN and SVA. Significance of the correlation is determined using a critical value of 5% for N -2 degrees of freedom (P < 0.05) based on the Student's t test and marked in the plot. We account for autocorrelation in SMB time series by reducing the degrees of freedom following ref. 27 . Categorization of mass-balance years (Fig. 2) is based on the anomalies of the annual mass balance after polynomial detrending of the MAR3.11 SMB time series. Temporal correlations between the regions' annual mass-balance anomalies are represented by Pearson's correlation coefficients of lag zero in moving windows of N = 15 years, where significance is determined using Student's t test for P < 0.05, considering autocorrelation of each region's SMB time series in the determination of the degrees of freedom 27 . Significance thresholds for the composite differences shown in Fig. 5 are determined on the basis of the Wilcoxon test (P < 0.05).

Data availability
Monthly datasets of summer (June-July-August, JJA) two-metre surface air temperature (t2m) and 500 hPa geopotential heights (z500) were obtained from the latest reanalysis of the National Centers for Environmental Prediction and National Center for Atmospheric Research (NCEP/NCAR Reanalysis 1, 1948-2020 60 ). NCEP reanalysis data are provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their website at https://www.esrl.noaa.gov/psd/. Monthly data of JJA vertically integrated northward heat flux were obtained from the latest reanalysis of the European Centre for Medium-Range Weather Forecasts (ERA5) for 1979-2019 13 , augmented with ERA20C for 1950-1978 61 (https://cds. climate.copernicus.eu/). GRACE/GRACE-FO gravity fields (Level 2 data) and supporting documentation are freely available from the websites of the Science Data Systems Centres and may be accessed at http://podaac.jpl.nasa.gov and http://isdc.gfz-potsdam.de. The gridded Level 3 data product CSR RL06 GRACE/ GRACE-FO mascon solutions can be downloaded from http://www2.csr.utexas. edu/grace. MAR3.11 output of the surface mass balance and its components are freely available from the authors upon request and without conditions. To submit a request, please contact Xavier Fettweis: xavier.fettweis@uliege.be. Source data for the reproduction of Figs. 1, 2 and 4 are available for download 62 , allowing the reproduction of Figs. 3 and 5 using public archives of NCEP (https://www.esrl. noaa.gov/psd/) and ERA (https://cds.climate.copernicus.eu/) reanalysis data. Additional mass time series for ACN and SVA from GRACE/GRACE-FO data are freely available without further conditions from the corresponding author.

Code availability
In the GRACE/GRACE-FO data processing, spherical harmonic functions were handled using the software package of Frederik J. Simons, kindly provided at http:// geoweb.princeton.edu/people/simons/software.html (last accessed 17 June 2021). Geographic maps were produced using the Matlab Mapping Toolbox based on publicly available digital datasets.  (1963-1977; purple bars) and anti-correlation (2004-2018; brown bars) between the mass balances in ACN and SVA. The mean and standard deviation of the number of years belonging to each time period and category is specified with matching colours in the diagram, together with the amount of overlap (%; darker shading of bars) between the early and late period. Without noise in the SMB estimates, each period would be associated with a distinct number of years in each category. With noise dominating the SMB data, the same number of years would belong to each category, impeding a meaningful categorization. Fig. 7 | ensemble distributions of categorization. a, Overlap of the probability distributions of random realisations (N=5000) with respect to their categorization splits (number of synchronous vs. asynchronous years) varying from just under 40% to well over 90% until the early 2000s. The overlap is expressed as fraction and number of members. An example of the distributions of the four categories for two specific periods is shown in Extended Data Fig. 6. b, Fraction and number of random realisations exceeding the number of asynchronous years in 2011. Similar to Fig. 4 the graph illustrates that phases of asynchronous behaviour have existed in the past, however, the portion of asynchronous members exceeding the value of 2011 (11.3 yr) remains below 5%.