Intense atmospheric rivers can weaken ice shelf stability at the Antarctic Peninsula

The disintegration of the ice shelves along the Antarctic Peninsula have spurred much discussion on the various processes leading to their eventual dramatic collapse, but without a consensus on an atmospheric forcing that could connect these processes. Here, using an atmospheric river detection algorithm along with a regional climate model and satellite observations, we show that the most intense atmospheric rivers induce extremes in temperature, surface melt, sea-ice disintegration, or large swells that destabilize the ice shelves with 40% probability. This was observed during the collapses of the Larsen A and B ice shelves during the summers of 1995 and 2002 respectively. Overall, 60% of calving events from 2000–2020 were triggered by atmospheric rivers. The loss of the buttressing effect from these ice shelves leads to further continental ice loss and subsequent sea-level rise. Under future warming projections, the Larsen C ice shelf will be at-risk from the same processes. The most intense atmospheric rivers to hit the Antarctic Peninsula induce extremes in temperature, surface melt, sea ice disintegration or swell that destabilize the ice shelves with 40% probability, suggest analyses of observations and regional climate model simulations.

I n the past 30 years, the large and dramatic collapses of the Larsen A, Larsen B, along with other major ice shelves along the Antarctic Peninsula (AP) raised fears for the fate of other ice shelves controlling the outgoing continental ice. The associated loss of the buttressing effect from ice shelves and consequent acceleration of continental ice loss makes a critical contribution to sea-level rise 1 , while global warming could push the Antarctic ice shelves beyond a stability threshold during the twenty-first century [2][3][4] .
Among the many preconditioning processes of ice-shelf collapses, that is any event that pushes the ice shelf closer to collapse, atmospheric forcings are seemingly considered as disparate sources of ice-shelf weakening compared to non-atmospheric related processes [5][6][7][8][9] . Intense surface melting, generally associated with downsloping foehn winds descending over the leeward Larsen ice shelves 10,11 , can drive firn air depletion allowing crevasses to later be filled with water generating a hydrostatic pressure in a process known as hydrofracturing 12 . Speed limits on hydrofracture-induced ice-shelf disintegration postulates that collapse occurs when a large surface area of the ice shelf experiences surface melt and melt pond formation as observed during the collapse of the Larsen B in 2002 13 . It is also observed that the ice shelf's final collapse can be triggered by storm generating swells causing calving on the ice-shelf front in the absence of a regional sea-ice cover 14 as observed in January 1995 during the sudden collapse of the Larsen A 14,15 and again with the Larsen B collapse 14 . The importance of these atmospheric related processes on ice-shelf collapse are understood on an individual basis. However, in question is the ice-shelf weakening potential if an atmospheric forcing triggered compound extreme events. Regarding potential ice-shelf destabilization risk, what is the likelihood that these extreme events co-occur with a singular atmospheric event sensitive to large-scale atmospheric circulation variability and does the risk increase when these extreme events are compounded together? Here, we show that atmospheric rivers (ARs) aggregate seemingly disparate extreme events and are linked to ice-shelf collapse via melt-induced hydrofracturing and swellinduced strain from sea-ice clearing.
ARs are narrow long bands of enhanced moisture fluxes originating from the mid-latitudes and sub-tropics 16 . They have been associated with particular instances of high temperature, moisture, and wind speed in the lower troposphere [17][18][19] along with sometimes co-occurring with foehn winds while causing major melt events 20,21 and being linked with instances of sea-ice decay 22 . However, a systematic long-term analysis of extreme events co-occurring with ARs still needs to be conducted along the AP to show their role in ice-shelf weakening and the eventual initiation of ice-shelf calving and disintegration events.
Here, we used an AR detection algorithm previously utilized for studying Antarctic AR climatology, Antarctic ice sheet precipitation, and melt events in West Antarctica 23 , a polar-specialized regional climate model, Modèle Atmosphérique Régional 24 (MAR), driven by ERA-5, and a series of satellite observations (see Methods) to retrieve calving events, surface melt occurrences, swell heights, and sea ice extent around the AP. We analyze the link between ARs and the Larsen A and B ice-shelf calving and collapses and the major influence of AR on the processes behind ice-shelf weakening (i.e., extreme temperature, surface melt, runoff, melt pond formation, sea-ice removal, ocean swell and sea-level height-induced ice-shelf flexure, foehn wind, and iceberg calving). Particularly we address the probability of intense ARs generating compound extreme events linked to ice-shelf collapse triggers.

Results and discussion
Ice-shelf calving and collapse co-occurrences with ARs. To examine co-occurrences between AR landfalls and iceberg calving along the Larsen ice shelves, we produced a climatology of AR landfalls along the AP from 1989 to 2020 using output from the AR detection algorithm with input from Modern-Era Retrospective analysis for Research and Applications, Version 2 25 (MERRA-2) reanalysis where instantaneous intensity is measured by integrated vapor transport (IVT) (see Fig. 1 and Methods for AR frequency and specifics on the AR detection algorithm, respectively). We also examined calving and collapse events along Larsen A, B and C (between 64.5°S and 68.5°S) using Moderate Resolution Imaging Spectroradiometer (MODIS)-Terra and Aqua visible imagery from August-March 2000-2020. We observe that AR landfalls along the AP are slightly more frequent than other coastal regions of the continent 23 (around 1-5 landfalls per austral summer, i.e., between December-March) and occur during strong downstream blocking conditions ( Supplementary Fig. 1) and the positive phase of the Southern Annular Mode thus agreeing with previous Antarctic AR climatology assessments 20,23,26 .
We found 21 major calving and collapse events between 2000-2020, all occurring before a decrease in AR activity after 2010 (Supplementary Table 1). Thirteen out of the 21 calving events were preceded by an AR landfall within 5 days prior (Table 1; the Larsen C calving in July 2017 is not considered here as it could not be retrieved during austral winter using MODIS although an AR was detected within 5 days prior of the supposed calving date). A statistical analysis based on thousands of artificial AR occurrences at randomly selected dates between August and March from 2000 to 2020 confirms that ARs and these calving/ collapse events were not independent at p = 5.3 × 10 −5 . The number of co-occurrences decrease if we consider AR landfalls within <5 days prior to calving events, but the co-occurrences are still statistically significant (except if we only consider cooccurrence on the same date) confirming that these ARs and calving/collapse events were not independent. In contrast with high-swell events (swells within the 85th percentile of height), ARs are more likely to cause a calving/collapse event with only 9 out of 21 events co-occurring with 5 days prior of a high swell.
As an example of AR impacts on ice-shelf weakening processes, an intense AR that made landfall from January 24-26, 2008, is particularly striking. MODIS satellite imagery shows that this historically intense AR detection (IVT~962 kg m −1 s −1 , third highest intensity of all AR landfalls in our AR climatology) disintegrated and fragmented nearly all the land-fast ice in the Larsen A and B embayments and generated 6.3 Gt (7.2 Gt) of runoff (meltwater) simulated by MAR likely leading to a calving event 6 days later (Fig. 2). This calving event (not included in our 15 co-occurring calving events as it appears to have occurred 1 day after our self-imposed 5-day selection window) led to a 11 km retreat of the eastern region of the remnant Larsen B (named Scar Inlet ice shelf) 27 . Using Flexible Particle Dispersion Model back-trajectories, we observed air parcels associated with this AR event traversing a region of sea surface temperatures in the southeastern Pacific 2-3°C above average before reaching the AP 2 days later (Fig. 2e). Over these anomalously warm waters, the temperature and specific humidity content of the air parcels were boosted before decreasing as the parcels were isentropically lifted upon being advected over cooler sea surface temperatures ( Supplementary Fig. 2). This process was observed in other important AR events on the AP during the study period (i.e., Larsen A and B collapses, 2005 Larsen C calving event, and recently during the breakup of land-fast ice in the Larsen B embayment in January 2022).
Other striking examples are the multiple AR events that preceded the collapse of the Larsen B in austral summer 2002 and the historically intense AR preceding the collapse of the Larsen A ice shelf and associated iceberg calving (~1700 km 2 ) from the Larsen B in 1995 15,28 (Supplementary Fig. 3). In all these cases, the AR generated several Gt of meltwater and runoff, large reductions of sea-ice concentrations in the Larsen embayment, strong winds, and large swells (Table 1; see Supplementary Note 1 for more examples of AR landfalls and their associated impacts on processes connected to ice-shelf weakening).
Atmospheric river impacts on processes detrimental to iceshelf stability. Here focusing on processes related to ice-shelf destabilization, we demonstrate that ARs generate co-occurring extreme temperature, melt, runoff, and swell events, by analyzing values within the 99.9th percentile that occurred 24 h before and after an AR landfall. Given AR landfall frequency on the AP, a result above 61% of co-occurrences is significant at the level of 1%. We also analyzed the rainfall percentage and composite seaice fraction change associated with ARs. Our focus is on the summer months (December-March), when the most intense ARs  i.e., maximum IVT above 900 kg m −1 s −1 ) and AR-related surface melt events have historically occurred ( Supplementary  Fig. 4).
First, ARs are observed to strongly affect extreme hightemperature events. Figure 3a shows that nearly all occurrences of 2 m air temperatures above the 99.9th percentile simulated by MAR (61-100% out of~40 at a 3-h timestep) and of most temperatures above the 99th percentile (61-70% out of~280 occurrences at a 3-h timestep) occur within 24 h before and after an AR landfall ( Supplementary Fig. 5a, b). In fact, the previous two Antarctic temperature records recorded at Esperanza Base in March 2015 and in February 2020 coincided with intense AR events that triggered foehn wind events 21,29 (see Supplementary  Fig. 6). Higher nighttime temperatures from enhanced downward longwave radiation were also observed meaning that the melted water during ARs is less likely to refreeze in the firn and thus contribute to runoff 30 .
Second, we observe (>61 %) of melt rates at or above the 99.9th percentile along portions of the Larsen and Wilkins ice shelves and higher percentages on the glaciers in the mountainous regions occurred within the 24-h period surrounding AR landfalls ( Fig. 3b and Supplementary Fig. 5c, d). To complement our analysis of surface melt and the overall presence of liquid water on the snow surface based on MAR simulated results, we used passive microwave radiometer data from Scanning Multichannel Microwave Radiometer (SMMR) and Special Sensor Microwave/ Imager (SSM/I). Results show that ARs typically engender neutral to largely positive melt extent anomalies in daily climatological melt extents along the Larsen C ice shelf especially outside of the peak firn saturation months of February and March (Supplementary Fig. 7). During these months, the presence of a saturated firn can obscure the occurrences of new melt detected by passive microwave remote sensing 31 . Still, from October to March, we observe that two thirds of AR landfalls were associated with positive melt anomalies, suggesting that ARs generally produce anomalously widespread melt.
Third, related with meltwater, rainfall events lead to periods of extremely high rates of firn saturation and runoff 32 . Runoff here refers to liquid water beyond the saturation point of the snowpack that does not necessarily drain into the ocean and could remain on the ice shelf. From 40 to 60% of total annual rainfall along the Larsen and Wilkins ice shelves are linked to ARs (Fig. 3d) whereas ARs days also produced a statistically significant higher amount of rainfall than non-AR days (p-value < 0.0025). Runoff rates thus show similar relationships to AR activity as melt rates, the only difference being that runoff is predominately constrained to the lower-elevation ice shelves ( Fig. 3c and Supplementary  Fig. 5e, f). High runoff values, produced by melting and rainfall when the firn is already saturated with liquid water, are detrimental for ice-shelf stability 33 . Indeed, these high melt and runoff rates lead to melt pond formation as demonstrated by the significant correlation between melt pond observations from MODIS and summer AR occurrences 34 (r = 0.78; Supplementary  Fig. 8).
Fourth, beyond processes related with melt and hydrofracturing, ARs are also observed to induce a rapid clearing of sea ice along the ice shelf fronts allowing co-occurring anomalously large swells to reach the ice-shelf and apply strain 14 . This swell-induced stress leads to the calving of large icebergs and could result in the retreat of the compressive arch like what occurred on the Larsen B in 2002 14,35 . On average, the strong ARs (those with integrated vapor transport (IVT) > 400 kg m −1 s −1 upon landfall) cause a 2-day sea-ice fraction decrease of nearly 10% (i.e., about 5 sigma of sea-ice variation in this region) on the northern edge of the Larsen A and B embayment (Fig. 3e). This average sea-ice fraction decrease is slightly larger when considering 4-day sea-ice change. Nearly all large and rapid sea-ice disintegration events of around a 10% decrease over 2 days east of the Larsen ice shelves (in the region labeled EAP) co-occurred with AR landfalls ( Supplementary Fig. 9). The average 2-day sea-ice change during intense AR days in this region (2.5% reduction) is statistically significantly higher than non-AR days (0.1% increase; p-value < 0.0025). At the same time, the winds found near the surface during an AR within the low-level jet are often extreme as noticed during the AR landfall coincident with the collapse of the Larsen A in 1995 ( Supplementary Fig. 10a) and can lead to large swell heights. When sea ice is depleted east of the Larsen ice shelves, 61-80% of swells in the 99.9th percentile of height co-occurred with an AR landfall (Fig. 3f). Consequently, AR landfalls prompt a state of high stress on the ice shelves through the combination of surface melt induced melt pond formation, leading to hydrofracturing and wave strain along the ice-shelf front. The tendency of ARs to create widespread melt pond formation makes them possible precursors of ice-shelf collapse via hydrofracturing cascades like what was observed on the Larsen B in 2002 13 . This combined with wind stress and radiative forcing leading to sea-ice clearing thus allowing swells to apply strain along the ice-shelf fronts 14,22 makes ARs a unique forcing of ice-shelf weakening. These conditions occur during periods of enhanced poleward heat advection from anticyclonic activity just east of the northern AP 36 ( Supplementary Fig. 1). Figure 4 provides a visualization summarizing all the processes during AR landfalls that are linked to ice-shelf weakening. The AR example in the schematic showcases an AR making landfall perpendicular to the ice shelf and generating a foehn wind along with the associated sensible heat, downward longwave, and downward shortwave radiative fluxes as the airflow adiabatically warms and dries over the leeward ice shelves while causing sea-ice disintegration and swells along the ice-shelf front.
We studied other possible connections between ARs and iceshelf weakening and did not find conclusive evidence for a direct relationship between AR activity and ice-shelf basal melting, one of the long-term precursors of ice-shelf collapse 37 . Indeed, looking at ocean hindcasts from 1979 to 2018 38 , we did not find any significant correlation between either yearly AR frequency or yearly cumulative maximum IVT and either the ocean temperature at the depth of the ice-shelf drafts or Ekman pumping (a driver of basal ice-shelf melt variations 39   shelf may have contributed to trigger a calving event 40 . However, looking at sea surface height anomalies in the global ocean simulations from 1979 to 2015 38 , we did not identify any noteworthy sea surface height anomalies associated with ARs on both sides of the AP.
Characterizing cumulative impacts of atmospheric river events.
The previously described ice-shelf weakening processes associated with AR landfalls become more impactful and likely as the AR intensity increases as a function of maximum IVT upon landfall. In fact, these most intense ARs often have the greatest surface melt potential. The relationship between maximum IVT upon landfall and both total melt and runoff along the Larsen ice shelves is exponential in shape with the average 3-hourly melt and runoff increasing from <0.05 Gt 3 h −1 for the weakest ARs tõ 0.4 Gt 3 h −1 for the strongest (Fig. 5; melt and runoff is calculated in the areas highlighted in Fig. 5c). These most extreme ARs (IVT > 800 kg m −1 s −1 ) only occur around once every five years but have the capability to alter the cryosphere dramatically and quickly like observed in the January 2008 example (Fig. 2), the collapse of the Larsen A in 1995 ( Supplementary Fig. 3), and in the March 2015 melt event/temperature record 21 . On a subseasonal scale, the positive relationship between maximum IVT and surface melt exists every month, but the melt quantities generated by ARs are much less in November and March compared to December-February ( Supplementary Fig. 11). Delving further into the AR characteristics, melting increases most uniformly during AR landfalls with a northwesterly orientation upon landfall ( Supplementary Fig. 12). This allows for moisture transport to be perpendicular to the AP mountains and efficiently generate a leeward foehn wind. Partitioning AR IVT into wind speed and integrated water vapor (IWV) reveals that surface melt responds much more clearly to increasing IWV compared with wind speed. In fact, the relationships between max IWV/melt and max IVT/melt are very similar while the magnitude of surface melt remains basically unchanged as 900 hPa wind speeds increase ( Supplementary Fig. 13). Knowing that the highest melt rates occur with northwesterly oriented ARs, which likely implies a foehn wind is necessary for intense melt, it appears that water vapor content is a better indicator for foehn melting potential than wind speed intensity. When AR strength and impacts are characterized by considering landfall duration in addition to the maximum IVT upon landfall, like the scale used to classify ARs in western North America 41 , consistent relationships between AR intensity and various impacts emerge. As the cumulative maximum IVT of the AR increases, the more likely a temperature extreme, melting extreme, sea-ice disintegration or high-swell event occurs (Fig. 4). This is visible in Fig. 6a, where AR intensity is assessed by summing the maximum IVT values throughout the day when an AR is detected. Taken as a whole, all these events represent a physical state that promotes ice-shelf weakening as confirmed by the co-occurring calving and even collapse events with AR landfalls. For the most intense ARs at landfall, there is a nearly 20% probability of a daily temperature, melt, or runoff extreme occurring over at least half the lower-elevation portion of the northern AP (blue line and area; see Fig. 6d for domain). There is an increase in probability when the same calculation is repeated but with the added option of a substantial sea-ice decline occurring (up to 40%, red line and area), and with the option of a high-swell event east of the AP (green line and area).
These are the effects from singular AR events, however ARs often provoke temperature and moisture extremes within weather regimes that favor warmth and melting like observed in East Antarctica 42 . In fact, summer mean temperature and runoff are significantly correlated with summer cumulative AR activity (Fig. 6c) showing that ARs are more prominent during already Fig. 4 Illustration of an AR landfall on the AP. An illustration of a typical intense atmospheric river over the northern Antarctic Peninsula and the associated observed meteorological features and impacts consequential to ice-shelf stability. Also, an example of a detected AR landfall on Feb. 6th, 2020, with the corresponding IVT values. The yellow, red, and green outlines are the shape of the AR as determined by the vIVT AR detection scheme, IWV AR detection scheme, and the original Antarctic AR detection algorithm, respectively 17,19 . The satellite imagery in Supplementary Fig. 21 is an accurate representation of the AR airmass adiabatically warming and drying as it passes over the AP mountains leaving scattered cloud cover over the leeward Larsen B ice shelf.  Fig. 14). During the summer seasons of 1994/1995 and 2001/2002 when the Larsen A and B collapsed, respectively, we observed that the usually steady rate of cumulative runoff would increase in conjunction with an increase in cumulative IVT when an AR landfall occurred (Supplementary Fig. 10).
We also analyzed if iceberg calving and ice-shelf collapse events were related with AR cumulative intensity within 5 days after the AR event. When defining AR intensity using the cumulative maximum IVT throughout a discrete continuous AR landfall event, we find around 60% of calving/collapse events (13 out of 21 total events over August-March 2000-2020 plus the collapse of the Larsen A in 1995 and the calving event from the Larsen C in July 2017), occurred when an AR with at least a cumulative IVT value of 500 kg m −1 s −1 was detected within 5 days prior. Although that percentage progressively decreases as the IVT threshold increases (see Methods for details on the statistical significance of the lag between AR and calving), thus suggesting that calving requires neither a large trigger nor large swell heights once hydrofracturing and other weakening processes have already weakened the stability of the ice-shelf front. However, the major collapses of the Larsen A and B were preceded by intense ARs like the early January-February 2002 AR that caused widespread melt pond formation across the Larsen B and contributed to large swell heights (Table 1; see Supplementary Fig. 15). Even more important in terms of summer AR activity, the 2002 summer had the second highest cumulative IVT over 1980-2020, coinciding with exceptional sea-ice-free conditions around the AP and with large wave-induced flexure on the ice shelves 14 (Table 1).
Foehn wind enhancement by ARs. When ARs make landfall perpendicular to the mountains of the AP, their large moisture transport typically engenders a large latent heat release. These are ideal conditions for foehn winds, which are commonly proposed to trigger intense surface melt along the leeward AP ice shelves 20,21,29,[43][44][45] . However, the enhanced moisture transport also triggers melt through clouds containing anomalously high liquid and ice water content resulting in high downward longwave radiation 20 . Thus, it is important to analyze the relationship between these distinct phenomena, ARs and foehn winds.
To compare AR and foehn frequency and surface melt potential, we compared AR landfalls events against instances of foehn winds detected by a dedicated algorithm, which is based on relative humidity observations from six automatic weather stations (AWSs, see  Table 2). The foehn events that are related to ARs appear to have a greater duration and extent over the ice shelves than non-ARrelated foehn events. Indeed, 43% of all foehn events detected at four or more AWSs are associated with ARs (compared to 22% only detected at one station) and 60% of foehn events lasting more than 3 days are associated with ARs (compared to 23% only lasting 0-2 days). These percentages for foehn detection extent and duration increase when considering only intense ARs (IVT > 400 kg m −1 s −1 upon landfall). Further east on the Larsen C, 80% of intense AR events prompted a foehn wind at AWS2/3 (130 km from the foot of the mountains), but decreases to 40% when looking at all ARs. These intense AR/foehn events were more likely to be observed further north at AWS1 than further south at AWS6, which is similar to foehn wind patterns in general. The greater impact area during AR-related foehn events can be explained by the higher windward flow velocity 10 . ARs with their associated high wind speed are more likely to create a linear flow across the ice shelves, mechanically mixing down warm air aloft and transporting sensible heat to the surface. Weaker foehn events are typically associated with nonlinear flow over the mountains, which results in hydraulic jumps shortly after descending thus limiting the warming extent over the ice shelf 10 . This explains why most AR/foehn events were detected at the Cole Peninsula station closest to the base of the mountains. Although greater windward flow velocity and foehn impact area does not necessarily translate to increased melt magnitude (Supplementary Fig. 13). Thus, standalone foehn events generally have a smaller warming extent across the Larsen ice shelves than those associated with intense ARs. This is further observed with intense ARs (IVT > 400 kg m -1 s −1 upon landfall). When we compared all intense ARs and foehn events from 2009 to 2012 using MAR, we observed that ARs generated smaller positive 2 m temperature anomalies along the base of the leeward AP mountains, but higher anomalies across the rest of the AP (Fig. 7a), with greater melt along regions in the northern AP (Fig. 7b). While ARs generally have less sensible heat flux across the base of the AP mountains, the latent heat flux partially balances out the heat fluxes ( Fig. 7c and Supplementary Fig. 16c). The higher moisture content within ARs also causes them to have more clouds, which results in less downward shortwave radiation across the AP, but more downward longwave radiation ( Fig. 7d and Supplementary  Fig. 16a, b). Higher melt is likely explained by ARs producing greater downward longwave radiation flux anomalies (associated with higher cloud liquid/ice water contents). We also see that ARs enhance the melting potential of the foehn wind when comparing foehn events co-occurring with an AR landfall against all foehn events. Supplementary Fig. 17).
Are other Antarctic ice shelves in trouble? Although ARs are an ever-present feature of the AP climate system that cause shortterm conditions, which are consequential to ice-shelf stability, their effects have mostly been noticed during the relatively recent period of warmth in the late twentieth century 46 . Temperature extremes along the AP also reached a maximum during this time 47 , but then decreased throughout the twenty-first century, which followed a decreasing trend in AR activity (Figs. 1 and 6c). As AR activity follows average temperature trends, ARs likely only become detrimental above a certain average temperature threshold.
After a record melt year along the AP ice shelves in 2019/2020 that brought an end to the decreasing melt trend that started in 1999/2000 31,36 , it appears the remaining AP ice shelves might be entering a new phase of vulnerability from ARs. Attention is now on the Larsen C ice shelf after a giant iceberg calving event on July [10][11][12]2017. This event has been attributed to basal melting from subsurface ocean warming 37,48,49 . Curiously, a relatively strong winter AR on July 5-6, 2017, may have applied a wind stress prior to the calving, but the impacts of winter ARs on ice-shelf stability are unclear as they generally produce lower melt values compared to the summer melt season. Other calving events on the Larsen C like in January and February 2005 are more directly related to ARs ( Supplementary Fig. 18). If the Larsen C recedes beyond its compressive arch or the passive shelf ice decays 35,50 , the shelf would possibly become vulnerable to collapse from further surface melting and hydrofracturing 51 . In fact, the Larsen C is even more sensitive to AR activity than the Larsen A and B as the highest 48-h 2 m temperature changes following AR landfalls occur along the central and southern part of the Larsen C likely due to the regional mountain topography 10,21 ( Supplementary  Fig. 19). Although melt is currently rare on these parts of the Larsen C, with predictions of an extended melt season 52 , a southward migration of surface and subsurface melting conditions 39,53 , and an increase in AR frequency and intensity across the Southern Ocean in future global warming scenarios 54 , the Larsen C may pass an average temperature threshold and become increasingly vulnerable to AR behavior like the Larsen A and B ice shelves already have experienced. Thus, tracking ARs by measuring IVT and duration is a useful complementary tool in assessing the future ice-shelf disintegration risk of the Larsen C.
Increasing temperatures may put other ice shelves with similar vulnerabilities to foehn winds, hydrofracturing, and swellinduced strain at risk of collapse 2,51 . It is important to identify, which ice shelves may become threatened by AR-related extreme processes in warming scenarios 2 . Previous research cataloging which ice shelves are susceptible hydrofracturing 51 , identifying a foehn melt mechanism on the Ross Ice Shelf 55,56 , and ice shelves that experience calving from ocean swells 40,57,58 are important examples of this.

Conclusions
Here, we demonstrated that a relationship exists between AR intensity (measured by duration and maximum IVT upon landfall) and the likelihood of an extreme state of conditions occurring on the northern AP (i.e., temperature extreme, melt extreme, sea-ice disintegration, swells). Strong ARs were precursors for more than 60% of the major calving events of the Larsen A and Larsen B ice shelves since 2000 and precluded their main collapses. On the Wilkins ice shelf, the relationship with ARs and calving is less clear. ARs induce extreme temperature and melt states on the Wilkins ice shelf, but do not appear immediately related to any calving or disintegration events. The prevalence of perennial firn aquifers on the Wilkins ice shelf may explain the disparity between AR occurrences and calving/collapse events since it can delay or reduce runoff 12,59,60 . Also, this disparity could be partly justified by the topography around the Wilkins ice shelf, which is less conducive for AR-triggered foehn events compared to the leeward ice shelves 10,11,45 .
ARs are linked to only a small percentage of the total summer melt on the AP 20 , but critically cause particularly widespread and intense foehn melt events like the event that triggered the collapse of the Larsen B via a hydrofracturing cascade 13,61 . This along with the connection to sea-ice clearance and swell-induced stress also linked to the Larsen B and A collapses 14 make ARs a destabilizing force for ice shelves. The AR detection algorithm used to attribute extreme events has a high moisture transport threshold for detection compared to other global AR detection algorithms 23 , while using multiple algorithms could quantify the range of uncertainty in ice-shelf risk introduced by the choice in AR detection 16,62,63 .Nevertheless, the processes that lead to the final collapse of an ice shelf are complex, not fully understood, and likely go beyond the influence of an event that is temporally minuscule in comparison to the age of the ice shelf and thus requires further exploration. In addition, the broader scale circulation patterns and local sea surface temperature perturbations that potentially influence AR transport and surface melt over the AP need to be better understood. Indeed, the positive phase of Southern Annular Mode significantly correlated with AR activity and general warming of the AP 23,64 , while enhanced deep convection in the central tropical Pacific has been shown to promote Rossby wave propagation, leading to circulation patterns conducive for AR transport and melt over the AP 65 . Still, the broader implication here is that ice dynamic/ice sheet stability models need to account for short-term extremes in atmospheric behavior beyond changes in the mean signal like understanding how an ice-shelf responds to an AR event.

Methods
The atmospheric river detection algorithm. The AR detection algorithm used for this study is the same version used to create a climatology of ARs across Antarctica 23 . This version of the algorithm has two detection schemes based on integrated water vapor (IWV) and the v-component of the integrated vapor transport (vIVT). In both schemes, the algorithm scans 3-hourly IWV and vIVT fields from 37.5°S-80°S for values within the 98th percentile of all monthly climatological values defined from 1980 to March 2020. If contiguous grid cells in the meridional direction satisfy this condition and extend at least 20 degrees in the meridional direction, then the shape is flagged as an AR. If any grid cell in the AR shape intersects with the mask of the Antarctic Peninsula shown in Fig. 1e, then the date and the maximum IVT within the AR mask are recorded.
Here, IWV (kg m −2 ) and vIVT (kg m −1 s −1 ) are calculated based on the specific humidity, q (kg kg −1 ), g (m s −2 ) gravitational acceleration, v meridional wind velocity (m s −1 ), and p atmospheric pressure (hPa). Using both an IWV and vIVT scheme in this study captures some of the variability between different approaches of AR detection. The IWV scheme is slightly advantageous for studying the effects of ARs on surface melting since excessive IWV leads to the development of mixed phase clouds containing anomalously high liquid and ice water paths with higher radiative fluxes 20 . However, the meridional wind component in the vIVT scheme captures some of the dynamical processes that lead to foehn wind generation and subsequent sensible heat fluxes on the AP ice shelves. While the vIVT generally detects more ARs than the IWV scheme, controlling for maximum IVT upon landfall greatly reduces this disparity and shows that differences in AR detection is mainly influenced by the depiction of weaker ARs (Fig. 1b, d).
The AR detection algorithm in this study was configured to be included in the Atmospheric River Tracking Method Intercomparison Project 16 (ARTMIP). This entails running the detection algorithm on fields from the MERRA-2 on all pressure levels and outputting the AR shapes in binary form. For the sake of comparison, the AR detection algorithm was also run on ERA-5 IWV and vIVT fields on a 1 × 1 degree grid. The previous studies of ARs over Antarctica confirmed that the AR detection frequency is similar amongst various modern reanalysis products 20,23 and it is confirmed again in this study between MERRA-2 and ERA-5 (Fig. 1). The analyzes presented in this manuscript were primarily computed using the vIVT scheme applied to MERRA-2 reanalysis unless explicitly stated otherwise. AR intensity is measured through different metrics depending on the variable being analyzed. Instantaneous intensity which is useful for measuring a corresponding instantaneous variable (i.e., 3-hourly temperature extreme, melt rate, and runoff rate seen in Figs. 3, 5, 7) is measured from the maximum IVT value simulated within the shape of the AR over the AP domain (see Fig. 1e). Throughout the manuscript, a maximum IVT threshold above 400 kg m −1 s −1 is used to describe intense ARs, which also reduces the spread in AR frequency estimation among the MERRA-2, ERA-5, vIVT, and IWV detection schemes (Fig. 1). Cumulative intensity which is useful for measuring the cumulative effects of an AR landfall (i.e., iceberg calving frequency seen in Fig. 6b and the relationship between annual cumulative AR intensity and melt/runoff seen in Fig. 6c) is the summation of maximum IVT over the course of a distinct AR landfall event thus incorporating the duration and IVT strength of the storm. Daily cumulative intensity is like the cumulative intensity except we make the summation of maximum IVT for daily periods when an AR landfall is detected. This is used for testing the probability of an extreme state (i.e., temperature extreme, melt extreme, runoff extreme, 2-day sea-ice removal, high swell seen in Fig. 6a) Regional climate model and reanalysis datasets. MAR is a regional climate model specifically designed for simulating polar climate 24,66 . MAR atmospheric dynamics are based on the hydrostatic approximation of the primitive equations 67 . The exchanges between the atmospheric part of MAR and the surface are handled by the complex energy-balance snow model SISVAT 68 , based on CROCUS 69 that explicitly simulates 30 layers resolving the 20 first meters of snow or ice. The surface module notably represents percolation of meltwater and its retention into the snowpack. Runoff occurs when the snowpack can no longer absorb additional liquid water (i.e., snowpack water content exceeding 5% or surface liquid water over bare ice or an ice-lense layer). MARv3.11 was run at a resolution of 7.5 km and was forced by 6-hourly outputs of the latest ERA-5 reanalysis between 1979 and March 2020. The first year was discarded as spin-up.
MAR evaluation. MAR has been thoroughly compared to observation from weather stations, satellite melt extent, AWS-forced melt estimates and SMB measurements over the Antarctic ice sheet 24,70-72 and more specifically over the Peninsula 11,73 . The model is also abundantly evaluated over the Greenland ice sheet 74 . The ability of MAR to represent melting over Antarctic ice shelves has already been discussed in previous research 71,73,75 , so we present here a short comparison with AWS observations and melt estimates. MAR correctly represents the variation in near-surface temperature ( Supplementary Fig. 20a) with a centered RMSE (i.e., RMSE where we removed systematic biases attributed to elevation difference) of 2.6°C over the ice shelves and 3.2°C over the grounded ice. The comparison also reveals a mean negative bias. MAR tends to slightly underestimate the high (positive) temperatures while overestimating the low temperatures that can be found over higher elevations or in winter. Despite the high resolution of these simulations (7.5 km), the altitude of pixels relating to AWS near the margins is overestimated in MAR by slightly less than 100 meters, likely explaining a major part of the negative bias over low-elevation areas and ice shelves. There is no direct observation of melt, so we compared MAR to melt estimates using an AWS-forced snow model 76 . Only AWS locations whose elevation difference with MAR does not exceed 250 m were selected. This comparison (Supplementary Fig. 20b) reveals that the model correctly reproduces the annual surface melt at the AWS situations (located on the Larsen C and Scar Inlet of Larsen B), except in 2017 at Larsen C West where MAR overestimated melt.
Sea ice, wave, and back trajectory data. ERA-5 sea-ice products are developed under the Ocean and Sea Ice Satelite Application Facilities (OSI-SAF) that include a reprocessed version (OSI-409-a), which is used by HadISST2 until 2007 and then switches to an operational version (OSI-SAF) afterwards 77 . We attempted to examine sea-ice thickness changes following AR events using SMOS Level 3 Sea Ice Thickness product, which provides daily estimates of SMOS-retrieved sea-ice thickness around the Antarctic continent. However, the possible range of retrievable ice thickness values decreases when the brightness temperature of the sea ice becomes saturated, which often occurs during AR events and could result in an underestimation of the true ice thickness (Kaleschke, L., personal communication).
ERA-5 includes a two-dimensional wave model coupled to the atmosphere. It solves the wave energy-balance equation, using wave spectra with 24 directions and 30 frequencies, and includes assimilation of altimeter wave height 78 . Although advanced methods to represent wave-sea-ice interactions have been tested at the European Centre for Medium-Range Weather Forecasts 78 (ECMWF), their representation in ERA-5 is very simple: the wave energy is set to zero where the sea-ice concentration exceeds 30%, and there is no wave attenuation for lower concentrations (Bidlot, J. R., personal communication).
The origin and trajectory of air masses were evaluated using the Lagrangian Particle Dispersion Model Flexpart 10 79 . A batch of 500 neutral inert air tracer particles are randomly released from a volume (0.1°x 0.1°x 100 m) around the −62°E −66.6°N coordinates at a 2000 m altitude. Flexpart is driven by 249 3-hourly meteorological fields at 1°x 1°horizontal resolution from ERA-5 (downloaded using the flexextract tool 80 ) to compute 10-day back-trajectories.
Satellite observations. MODIS-Terra Aqua images were acquired at the (https:// worldview.earthdata.nasa.gov) over the entire period of availability, from the end of August to the beginning of April from between February 25, 2000 and April 28, 2020. The worldview interface allows comparing images between two different dates, making it possible to retrieve the exact occurrence of the main calving events and collapses of the ice shelves when cloud interference is limited. Since the foehn effect induces a reduction of the cloud cover at the east of the Antarctic Peninsula, it was possible to retrieve the approximate occurrence of most events of the Larsen A, B and C ice shelves. Accurate examination of images before winter (end of April) and after winter (beginning of August) allowed us to verify that most large ice-shelf calving/collapse on the Larsen ice shelves occurred during the summer in the 20 years of observation (see images in supplementary material for examples). It is important to note that the dates we link with calving/collapse events are subject to uncertainty based on when sky conditions allow us to observe these events.
Since daily MODIS images were not available before February 25, 2000, collapses and large calving events occurring before this date were not considered in our analysis. Nevertheless, we analyzed the well documented collapse of the Larsen A ice shelf 15 using AVHRR satellite images available on the National Snow and Ice Date Center (NSIDC) website 81 surrounding the date of the Larsen A collapse.
Winter calving events were not retrieved with the same accuracy because even though thermal bands could be used during polar night, this imagery technique is more sensitive to cloud cover (even thin cloud are opaque). This lowers the resolution even given favorable cloud cover, thus making the acquisition of data at a daily time scale not straightforward. Even if we present the Larsen C calving event of July 10-12, 2017, we did not consider this event in the final statistics. This prompted us to consider the iceberg retrieval database from Stuart and Long, 2011 82 . Here, we extracted the initial location of icebergs around the Larsen ice shelves, but difficulty was encountered in tracing with confidence the origin of these icebergs to a given ice-shelf because the first detection of an iceberg can occur days after calving and its position can be far from its parent ice-shelf. Nevertheless, the information of proximity with an ice shelf provides some hints on calving in the Larsen area. In the latter database, winter occurrences are also largely less frequent.
Similar examination was more limited over the Wilkins ice shelf because cloud cover is generally denser, so collapses and calving events are generally observed a few days after their occurrence. Moreover, contrary to the eastern side, several large events that occurred in the darker autumn and winter periods were not visible in MODIS datasets. This made the attribution of ARs to calving occurrences more difficult on the western side of the AP.
Surface melt surface anomalies over the Larsen ice shelves were also retrieved between October 1 and March 31 every year over 1979-2020 83 . A series of microwave radiometer observations (namely Scanning Multichannel Microwave Radiometer (SMMR), the Special Sensor Microwave/Imager (SSM/I) and Special Sensor Microwave Imager/Sounder (SSMIS)) were processed to test our conclusions on AR impact on snow melt (see Supplementary Fig. 7). First, a climatology of daily melt surface extent based on mean values for each Julian day mean was computed over the full 41 year-long period. The daily melt anomalies were computed by comparison with this climatology and compared with AR occurrences. Data obtained before 1993 must be considered with caution, in particular with summer 1987-88 31 Finally, to complement our analysis on the relationship between melt pond observations and AR activity obtained using the MODIS image interpretation 34 , we also studied the Landsat8 image database 84 . Information on lake extent on the Larsen ice shelves were retrieved and merged over periods of 16 days (corresponding to the Landsat revisit time). Far larger melt pond areas were usually observed at the end of February and in March, but after a comparison with Sentinel 2 images over the same periods, we concluded that this very likely reflects artifacts resulting from the misclassification of scattered cloud shadows as lakes by the algorithm 84 . This issue is more prominent at the end of the summer when the zenith angle is high and the shadows of scattered clouds are far away from the cloud vertical, maximizing the areas of misclassification. This artifact impedes an accurate estimation of melt pond areas in February and March, i.e., when melt ponds should be most impactful on hydrofracturing processes. While improving the algorithm to address this issue seems worthwhile, it is beyond the scope of this study as it would require reprocessing the full Landsat8 archive on Larsen ice shelves.
Statistical significance of AR-extreme event co-occurrences. To confirm that AR and collapses were not independent, we considered various 1000 random selections of dates between August and March from 2000 to 2020, and we find that the probability of having an AR within the 5 preceding days is 24.5 %. If ARs and collapses were independent, then the number of collapses (among the 21) that would occur after ARs (within the 5 preceding days) would follow a binomial distribution with parameters n = 20 and a probability p = 0.245. Considering that 13 out of 21 events co-occurred with an AR within the 5 preceding days, this leads to a p-value of p = 5.3 × 10 −5 . In other words, there is less than 1/1000 chance to experience 13 AR events or more if AR and collapses were independent. We applied the same approach to retrieve the level of significance of co-occurrences between AR and extreme temperature, melt, runoff, and swell height. For composite AR-related rainfall and AR-related 2-day sea-ice change, we tested the equality of the means of AR days and non-AR days with a Student's t-test (see main text for results).

Data availability
The Modèle Atmosphérique Régionale (MAR) data is available at https://doi.org/ 10.5281/zenodo.6347190. We acknowledge the use of imagery from the NASA Worldview application (https://worldview.earthdata.nasa.gov), part of the NASA Earth Observing System Data and Information System (EOSDIS). NOAA AVHRR imagery was acquired from the NSIDC 81 . The output from the AR detection algorithm can be downloaded from the ARTMIP database at https://www.earthsystemgrid.org/dataset/ ucar.cgd.ccsm4.artmip.tier1.html.

Code availability
The code for the AR detection algorithm discussed in this paper is available at https:// doi.org/10.5281/zenodo.4009663.