Climatology and decadal changes of Arctic atmospheric rivers based on ERA5 and MERRA-2

We present the Arctic atmospheric river (AR) climatology based on twelve sets of labels derived from ERA5 and MERRA-2 reanalyses for 1980–2019. The ARs were identified and tracked in the 3-hourly reanalysis data with a multifactorial approach based on either atmospheric column-integrated water vapor (IWV) or integrated water vapor transport (IVT) exceeding one of the three climate thresholds (75th, 85th, and 95th percentiles). Time series analysis of the AR event counts from the AR labels showed overall upward trends from the mid-1990s to 2019. The 75th IVT- and IWV-based labels, as well as the 85th IWV-based labels, are likely more sensitive to Arctic surface warming, therefore, detected some broadening of AR-affected areas over time, while the rest of the labels did not. Spatial exploratory analysis of these labels revealed that the AR frequency of occurrence maxima shifted poleward from over-land in 1980–1999 to over the Arctic Ocean and its outlying Seas in 2000–2019. Regions across the Atlantic, the Arctic, to the Pacific Oceans trended higher AR occurrence, surface temperature, and column-integrated moisture. Meanwhile, ARs were increasingly responsible for the rising moisture transport into the Arctic. Even though the increase of Arctic AR occurrence was primarily associated with long-term Arctic surface warming and moistening, the effects of changing atmospheric circulation could stand out locally, such as on the Pacific side over the Chukchi Sea. The changing teleconnection patterns strongly modulated AR activities in time and space, with prominent anomalies in the Arctic-Pacific sector during the latest decade. Besides, the extreme events identified by the 95th-percentile labels displayed the most significant changes and were most influenced by the teleconnection patterns. The twelve Arctic AR labels and the detailed graphics in the atlas can help navigate the uncertainty of detecting and quantifying Arctic ARs and their associated effects in current and future studies.


Introduction
The Arctic region has warmed at rates double that of the global average over the past decades (Polyakov et al 2002, Serreze and Francis 2006, Graversen et al 2008, Serreze et al 2009, Screen and Simmonds 2010, Serreze and Barry 2011, Cohen et al 2014. This phenomenon is widely known as Arctic amplification. The rapid warming and associated sea ice decline, as well as glacier retreat, have had significant ecological and economic impacts locally and globally (Francis and Vavrus 2012, Post et al 2013, Cohen et al 2014, Stroeve et al 2014, Overland et al 2016, Crépin et al 2017, Francis et al 2017. There is a wide range of physical processes underlying the Arctic amplification. Among them, the most often cited mechanism is the Arctic sea ice decline and associated albedo feedback (Hall 2004, Serreze et al 2009, Screen and Simmonds 2010, through which initial warming tends to melt sea ice, exposing more open water that readily absorbs more solar radiation, leading to further warming. However, it has been shown that surface albedo feedback is not a dominating mechanism (Winton 2006). There are other important dynamical and radiative feedbacks (e.g. Alexeev et al 2005, Graversen et al 2008, Abbot et al 2009, Chylek et al 2009, Graversen and Wang 2009, Shindell and Faluvegi 2009, Serreze and Barry 2011, Graversen and Burtu 2016.
In particular, recent works on poleward moisture transport (e.g. Woods et al 2013, Woods and Caballero 2016, Johansson et al 2017, Liu et al 2018, Gimeno et al 2019 have shown it to be one of the factors that could influence Arctic warming and sea ice decline in the boreal winter through enhanced downward longwave radiation and cloud radiative effects. Atmospheric rivers (ARs), the narrow filaments of enhanced moisture transport typically associated with a low-level jet and extratropical cyclone , is a conspicuous pathway for poleward moisture transport (e.g. Neff et al 2014, Liu and Barnes 2015, Hegyi and Taylor 2018, Komatsu et al 2018, Mattingly et al 2018. The roles of ARs in Arctic and high-latitude climate have received growing attention in recent years (e.g. Baggett et al 2016, Guan et al 2016, Hegyi and Taylor 2018, Liang and Sushama 2019, Xiao and Lettenmaier 2021 due to the enhanced moisture that increases the downward infrared radiation, cloud-radiative effects, precipitation and precipitation process, and changes of surface energy balance, potentially leading to lasting surface warming in the Arctic.
We can therefore hypothesize that the changes in Arctic AR's frequency of occurrence, strength, and location are integral to the warming and surface melting in Arctic amplification. However, the Arctic is an observational data-sparse region, which has posed significant uncertainty in constructing the Arctic AR climatology. For example, the changes in the Arctic atmospheric moisture content and moisture transport have remained ambiguous in the literature. They appeared dependent on the studied region, time, elevation, and the reanalysis data (e.g. Zhang et al 2013, Dufour et al 2016, Woods and Caballero 2016, Villamil-Otero et al 2018, Gimeno et al 2019, Rinke et al 2019, Nygard et al 2020. The disagreements prompted us to answer these questions: what is the feature of climate change in Arctic AR's frequency of occurrence? Is there a change of moisture transported by ARs into the Arctic? Furthermore, is there a change in the proportion of poleward moisture transport attributed to ARs? Finally, are there recurring spatial patterns of ARs? In probing these questions, we created twelve sets of Arctic AR labels with ERA5 and MERRA-2 reanalysis from 1980 to 2019. Essentially, the Arctic AR detection and tracking algorithm was designed to be consistent with that developed in the mid-latitudes (Zhang et al 2021) since the AR labels were used to study the moisture transported by the ARs into the Arctic, while there is no physical boundary between the Arctic and the lower latitudes.
Even more, beyond the current work, the labels should be useful for studying the sensible and latent heat transports as well as the radiative and cloud-radiative effects associated with the Arctic ARs. These considerations are different from some of the previous AR identification and tracking works. For example, (Guan and Waliser 2015) applied an IVT minimum threshold of 100 kg m −1 s −1 to supplement the seasonal 85th percentile IVT climate thresholds when tracking global ARs. The lower limit of 100 kg m −1 s −1 was selected to isolate AR-like features in the polar regions and to be consistent with the AR detection results in East Antarctica by Gorodetskaya et al (2014) that captured Antarctic extreme precipitating events. Moreover, (Gorodetskaya et al 2014) integrated IWV and IVT from 900 to 300 hPa. The exclusion of the layer between the surface and 900 hPa was to avoid the katabatic flow near the continental edge. Similarly, (Wille et al 2019) also removed the lower boundary layer condition of the atmosphere (i.e. IWV and IVT) below 900 hPa to avoid the katabatic wind effects from the ice sheet plateau when constructing AR climatology over Antarctica. The seminal Arctic moisture instruction studies by Woods et al (2013) and Woods and Caballero (2016) defined Arctic moisture intrusions as intense northward moisture transport within the winter season. (Wille et al 2021) and (Maclennan et al 2022) also used the meridional component of the integrated water transport to detect continuous filament of extreme AR events with long meridional spans.
However, the physical considerations in Antarctica might not apply to the Arctic. Over the Arctic, AR-related poleward sensible and latent heat transport plays an essential role in the Arctic climate , Liu and Barnes 2015, Baggett et al 2016, Komatsu et al 2018, Mattingly et al 2018. Moreover, Arctic ARs increased cloud amount and higher cloud liquid and ice water contents, affecting subsequent surface energy budgets via cloud-radiative effects (Hegyi and Taylor 2018, Liu et al 2018, Box et al 2022. Besides, the topography of the Arctic is different from the Antarctic. The elevated Greenland ice sheet only constitutes a small fraction of the Arctic, which consists mostly of ocean and sea ice. Thus, excluding the lowest portion of the atmosphere in the AR detection algorithm may be inappropriate for detecting Arctic ARs, characterizing associated AR climatology, and examining AR-related clear-sky radiative and cloud-radiative effects. Similarly, the extent of the AR effects on surface energy budgets should be examined beyond strong poleward-only moisture transports. Last but not least, the radiative forcing of water vapor scales logarithmically with the specific humidity (Raval andRamanathan 1989, Huang andBani Shahabadi 2014). So, water vapor's greenhouse effect increases the most with a slight uptick of specific humidity in originally low-moisture places such as the Arctic, where cold sea surface and ice packs dominate. Setting a fixed minimum threshold on IVT or IWV could preclude a comprehensive study of AR's radiative forcing.
Henceforth, we detected and tracked Arctic AR events that might be associated with significant heat and moisture transports or distinctly accompanied by extensive radiative effects using the multifactorial approach in Zhang et al (2021). We then examined and discussed the AR climatology, trends, and decade-by-decade changes in the Pan-Arctic region (figure 1). In the remainder of the paper, the data and methods are described in section 2; the analysis and discussions of AR occurrence frequency, moisture, moisture transport, trends, and modulation by the teleconnections with internal climate variability are in section 3; section 4 is the conclusion. Appendix A details the AR detection and tracking algorithm, while appendix B provides additional figures supporting the discussions in the main text. An extensive set of figures that visualizes and compares these AR labels is presented in the Supplementary materials.

Data
Reanalysis has been considered as a comprehensive data in the Arctic due to the sparse in situ observations. For the atmospheric moisture transport, reanalyses have been found to be qualitatively consistent with the radiosonde data in terms of spatial and temporal patterns, even though a slight overestimation by the reanalyses (Dufour et al 2016). In addition, because (Nayak and Villarini 2017) recommended high-resolution products for AR impact assessments, our analysis focused on the two highest-resolution reanalysis data, the fifth generation of ECMWF atmospheric reanalysis (ERA5, (Hersbach et al 2020)) and NASA Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2, (Gelaro et al 2017)).
At the start of this project, an ERA5 source dataset for the AR Tracking Method Intercomparison Project (ARTMIP; (Shields et al 2018, O'Brien et al 2020, Collow et al 2022) was not available. Hence, we calculated the two conventional moisture measurements, integrated water vapor transport (IVT) and integrated water vapor (IWV), for AR detection and tracking in ERA5, specifically following the definition used by the ARTMIP in their creation of the MERRA-2 source data for ARTMIP (NCAR CDG 2019). This MERRA-2 source data was calculated by the University of California, San Diego (UCSD) Center for Western Weather and Water Extremes according to the equations (e.g. Shields et al 2018): The three variables, horizontal wind (V h = (u, v) where u is the zonal and v the meridional winds in m s −1 ), specific humidity (q in kg kg −1 ), and pressure (p in hPa), used in the equations were from NASA MERRA-2 (Gelaro et al 2017). MERRA-2 is the latest atmospheric reanalysis of the modern satellite era produced by NASA's Global Modeling and Assimilation Office. It made progress on the assimilation of modern hyperspectral radiance and microwave observation, as well as the GPS-Radio Occultation datasets (Gelaro et al 2017). The horizontal spatial resolution and the temporal resolution of the vertically integrated fields in the MERRA-2 source data for ARTMIP are 0.625 • longitude by 0.5 • latitude and 3 hours, respectively. The ARTMIP used the 3-hourly reanalysis data as a compromise between computational efficiency and sufficient temporal information for AR detection (Shields et al 2018).
ERA5 is the latest version released by ECMWF in 2019 with high spatial (0.25 • longitude by 0.25 • latitude) and temporal (hourly) resolutions. ERA5 is based on ECMWF Integrated Forecasting System (IFS) CY41r2 by 4D-Var data assimilation and model forecasts with 137 hybrid sigma/model levels to 1 Pa in the vertical (Hersbach et al 2020). Besides, Atmospheric data in ERA5 on the standard pressure levels from 1000 hPa to 1 hPa (i.e. a total of 37 vertical pressures) are interpolated from the hybrid sigma/model levels in the IFS.
However, surface pressures over the Arctic may be lower than 1000 hPa at high altitudes (e.g. Greenland, East Siberia, Alaska). Therefore, considering Arctic topography, we integrated the IVT and IWV from the surface pressure, P sfc , to 200 hPa using the modified (1) and (2) at the same 3-hourly time steps as those in the MERRA-2 source data. The corresponding equations are shown in (3) and (4) as follows: We generated the twelve Arctic AR labels from January 1980 to December 2019 using the IVT and IWV derived from either ERA5 or MERRA-2 with the algorithm described in the next section. In the subsequent analysis of the ARs detected by these labels, we incorporated the analysis on the surface 2-m temperature (T2m) from both reanalyses. To facilitate comparison between the AR detected in ERA5 and those in the MERRA-2 source data, we down-sampled the ERA5 data from hourly to 3-hourly at 00, 03, 06, 09, 12, 15, 18, 21 UTC each day from January 1980 to December 2019. Moreover, the fields from the 0.625 • × 0.5 • MERRA-2 were interpolated bilinearly to ERA5's 0.25 • × 0.25 • horizontal mesh for comparison. In addition, we extracted the 25 km×25 km monthly sea ice concentration (SIC) data from 1980 to 2019 from the National Snow and Ice Data Center (National Snow and Ice Data Center 2020). The daily indices of Arctic Oscillation (AO), North Atlantic Oscillation (NAO), and Pacific/North America (PNA) pattern were retrieved from (NOAA National Weather Service Climate Prediction Center 2021). The monthly indices of Interdecadal Pacific Oscillation (IPO) and AMO (Atlantic Multidecadal Oscillation) were both extracted from NOAA Physical Sciences Laboratory, specifically at (NOAA Physical Sciences Laboratory 2023b) and (NOAA Physical Sciences Laboratory 2023a), respectively.

AR detection and tracking algorithm
The AR detection and tracking algorithm followed a similar multifactorial design as (Zhang et al 2021), in which we found that the choices of moisture fields and climate thresholds strongly dictated the summary statistics of the identified ARs and AR-related surface hydrometeorological effects. In this work, we created . Schematic diagram illustrating the Arctic AR detection and tracking algorithm using factors including two moisture fields (IVT and IWV) and three levels of climate thresholds (75th, 85th, and 95th percentiles) with ERA5 and MERRA-2 source data for ARTMIP. The two light blue frames in the flowchart indicate the two main steps in the algorithm. The first step is the AR object detection with the three levels of climate thresholds and the geometry threshold (based on Zhang et al (2021)). The second light blue frame illustrates the AR event tracking with a Lagrangian AR tracking routine (modified based on Guan and Waliser (2019)) as well as the event penetrating latitude and duration threshold. For each set of AR labels, a total of five AR events (red contour) were detected over the Arctic at this time. The five AR events are marked as '1' , '2' , '3' , '4' , and '5' in each subplot. six sets of Arctic AR labels with each reanalysis dataset from January 1980 to December 2019. That is, a total of twelve sets of labels were created based on the combinations of conditions depicted in figure 2.
The schematic diagram in figure 2 summarizes the Arctic AR detection and tracking algorithm. The algorithm is composed of two steps: AR object detection (shown in the first light blue frame) and AR event tracking (the second light blue frame). The first step of AR object detection was based on the IVT or IWV field exceeding one of the three levels of monthly climate thresholds (75th, 85th, and 95th percentiles) at each grid-point location in each reanalysis dataset. Moreover, to be qualified as an AR event, the extracted objects must satisfy a geometric constraint (length ⩾ 1500 km and length: width ratio ⩾ 2). In the second step, we examined the time sequence of identified AR objects to construct AR events with the Lagrangian tracking framework based on the spatial proximity and morphological similarity between two consecutive time steps (Guan and Waliser 2019), as schematically shown in the second light blue frame in figure 2. Particularly, we kept only AR events penetrating the Arctic (the boundary of the Arctic is defined as the 60 • N) and continuously existed for at least 18 hours in the Arctic. In the end, a total of twelve sets of Arctic AR labels expressing spatial and temporal AR information based on the factor of moisture fields and climate thresholds from January 1980 to December 2019 were completed with ERA5 and MERRA-2 reanalyses. The algorithm is further detailed in appendix A.
Figures 3(a) and (b) are examples of the Arctic ARs identified by the 85th_IVT-and the 85th_IWV-based AR labels, respectively, at the same time on 30 May 2017. Both AR labels identified five AR events. The locations of four AR events were roughly consistent between two AR labels, corresponding to the respective IVT or IWV maximum regions. Both labels detected the most prominent Arctic AR event, the long narrow filaments of enhanced IVT or IWV from the Eurasian continent across the Barents-Kara Seas and towards the Greenland Sea, marked as Event 1. According to the 85th_IVT-based labels (figure 3(a)), this AR event persisted over the Arctic for 2.5 days from 29 May to 1 June 2017. The 85th_IWV-based labels (figure 3(b)) detected the event in the Arctic 2 days earlier since 27 May for 4.5 days. The second AR (Event 2 in figure 3) identified by both AR labels occurred from Western Canada to the Beaufort Sea, with more extensive AR coverage determined by the 85th_IWV-based AR label. The IWV-based labels (figure 3(b)) determined a more extensive coverage than the IVT-based labels (figure 3(a)). Also, Event 2 identified by the IVT-based labels lasted about four days in the Arctic, in contrast to more than six days according to the IWV-based label.
In the Atlantic sector, the 85th_IVT-based AR labels captured one AR over the Denmark Strait (Event 3 in figure 3(a)) and the other (Event 4) from Western Europe, Scandinavia, to the Arctic. However, in the IWV field, these two AR events merged as one AR event (Event 3 in figure 3(b)) according to the 85th_IWV-based AR label. This merged AR event stemmed from Western Europe, penetrated Scandinavia, and moved northward towards the Denmark Strait. Moreover, at the same time, the 85th_IVT-based labels picked up a small area of an AR over East Siberia (Event 5 in figure 3(a)), which had a much more substantial presence with the IWV-based label. In both labels, this was the end of a long event persisting over East Siberia since 21 May 2017. One additional AR event (Event 4 in figure 3(b)) was detected only by the 85th_IWV-based labels at this time step. The event came from North America and ascended into the Arctic over Hudson Bay.
AR-related IVT and IWV magnitudes over the Arctic can be comparable to those over the mid-latitudes. As shown in figure 3, the AR-related maximum IVT and IWV values were greater than 350 kg m −1 s −1 and 27 kg m −2 , respectively. Furthermore, the 85th_IWV-based AR labels identified larger, broader, and longer persistent AR events than those using the 85th_IVT-based label. Interestingly, from 1980 to 2019, the largest areas ever covered by an Arctic AR instantaneously were approximately 6.5 and 8.5 times the size of Greenland according to the 85th_IVT-and the 85th_IWV-based label, respectively, on 15 November 1996, at 06 UTC and 31 May 2012, at 21 UTC. Figure 4 visualizes the time series of annual AR event counts in the Arctic, north of 60 • N, based on the 12 AR labels. Overall, in both ERA5 (figure 4(a)) and MERRA-2 (figure 4(b)) reanalyses, regardless of the IVT-or IWV-based moisture field, the AR annual events decreased dramatically as the climate threshold increased from the 75th to the 95th percentile. With the IVT-based field for AR detection, the 75th percentile climate threshold-based AR labels from ERA5 and MERRA-2 yielded a median of 559 counts year −1 . The corresponding 85th-and 95th-based labels resulted in median values of 476 counts year −1 and 267 counts year −1 . The IWV-based labels led to median values of 428, 406, and 246 counts year −1 from the 75th, the 85th, to the 95th percentile climate thresholds, respectively. Nonetheless, slightly more AR events were identified in ERA5 than in MERRA-2 when factors of moisture field and climate threshold are controlled.

AR event statistics
Despite the different magnitudes, all labels except for the 75th_IWV-based MERRA-2 labels showed consistent upward trends of AR events from the mid-1990s to 2019, as illustrated by the dashed lines of the loess fits (locally-weighted regression, (Cleveland 1979)) with a window size of 30 years. The trends are especially prominent in recent decades. Such inflection is most evident in the ERA5 IVT-based labels. It is noted that with the IWV-based labels from both reanalyses, the upward trends were increasingly steeper as the climate threshold increased from the 75th to the 95th percentile, suggesting a more rapid increase of extreme AR events within this group. Accompanied by the increasing trends, all of the AR event count time series display quasi-biennial oscillatory signals over the 40 years. As will be discussed in section 3.4, the Arctic ARs activities were strongly influenced by the Arctic Oscillation. Since Quasi-Biennial Oscillation (QBO) has been known to impact on the spatial structure of the AO Yong 2021, Cai et al 2022), the oscillatory signals in AR event time series might be associated with QBO. Figure 4 shows that IWV-based labels resulted in fewer Arctic AR events than IVT-based ones. However, AR events identified with IWV tend to last longer. The boxplots in figure 5(a) summarize the sampling distributions of the AR per-event duration (in days) in the Arctic across the entire 40 years based on the 12 AR labels. Figure 5(b), similarly, shows the distributions of the areas swept by AR events resulted from the 12 AR labels (in m 2 ). An area swept by an AR event was determined with all the grid points to the north of 60 • N visited at least once by the AR throughout its event lifetime, consistent with the definition used previously in Zhang et al (2021), c.f. their figure 9. The values of AR per-event duration (figure 5(a)) and swept-area (figure 5(b)) shown in the boxplots are the original values transformed with base-2 logarithm to accommodate the wide range.
It is noted that in figure 5(a), more restrictive climate thresholds did not always yield shorter per-event duration. Regardless of the moisture field or the reanalysis data used, the 85th percentile climate threshold-based labels surpassed the 75th percentile-based ones and presented the longest median and Q3 of per event duration. This may be due to the 75th percentile thresholds permitting the detection of more splitting of ARs into shorter events. Nevertheless, the effects of climate thresholds appear to be secondary to the moisture detection fields for AR detection. The IWV-based labels identified longer per event duration than IVT-based ones when holding the climate thresholds constant. This conspicuous feature can be seen as x-axis). Each figure of the data sets has two packets, or subsets of values conditional on different moisture field for detection (IVT or IWV). Each packet was grouped by color into three levels of climatological thresholds (75th, 85th, and 95th percentiles) shown respectively in blue, red, and purple. The dashed black lines show the loess fits (locally-weighted regression) with a window size of 30 to the annual Arctic AR events across the 40 years detected by each AR label.

Figure 5.
Boxplots of base-2 logarithmic transformation of (a) AR per-event duration (day) and (b) AR per-event swept area (m 2 ) according to the IVT and IWV-based AR labels conditional on reanalysis data (ERA5 or MERRA-2) and climatological thresholds labeled as percentile in 75th, 85th, or 95th. Each boxplot includes the colored box spanning from the first quartile (Q1) to the third quartile (Q3) of the distribution, a black dot marking the second quartile (median, Q2), and the whiskers. The whiskers extend to the most extreme data point that is no more than 1.5 times the length of the box interquartile ranges (IQR) away from the box. Any data points outside the whiskers are marked as potential outliers in light blue.
the larger Q1, Q2, Q3, and the maximum values of each subset with IWV-based labels in figure 5(a). Furthermore, for each set of AR labels in the two reanalyses, the distribution of per-event duration across each year did not appear to have a distinguished upward trend (figures not shown). Moreover, as expected, the areas swept by ARs decreased with rising climate thresholds on both moisture fields in both reanalyses. Analogously, the median per-event average swept-area corresponds to about 141%, 113%, and 80% of the size of Greenland using the 75th-, 85th-or 95th-percentile climate thresholds in both reanalyses (figure 5(b)). ARs detected by the 75th-percentile-based labels influenced discernibly broader areas in the recent decades (supplementary figure S1), reflected in a positive shift in the sampling distribution over time. The figure S1 illustrates this change with annual boxplots. The 75th_IWV-based labels exhibit the most obvious AR broadening trends, in contrast to their lack of trend in event counts (figure 4). To a lesser degree, the broadening is also seen in the ARs identified by the 85th_IWV-based labels in figure S2. There is no apparent broadening of AR-affected areas with 85th_IVT-, 95th_IWV-, or 95th_IVT-based labels in figures S2 and S3.
Overall, the statistics of AR event occurrences presented in this section were more sensitive to the selection of climate threshold and moisture field than to the choice of reanalysis data. More restrictive climate thresholds, such as the 95th percentile, emphasized extreme AR events with limited occurrences, smaller affected areas, and shorter duration in the Arctic. These results were mostly consistent with the conclusions derived from the ARs making landfall on the US continents (Zhang et al 2021). Furthermore, even though ERA5 had more AR events than MERRA-2. Given the same climate threshold and moisture field (figure 4), both reanalyses showed an increasing trend in AR event counts in the recent decades. The recent trend was the most pronounced when the 95th_IWV-based labels were used. The 75th percentile IVT-and IWV-based labels, as well as the 85th IWV-based labels, detected some broadening of AR-affected areas over time, while the rest of the labels did not. It is reasonable to hypothesize that the 75th percentile labels and the 85th IWV-based labels are more responsive than the rest to Arctic warming manifested in the moisture field through the Clausius-Clapeyron relation.

Spatial distribution of AR frequency
At each grid point, AR occurrence frequency based on each set of AR label was calculated as the accumulated 3-hourly time steps in AR events during four decadal segments (1980-1989, 1990-1999, 2000-2009, and 2010-2019) and throughout 1980-2019, expressed as a percentage of total time steps in the designated time period. Figure 6 visualizes the occurrence frequency derived from the 85th percentile climate threshold-based labels in ERA5 and MERRA-2. The results derived from 95th-percentile labels are in appendix figure B1, while those from 75-percentile labels are in supplementary figure S4. The zonal mean of AR occurrence frequency as a function of latitude for all time segments are shown in figures 6(e) and (f) for the 85th percentile labels (cf 95th-and 75th-labels in figures B1(e) and (f) and S4(e) and (f), respectively). In both reanalyses, the climate threshold had a first-order influence on the AR occurrence. Regardless of the moisture field, the mean AR occurrence frequency decreased dramatically as the climate percentile threshold increased from 75th (10%-29% in figure S4), 85th (4%-18% in figure 6), to 95th (0.7%-6% in figure B1) during each analysis period. The choice of moisture field in AR detection rendered slightly different outcomes between the two reanalyses. For ERA5, the IWV-based labels detected overall more AR occurrence than the IVT-based labels for all analysis periods. In MERRA-2, the IWV-based labels resulted in slightly lower occurrence during 1980-1999 but much higher occurrence during 2000-2019 than the IVT-based labels. This rapid increasing trend is consistent with that in the event counts based on the 85th-and 95th-percentile labels in figure 4 The spatial distributions of ARs identified by the 12 labels were qualitatively similar (figure 6, cf figures B1 and S4). The maximum AR frequencies were over the Eurasian and North American continents from 1980-1989. In the following decade, the North Atlantic ARs increased. From 2000-2009, North Atlantic (Norwegian-Barents Seas) and northeast Canada, Davis Strait, and southern Greenland continued witnessing more ARs. In North Pacific, Bering Strait had more ARs according to the IVT-based labels, while the IWV-based labels detected more ARs over the Chukchi Sea. From 2010 to 2019, the maximum AR occurrence regions moved poleward to the Chukchi Sea, Barents-Kara Seas, and central Greenland with an increased magnitude and greater coverage than in previous decades. The maxima over the Arctic Ocean in the last two decades appeared to coincide, although without clear causality, with the retreat of September sea ice margins, an indicator of annual minimum sea ice concentration, in figure 6(g). Across the four decades, more ARs occurred over the Bering Strait, Davis Strait, southern Greenland, Norwegian-Barents-Kara Seas, and Siberia, which were roughly consistent with the main pathways for moisture transport into the Arctic found in the previous studies (Dufour et al 2016, Nygard et al 2020. Again, the latitudinal profile of zonal mean AR frequency (figures 6(e) and (f), cf (e) and (f) in figures B1, S4) suggested AR's poleward-moving and increasing trends: the peak occurrence frequency shifted into higher latitudes and increased from the first to the last decade. Most of the Arctic experienced increasing AR occurrence in the latest decade. The Figure 6. Maps of mean AR occurrence frequency (in %) in the Arctic during each decade (1980-1989, 1990-1999, 2000-2009, or 2010-2019) and the entire time period  identified by the 85th-percentile climate threshold applied to either IVT or IWV moisture field. Shown are results from (a) 85th_IVT-and (b) 85th_IWV-based AR labels with ERA5, and (c) 85th_IVT-and (d) 85th_IWV-based labels with MERRA-2. The corresponding latitudinal profiles of zonal mean frequency from (e) ERA5 and (f) MERRA-2 labels are shown on the right, where the grey area is bounded by the 2.5% and 97.5% quantiles of the zonal mean of AR occurrence frequency over 1980-2019. The subplot (g) shows the regions where the mean September sea ice concentration (SIC) was ⩾15% in each analysis period. extreme ARs represented via the 95th-percentile labels exhibited the largest increasing rates relatively, especially in the latest decade (e.g. figures B1 and B2). Figure 7 further verifies the conspicuous poleward shift of the maximum AR occurrence frequency from overland in the first two decades to the ocean in the last decade within the background trend of increasing AR frequency in the Arctic. It shows the longitudinal distribution of mean AR frequency of occurrence identified by the 6 ERA5 AR labels (top row) and the 6 MERRA-2 labels (bottom row) during each analysis decade at 60 • N and 75 • N, respectively. The land-ocean contrast of AR occurrence during each individual decade was observed at both latitudes. Along 60 • N, during the first two decades, AR occurred more frequently overland, i.e. the North American and Eurasian continents indicated by the grey segments along the x-axis in figure 7(a). In contrast, over the ocean, marked by the black segments in figure 7(a), lower AR occurrence frequency was observed. Within the following decade of 2000-2009, the maximum AR occurrence shifted eastward to the Labrador Sea and Greenland. In the last decade, the longitudinal maximum AR occurrence moved to the ocean with enhanced magnitudes, while the minima moved to the land. AR's different behaviors over land and ocean across the four individual decades were also observed at the 75 • N. Along 75 • N, the maximum AR occurrence was over the Eurasian continent (centered around 100 • E in figure 7(b)) in 1980-1989, then moved to the Chukchi-Beaufort Seas (around 180-120 • W in figure 7(b)), north of the North American continent (cf the peak AR occurrence in figure 7(a)) in 1990-1999. In recent two decades of 2000-2019, the peak AR occurrence shifted to the Arctic Ocean and its marginal seas (especially Chukchi, Barents-Kara, and Greenland Seas) with an even greater magnitude than those at 60 • N.
In figures 6 and 7, results from both reanalyses show that IWV-based labels performed differently from IVT-based labels in recent decades over the ocean and especially at 75 • N (figure 7(b)). The IWV-based labels detected more AR occurrences over and around the Arctic Ocean in recent decades, implying that IWV-based labels are not only more responsive to Arctic surface warming (as in figure 4) but specifically more sensitive to the surface warming over the ocean. Moreover, consistent with the results in section 3.1.1, Each subplot has 8 packets, or subsets of values conditional on 8 combinations of reanalysis data (ERA5 or MERRA-2 in the green boxes), and analysis period (1980-1989, 1990-1999, 2000-2009, 2010-2019 in yellow boxes).
Each packet has 6 paired values of mean AR occurrence frequency grouped by color into the 6 AR labels. To be consistent with the samples used in figure 8, data at each latitude were randomly sampled according to the sample number determined by twice the cosine value of the latitude to preserve the sample size:distance ratio. The grey and black segments at the base of each figure denote where the land and ocean regions are at each latitude.
Figure 8. Q-Q plots for AR occurrence frequency (in %) showing pair-wise comparison of AR labels in each of the eight Arctic sectors. Specifically, for each sector, (a) compares the probability distribution of AR occurrence frequency (%) between the 75th_IVT-and 85th_IVT-based AR labels by plotting their quantile values against each other in ERA5 (in red) and MERRA-2 (in blue). The 5%, 50% (median), and 95% quantile values were marked with the dotted red and solid blue lines for ERA5 and MERRA-2, respectively. The grey line on each panel represents a reference line with a slope of 1 and a y-intercept of 9%. Plots in (b) and (c) are similar to (a), except that (b) compares the distributions of AR occurrence frequency identified by the 95th_IVTagainst the 85th_IVT-based AR labels, where the reference lines have a y-intercept of −9%; (c) compares the distributions by the 85th_IWV-against the 85th_IVT-based AR labels, where the reference lines have a y-intercept of 0.
the spatial features of AR occurrence frequency were more sensitive to the selection of climate threshold and moisture field than the reanalysis data. However, the spatial analysis reveals further details about the quantitative differences between the two reanalyses. For example, among the IVT-based labels, more ARs were detected over a large portion of the central Arctic over the first 30 years in MERRA-2; among the IWV-based labels, ERA5 had more ARs during 1980-1999 and 2010-2019 (cf figures S5 and S6 for differences between pairs of ERA5 and MERRA-2 labels).

Empirical cumulative distribution comparisons in the most recent decade
We further quantified the relationship between the selection of climate threshold and moisture field in AR label generation and the resulted AR frequency of occurrence in the most recent decade of 2010-2019, when the maximum AR occurrences were observed. We divided the Arctic into eight equal sectors with an angle of 45 • , i.e. from 0 • -45 • E to 45 • W-0 • , to partition the spatial features of AR occurrence frequency seen in figures 6 and 7. For each Arctic sector, the quantile-quantile (Q-Q) plots in figure 8 compare the empirical cumulative distributions of AR occurrence frequency detected by the 85th-percentile IVT-based AR labels (85th_IVT) against those detected by the 75th-percentile IVT-based labels (75th_IVT, figure 8(a)), the 95th-percentile IVT-based labels (95th_IVT, figure 8(b)), and the 85th-percentile IWV-based labels (85th_IWV, figure 8(c)) by pairing the distributions according to their common quantiles. The comparisons were performed for labels from both ERA5 and MERRA-2. Within each octant, we randomly sampled along each latitudinal arc between 60 • and 85 • N with a sample size proportional to the cosine of the latitude. This way, the Q-Q plots of each Arctic sector in figure 8 do not over-estimate the AR's occurrence at higher latitudes and exclude the values above 85 • N. In figure 8, the grey solid line in each subplot has a slope of 1 and serves as a reference line for comparison. If two distributions have a similar spread, or standard deviation, the pairwise comparison of their common quantile values should be on or in parallel to the reference line. Due to the widely different magnitudes of AR occurrence frequency depending on the climate thresholds in figures 8(a) and (b), we vertically shifted the y-intercept of the reference line upwards or downwards from zero by 9%, respectively, to facilitate visual inspection. Each subplot also marks the 50% quantile values to show the median, along with the 5% and 95% quantiles to inspect the spread and the tails of a distribution.
Among the AR occurrence frequency detected by the 75th_IVT, 85th_IVT, and 95th_IVT labels, the higher the climate threshold, the smaller the median of the sample distribution, the narrower the distribution spreads, and the thinner its lower tail becomes. The Q-Q plots in figure 8(a) show that in all Arctic sectors, the 75th_IVT labels resulted in AR occurrence frequency with much larger values than the 85th_IVT labels, as the reference lines intercept the y-axis at -9%. The distributions of the 75th_IVT AR frequency have a larger spread than those of the 85th_IVT AR frequency since the majority of the 5%-95% quantile values form lines slightly steeper than the reference lines. The differences are minor but can be easily seen in the sectors including the Nordic Seas (0 • -45 • and 45 • -90 • E), Hudson Bay to Baffin Bay (45 • -90 • W), and North Atlantic Ocean and Greenland Sea (0 • -45 • W). Moreover, in figure 8(b), the reference lines intercept the y-axis at 9%, reflecting the smaller AR occurrence frequency from the 95th_IVT labels compared against that from the 85th_IVT labels. The most conspicuous feature in these Q-Q plots are the slopes formed by the 5%-95% quantile values. They are clearly smaller than 1, indicating a larger spread of the 85th_IVT AR frequency distributions than that of the 95th_IVT ones. Between the ERA5 labels, in all sectors except for 90 • -135 • W, the 85th_IVT AR frequency distributions had significantly thicker and longer lower tails than the 95th_IVT AR frequency distributions. This is seen with the lower tails rising toward the y-axis.
Furthermore, as seen in figure 6 with 85th_IVT labels and in the supplementary figure S4 with 75th_IVT, the AR occurrence frequency maxima in 2010-2019 were embedded in areas of large variances over the various Arctic oceanic regions. Specifically, in figure 8(a), the 75th_IVT labels resulted in maximum occurrence frequency over the Nordic Seas (0 • -45 • E in both reanalyses and 45 • -90 • E in MERRA-2) and the Chukchi-Beaufort Seas (135 • -180 • W). The 85th_IVT labels behaved similarly, although its maximum in the Chukchi-Beaufort Seas surpassed those in the Nordic Seas. However, with the 95th_IVT labels, the maxima were limited over the Pacific side from East Siberia to Chukchi-Beaufort Seas (135 • -180 • E and 135 • -180 • W) in figure 8(b) (cf figure S1). This is consistent with figure 4 discussed in section 3.1.1 and suggests that the lower climate thresholds are more sensitive to the Arctic surface warming and the ensuing moistening effects. And, the resulted empirical distributions of AR occurrence frequency have a wider spread and longer tails than the restrictive 95th-percentile climate threshold. It is however intriguing that, as the climate threshold increased, the maximum AR occurrence frequency shifted away from the Atlantic side to focus on the Pacific side of the Arctic. The signal on the Pacific side will be further explored in section 3.3.2. Figure 8(c), on the other hand, compares the quantile values between the 85th_IWV and the 85th_IVT AR occurrence frequency from each reanalysis in each Arctic sector. Here, the y-intercepts of the reference lines are zero. Hence, the upward displacements of the quantile values with respect to the reference lines indicate that the distributions resulted from the 85th_IWV labels are consistently larger than those from the 85th_IVT labels in both reanalyses and in all Arctic sectors. The majority of the 5%-95% quantile values form lines slightly steeper than the reference lines, indicating the AR occurrence frequency from the 85th_IWV labels had larger variance than those from the 85th_IVT labels. This is likely due to the IWV-based labels being more sensitive to the Arctic surface warming. There are however exceptions: the AR occurrence frequency from the ERA5 85th_IWV labels spreads less than those from the 85th_IVT labels in the Pacific sectors of 135 • -180 • E and 135 • -180 • W.
Figures 8 also compare the distributions of AR frequency derived from ERA5 against those from MERRA-2. Several differences between the two reanalysis products are nontrivial. The 75th, 85th, and 95th IVT-based labels from ERA5 all yielded much longer lower tails with smaller values than those from MERRA-2 in all sectors. In contrast, in figure 8(b), 75th_IVT labels from MERRA-2 resulted in longer upper tails with larger values than those from ERA5 in all sectors except for those covering East Siberian Sea (135 • -180 • E) and eastern Beaufort Sea and Canadian Arctic Archipelagos (90 • -135 • W). To a lesser degree, similar upper-tail relationships were observed between the 85th_IVT MERRA-2 and ERA5 labels in 45 • -135 • E and 45 • -90 • W. This disparity in the upper tails diminishes in the 95th_IVT labels. There is a prominent difference between the lower tails of the two reanalyses in figure 8(c). In ERA5 the 85th_IVT labels resulted in longer lower tails than the 85th_IWT, while in MERRA-2 the opposite was observed. It is noted that MERRA-2 fields were interpolated from 0.625 • × 0.5 • to ERA5's 0.25 • × 0.25 • horizontal mesh for comparison. The IWV-based labels are spatially broader than the IVT-based labels. The interpolation can not represent the fine spatial details of IVT unresolved by MERRA-2's coarser mesh. This can explain the longer lower tails of the AR occurrence frequency distributions from ERA5. Therefore, the ERA5 IVT are more precise and desirable choice than the MERRA-2 IVT for detecting and tracking Arctic ARs.
3.2. AR contributions to Arctic moisture and moisture transport 3.2.1. Arctic poleward moisture transport From 1980 to 2019, the maxima of high-latitude poleward integrated water vapor transport (IVT P ) persisted over the North Pacific and the North Atlantic, at latitudes corresponding to the mid-latitude storm tracks (cf figure B3). In recent decades, a large part of the Arctic received increasing poleward moisture transport. As seen in the dark blue areas of the decadal poleward IVT anomalies, IVT ′ P , in ERA5 (top) and MERRA-2 (bottom) in figure 9(a), the maximum positive IVT ′ P were over Bering Sea-East Siberia, Laptev-East Siberian Seas and the adjoining Arctic Ocean, Greenland-Norwegian Seas, and Baffin Bay-Davis Strait in the latest decade. Figure 9(b) shows the longitudinal transects of IVT ′ P in both reanalyses, decade-by-decade, along 60 • N (top) and 70 • N (bottom), respectively. An overall increasing trend of poleward moisture transport from the second decade was observed along both latitudes in the two reanalyses. The longitudinal transects along the two latitudes evolved roughly in sync, which is consistent with the spatial patterns in figure 9(a). The strongest positive IVT ′ P were at the lower latitudes in the Arctic, such as 60 • N in figure 9(b). However, the local maxima of positive anomalies were not just confined to where the mid-latitude storm tracks resided, suggesting that other synoptic weather systems, presumably regional Arctic ARs, contributed to the rising poleward moisture transport.
Therefore, we divided the Arctic into eight sectors to tease out the regions where the increasing poleward moisture transport occurred the most in the recent decade. By inspecting the Q-Q plots that compared the  figure B3). Therefore, we further estimated the role of ARs in transporting moisture poleward in the Arctic to test the hypothesis that the increasing transports can be attributed to ARs.

AR contributions to total poleward moisture transport
Past studies have shown that ARs contributed to a significant amount of total poleward moisture transport to the Arctic during the boreal winter. However, the estimates varied widely from 28% or 36% (Woods et al 2013), 38% (Liu and Barnes 2015), to up to 60% (Nash et al 2018), depending on the detection algorithm, estimate method, and study location. We calculated the ratio of AR-related poleward IVT, AR_IVT P , and the total poleward IVT, IVT P , with each AR label. Figure 10 shows the resultant decadal mean AR_IVT P /IVT P derived from the 85th-based labels in ERA5 and MERRA-2. The complete maps from all 75th-, 85th-, and 95th-percentile AR labels can be see in the supplementary figures S7 (IVT-based ERA5 labels), S8 (IWV-based ERA5 labels), S9 (IVT-based MERRA-2 labels), and S10 (IWV-based MERRA-2 labels). In general, the maximum ratios were consistently located over the regions from North Pacific to Bering-Chukchi Seas, from Canadian Arctic Archipelagos to southern Greenland, and the Nordic Seas across different analysis periods, reflecting the connection between mid-latitude storm tracks and Arctic ARs. These three regions of maximum AR contribution to the poleward moisture transport coincided with the three IVT ′ P hotspots in figure 9(c), deviating from the climatological poleward moisture transport (cf IVT P in decade (1980-1989, 1990-1999, 2000-2009, 2010-2019) -1989, 1990-1999, and 2000-2009  quantile values are marked with dotted red and dashed blue lines for ERA5 and MERRA-2, respectively. To be consistent, data in (b) and (c) were randomly sampled according to the sample number determined by twice the cosine value of the latitude to preserve the sample size: distance ratio.

Figure 10.
Maps showing the fraction of AR contribution to the total poleward moisture transport, AR_IVTP/IVTP, during each decade (1980-1989, 1990-1999, 2000-2009, 2010-2019) and the entire time period  identified by the 85th-percentile climate threshold applied to either IVT or IWV moisture field, with (a) 85th_IVT-and (b) 85th_IWV-based AR labels in ERA5, and (c) 85th_IVT-and (d) 85th_IWV-based labels in MERRA-2. figure B3). This suggests that ARs were indeed crucial in transporting the moisture poleward to the three hotspots. Particularly, in Arctic Archipelagos to southern Greenland (45 • -90 • W), AR's contribution to the poleward moisture transport was the most prominent feature in each decade. Consistent with the increasing poleward moisture transport anomalies over a large part of the Arctic (figure 9(a)), most grid points in the Arctic experienced increasing AR contribution to the poleward moisture transport over the decades, as seen in the deepening blue areas over time in figure 10. That means, ARs were increasingly more influential to the rising moisture transport into the Arctic, particularly to the three hotspots. Furthermore, figures 11(a) and (b) show the decadal mean fractional poleward moisture transport attributed to ARs (AR_IVT P /IVT P ) and the surface temperature (T2m) along 60 • N and 70 • N, respectively. The results derived from the 6 ERA5 AR labels (top rows) and the 6 MERRA-2 AR labels (bottom rows) during each analysis decade are largely consistent. Different from the significant poleward shift of the maximum AR occurrence frequency from land to ocean (figure 7), the highest percentages of AR contribution to poleward moisture transport persisted around Bering-Chukchi Seas (centered about 180 • ) and Canadian Arctic Archipelagos to southern Greenland (60 • W) along 60 • N; around the Nordic Seas (60 • E) and the Arctic Archipelagos to southern Greenland (60 • W) along 70 • N. And, the magnitudes of AR_IVT P /IVT P along 60 • N appear to be slightly larger than those along 70 • N.
There were striking differences in the locations of the maximum AR_IVT P /IVT P between IVT-and IWV-based labels over the Arctic Pacific sector (around 180 • in figure 11). For example, at 60 • N (figure 11(a)), the local peak ratio from the IVT-based labels appears to be larger and further east to around the Gulf of Alaska, compared against the peak ratio from the IWV-based labels. Indeed, the Gulf of Alaska was the region where the maximum climatological poleward moisture transport was observed over the Pacific (cf figure B3). This implies that IVT-based labels are better associated with mid-latitude storm-track activities and tracking more poleward moisture transport to the Arctic than the IWV-based ones. At 70 • N ( figure 11(b)) and around 180 • , the divergence of the peak ratios between IVT-and IWV-based labels diminished, and their local maxima were primarily consistent with the local maximum T2m. As shown in figure 10, and in the supplementary figures S7 and S8, the local maximum of AR_IVT P /IVT P was shifted northward to the Chukchi Sea in IWV-based labels (particularly in the last decade), implying IWV-based labels are better related to the local surface warming and moistening in the Arctic. Naturally, ARs identified by IVT-based labels accounted for more poleward IVT (figures 10 (a), (c) and S7, S9) compared with those by IWV-based labels (figures 10(b), (d) and S8, S10), with the maximum AR_IVT P /IVT P values up to 75%, 65%, and 35% for the 75th-, 85th-, and 95th-percentile labels in both reanalyses. It is noted that, among all AR contributions to poleward IVT, extreme events based on the 95th-percentile labels exhibited the greatest relative growth and contribution over time ( figure B4). The spatial pattern of the AR's proportion of total IVT (AR_IVT/IVT, in supplementary figures S11-S14) roughly matched that of AR's proportion of poleward IVT (AR_IVT P /IVT P , figures S7-S10). However, the magnitudes of AR contribution to total IVT were much smaller. The maximum ratios of total IWV attributed to ARs (AR_IWV/IWV, figures S15-S18) were in northeast Canada, southern Greenland, and western Siberia in the first two decades, then moved to Arctic Pacific, Arctic Atlantic, and Greenland in the last two decades. This is consistent with the migration of the maximum AR frequency regions seen in figures 6, B1, and S4. ARs identified by the IVT-based labels contributed more to the total IVT (figures S11-S14) than those identified by the IWV-based labels, with the maximum percentages more than 60%, 47%, and 23% for the 75th-, 85th-, and 95th-percentile labels in both reanalyses. Not surprisingly, the IWV-based labels showed a higher AR-related IWV to total IWV ratio than the IVT-based labels (figures S15-S18). It is noted that the maximum contributions of AR to total IVT and total IWV were both over Greenland, as a result of a local climatological minimum year-round.

Seasonality of AR contributions to Arctic moisture transport
Atmospheric moisture and moisture transport over the Arctic have been known to display strong seasonality (Dufour et al 2016, Nygard et al 2020, Fearon et al 2021, as indicated in the maps of three climate thresholds for IVT and IWV fields (figures A1 and A2). The global AR detection study (Guan and Waliser 2015) according to the seasonal 85th percentile of the IVT-based AR labels supplemented with the lower limit of 100 kg m −1 s −1 reveals that in the Arctic, AR frequencies were increased in the extended summer season (May-September) but decreased in the extended winter (November-March). The two extended seasons were selected based on the hydrological importance to characterize distinct dry and wet conditions in the west coast of mid-latitudes. Therefore, we further divided the AR occurrence frequency into the four seasons: Spring (March-May), Summer (June-August), Fall (September-November), and Winter (December-February) to examine their seasonal variation. Figure 12 illustrates the spatial distribution of the average annual percentage of AR occurrences during these seasons across each decade and over the 40-year study period, based on the 85th percentile IVT-based AR labels in ERA5. Results derived from the other 11 sets of AR labels are shown in figures S19-S29. Our study consistently reveals that, irrespective of the combination of the IVT or IWV field with the 75th-, 85th-, or 95th-percentile climate thresholds, a greater number of Arctic ARs are observed annually during the winter and spring seasons, with fewer ARs observed in the summer and fall seasons (figures 12 and S19-S29).
The distribution of AR occurrence exhibited distinct seasonal variations. In spring, the highest AR frequency was observed over the central Eurasian Continent, East Siberia, and Labrador Sea during the first decade (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989), shifting towards the Atlantic sector in the 1990s, while reducing over the Pacific sector (e.g. Chukchi Sea). The final decade saw a move of maximum AR occurrence towards the Chukchi and Kara Seas. In summer, the area of maximum AR occurrence shifted from continental regions in the first three decades  towards Arctic oceanic regions, specifically the Beaufort and Kara Seas, in the last decade. During fall, AR activity remained relatively low until the final decade (2010-2019), when the Eurasian continent, Beaufort Sea, and Baffin Bay saw the most ARs. In winter, the area of maximum AR occurrence moved from the Eurasian continent in the first two decades to the Arctic, Pacific, and Atlantic sectors with increased magnitude in the last decade.
In summary, the highest AR occurrence shifted from the Eurasian continent in spring and winter of the initial decade, towards the Arctic Pacific and Atlantic sectors in the same seasons by the last decade. This shift likely reflects the intensification of storm tracks in colder seasons, driving ARs further into the Arctic. Of particular note is the northward shift of maximum AR occurrence over the Eurasian continent during summer in the 1990s compared to winter, possibly indicating the seasonally shifting jet streams. In the most recent decade (2010-2019), the highest AR occurrences during summer were observed over the Arctic oceans, potentially tied to increased sea ice melting and subsequent evaporation.
In line with the findings on the average AR occurrence frequency outlined in section 3.1.2, the IWV-based AR labels (figures S20, S21, S23, S25, S27 and S29) identified a greater number of ARs overall  (1980-1989, 1990-1999, 2000-2009, or 2010-2019) and the entire time period  identified by the 85th-percentile climate threshold applied to the IVT moisture field, with 85th_IVT-based AR labels in ERA5. than the IVT-based labels (figures 12, S19, S22, S24, S26 and S28) for both reanalyses. However, the IVT-based labels emphasized ARs associated with storm tracks more distinctly. For instance, in the winter of the last decade, the IVT-based labels pinpointed maximum AR occurrences in the Arctic Pacific and Atlantic sectors where the storm tracks resided. Conversely, when employing the IWV-based labels, regions of maximum AR occurrence showed a larger magnitude and expanded across the entire Arctic Ocean.
The two reanalyses showed qualitatively similar seasonal variability during the 40 years (figures 12, 13, and supplementary figures S19-S40). There was a tendency for more fractional poleward IVT attributed to ARs during the fall, winter, and spring but less so in the summer (figures 13 and S30-S40). In the winter, the local maxima of AR contribution were pronounced over the Arctic Pacific and northeast Canada-southern Greenland. Yet, a smaller fraction of poleward IVT over Siberia and the Laptev-Kara Seas can be attributed to ARs. This distinct pattern of AR contribution to poleward IVT remained the same but with a weaker magnitude in the spring. Finally, each Arctic location experienced a smaller AR contribution to poleward IVT in the summer. Throughout the records, ARs' poleward moisture transport was the largest over the Arctic Pacific in the winter during 1980-1989 and 2010-2019. The northeast Canada-southern Greenland also exhibited local maxima of poleward IVT attributed to ARs in the winter and spring of 1990-1999 and 2000-2009.

AR trends in climate 3.3.1. Spatial patterns of AR linear trend
We further examined the long-term linear trends of AR occurrence detected by each Arctic AR label. There is a caveat, however. As studies have shown that the Arctic amplification is the most prominent in recent decades Francis 2006, Serreze et al 2009), calculation of climate trends of moisture or moisture transport in the Arctic has been known to sensitively depend on the definition in the temporal domain (Dufour et (1980-1989, 1990-1999, 2000-2009, or 2010-2019) and the entire time period  identified by the 85th-percentile climate threshold applied to the IVT moisture field, with 85th_IVT-based AR labels in ERA5.
years (e.g. figures 6, 7, B1 and S4). Overall, the two reanalyses yielded comparable spatial patterns derived from the IVT-or IWV-based labels. Consistent with the magnitude of the AR frequency (figures 6, B1 and S4), increasing the climate threshold from the 75th to the 95th percentiles decreased the amplitude of AR trends. Still, scaled by the individual climatology of 30-or 40-year mean, the extreme events represented by the 95th-labels exhibited the largest increase of occurrence (figures S43-S46).
Figures 14 ab show the 1990-2019 AR occurrence trends estimated from the 85th_IVT-and 85th_IWV-based AR labels in ERA5. In figure 14(a), the 85th_IVT-based AR single-sign trends within 95% confidence intervals were all positive and extended from the Chukchi Sea, across the Arctic Ocean, to northeast Greenland. Similarly but with stronger signals, in figure 14(b), the 85th_IWV-based ARs occurred increasingly more frequently over most of the Arctic Ocean, extending to the Greenland Sea and Barents-Kara Seas. The strongest trending regions were over central Arctic Ocean margins, such as the East Siberian Sea and the Laptev Sea. The intensified trends of AR occurrence in the recent 30 years in figures 14(a), (b), B5 and B6 (cf 40-year trends in figures S41 and S42) were mostly located over the Chukchi Sea for the IVT-based labels, and the East Siberian-Laptev Sea and Barents-Kara Seas for the IWV-based labels.

Trend analysis by scaling decomposition
It is possible to infer the possible mechanisms leading to the different AR occurrence trends between the IVT-and IWV-based labels. As an example, we decomposed the AR trends obtained by the 85th_IVT-based labels ( figure 14(a)) into the long-term trends, most likely directly associated with Arctic surface warming, and the 'remainder' trends, where the changes in atmospheric dynamics more likely manifest. We used a scaling approach similar to that used in Ma et al (2020) to preserve the non-negative property of IVT while detecting ARs. At each 3-hourly time step, for each unique longitude, latitude, and pressure-level coordinates, a factor of q(p) c /q(p) m was applied to scale the specific humidity, q(p), in equation (3). Here, q(p) c is the 40-year climate mean specific humidity at each location. q(p) m is the annual mean specific humidity at each location for any year between 1980-2019. The approach is intended to suppress the long-term variability and trends of the specific humidity that roughly scaled with the variability and trends  (g) are the correlation coefficients between the T2m and the AR occurrence frequency derived from the 85th_IVT-and 85th_IWV-based AR labels in ERA5, respectively, during 1990ERA5, respectively, during -2019 shows The correlation coefficients between AR occurrence frequency using the 85th_IVT-based labels and the AR occurrence frequency using the 85th_IWV-based label. The grey dots are plotted over regions of values outside of the 95% confidence intervals for trends (a)-(e) or correlation coefficients (f)-(h).
of surface temperature as inferred by the Clausius Clapeyron relation. Then, we further applied the AR detection and tracking method described in section 2.2 with the re-calculated IVT using the scaled moisture and the corresponding scaled-IVT 85th percentile monthly climate threshold. At each location, the difference between the original 'total' AR occurrence trend and the scaled trend represents the long-term climatological effects of surface warming. In (Ma et al 2020), the scaled trend was considered as the dynamically driven component, while the long-term trend was driven by the thermodynamic changes.
The scaled trends displayed a narrow region of increased AR occurrence confined to the Chukchi Sea on the Arctic Pacific side ( figure 14(c)). This was in the same region with the conspicuous increase of AR occurrence in the latest decade in figure 8. The results imply that changes in atmospheric circulation, such as the Arctic cyclones or teleconnection patterns, that operated in time scales deviating from that of the long-term surface warming might be the driver of increasing AR occurrence over the region.
In contrast, the scaled trends also suggested decreased AR occurrence over East Canada-northern Hudson Bay, Norwegian Sea, and East Siberia. However, the long-term increasing trends of the atmospheric moisture field largely contributed to increases in total AR trends over almost the entire Arctic region ( figure 14(d)). The larger trends attributed to the long-term atmospheric moisture were over the central Arctic Ocean margins, from the Pacific sector, across the East Siberian Sea, Laptev Sea, Barents Sea, and Kara Seas, towards the Greenland Sea-Baffin Bay, which reflects the long-term increase of atmospheric moisture scaled with increasing T2m (figure 14(e)) following the Clausius-Clapeyron relationship. Furthermore, the larger long-term trending area roughly coincided with the intensified AR trends obtained by the 85th_IWV-based labels in figure 14(b). The coincidence of the increased AR trending area in figure 14(b) and (d) further supported that AR occurrence obtained by the IWV-based labels, especially with the 75thand 85th-percentile climate thresholds, was more sensitive to the Arctic increasing moistening and surface warming previously discussed. AR occurrence estimated using the IVT-based labels was more easily affected by changes in the moisture transport mechanisms, which leads to localized AR trends that stood above the long-term Arctic warming and moistening trends.

Relationship with surface warming
As discussed above in figures 14(a), (b), B5 and B6, the large trends of AR occurrence over the Chukchi Sea from the IVT-based labels, as well as those in the East Siberian-Laptev Sea and Barents-Kara Seas from the IWV-based labels, coincided with the strong warming areas in T2m (figures 14(e), S47(b) and (d)), and were consistent between the two reanalyses (figures B5, B6, S41, S42 and S47). To quantify the association between AR and T2m, figures 14(f) and (g) show the spatial distribution of linear correlation coefficients in the recent 30 years between T2m and AR occurrence frequency identified by the 85th_IVT-and the 85th_IWV-based label, respectively. For both labels, the highest positive correlation values existed over Greenland, a region with maximum AR occurrence (figure 6) and highest fraction of AR contribution to poleward moisture transport (figure 10). As stated previously, Greenland had the minimum climatological mean IWV (figures A1(b) and A2(b)) and poleward moisture transport (figure B3) but received the maximum positive poleward moisture transport anomalies in the most recent decade (figure 9). The highest positive correlation region in Greenland implies that Greenland ARs were more sensitive to the surface warming effects than other Arctic regions. As seen in figure 14(e), Greenland was not the strongest T2m warming area, yet a change in T2m could influence the local moisture content and the subsequent AR occurrence that further transported the moisture poleward.
For both labels, the second highest positive correlation regions were over the Atlantic sector, extending to the Barents-Kara Sea and the Pacific sector (i.e. Bering Land Bridge-Chukchi Sea), which coincided with the strong warming areas in T2m (figure 14(c)). Over the Atlantic sector, areas with the highest correlation collocated with the largest AR occurrence trends from the 85th_IWV-based labels ( figure 14(b)). Over the Pacific sector (i.e. Chukchi Sea), the highest correlation area coexisted with the largest AR trends using both AR labels (figures 14(a) and (b)). In addition, with the 85th_IWV-based label, the East Siberian Sea experienced the most prominent AR trends ( figure 14(b)) and amplified surface warming ( figure 14(e)). However, the linear relationship between the AR occurrence and T2m was not distinct over the East Siberian Sea ( figure 14(g)), implying that other factors besides surface warming might contribute to higher AR occurrence trends.
Overall, the surface temperature positively correlated with the AR occurrence frequency over the entire Arctic for all 12 AR labels in the recent 30 years and across the entire 40 years (figures S48-S53). However, their relationships were not distinctly strong. Even over Greenland, the largest linear correlation coefficients were about 0.4 for the 75th percentile-based AR labels, which means that T2m can only explain 16% of the variation in AR occurrence in linear regression. As expected, AR frequency detected by IWV was more correlated with T2m than that detected using the IVT-based labels. Moreover, the relationship between T2m and AR occurrence frequency became weaker with increasingly restrictive climate thresholds from the 75th to the 95th percentiles in both reanalyses. Regardless, the relationships between the T2m and the AR frequency do not appear to be strikingly different between the recent 30 years (figures S48-S53(d) and (e)) and across the entire 40 years (figures S48-S53(a) and (b)). Lastly, as expected, AR occurrence frequencies identified by the IVT-and the IWV-based labels were highly correlated for all three climate thresholds during 1980-2019 (figures S48-S53(c)) and 1990-2019 (figures 14(c) and S48-S53(f)). The maximum coefficients were greater than 0.7 over the Greenland-Baffin Bay, central Siberia, and the Bering-Chukchi Seas.

Relationship with teleconnections
Existing research has demonstrated significant modulation of ARs by climate variability in the ocean and the atmosphere (e.g. Guan et al 2013, Lavers and Villarini 2013, Guan and Waliser 2015, Ma et al 2020, Sharma and Déry 2020. The correlation between Arctic ARs and climate variability modes, however, remains largely unexplored.  suggested through ensemble CESM2 model experiments that Arctic ARs in the Barents-Kara Seas, and their surroundings, increased from 1979-2021 due to anthropogenic warming and teleconnection with internal climate variability in the tropical Pacific. This led to surface melting that delayed winter-time sea-ice recovery in recent decades. Hence, the conjecture regarding internal climate variability's influence on Arctic ARs can be inferred from previous studies on sea ice variability (e.g. Hurrell 1995, Hilmer and Jung 2000, Rigor et al 2002, L'Heureux et al 2008, Kwok et al 2013, Hegyi and Taylor 2017, Caian et al 2018. (Zhang et al 2020), utilizing quantile regression, observed that teleconnection patterns (AO, NAO, and PNA) exerted a more profound influence on Arctic sea ice than internal climate variability (such as ENSO, AMO, and PDO the Pacific Decadal Oscillation). However, the positive phase of AMO and PDO augmented the effects of these teleconnection patterns. It is plausible that AR modulation by internal climate variability could mirror these spatial-temporal patterns, albeit compounded by interactions with anthropogenic warming. For instance, (Nakamura et al 2015) discovered that sea-ice reduction in the Barents region triggered a Rossby wave response, inducing an anomalous circulation pattern akin to the negative phase of AO/NAO throughout most of the 1980s. Furthermore, (Hamouda et al 2021) suggested in their future projection that AO and NAO would respond differently to a warming climate, effectively decoupling the two. AO variability would lessen over the Atlantic, strengthen over the Pacific, while the NAO pattern would remain relatively unchanged.
Here, we computed for each set of AR labels the anomalies of Arctic AR occurrence in each decade by subtracting the decadal mean at each grid point. We defined the positive or negative phase of AO, NAO, or PNA with index values > 0.5 or < −0.5. Then, we composited the spatial distributions of anomalous AR occurrence for each phase. Figure 15 shows the results with the 85th_IVT-based AR labels from ERA5, while the remaining results are in supplementary figures S54-S71. It is noted that the extreme ARs based on the 95th percentile thresholds were the most influenced by the teleconnections (figures S74-S85).
In order to elucidate the decadal difference of the composited spatial patterns in figure 15, we have included time series of AO, NAO, PNA, ENSO, IPO, and AMO indices in figure B7. A notable shift in internal climate variability is AMO's transition from negative to positive around 2000, paralleled by a shift from predominantly positive ENSO and IPO phases in the first two decades to more negative phases in the subsequent two. Nonetheless, as further examined in the subsequent subsections, these transitions do not seem to fully explain the decadal changes of AR occurrence conditioned on the phases of AO, NAO, and PNA (figures 15 and S75-S85). This is similar to the previous research of sea ice variability. It is however possible that in the last decade the anthropogenic warming was caught up and amplified by the positive AMO phase, resulting in a significant and distinct increase of AR activities over the ocean.

Arctic oscillation (AO)
During the positive phase of AO, the jet stream is shifted farther north, and so are the storm tracks (Thompson and Wallace 1998). Under these conditions, all AR labels consistently displayed that Arctic ARs occurred more frequently over a large area from the Norwegian Sea across the central north Eurasian continent; however, ARs decreased over a relatively smaller area including the entire Greenland and the Baffin Bay/Davis Strait (figures 15(a) and S54-S59(a), (c)). Moreover, the magnitudes of AR anomalous frequency during the positive AO were the strongest with the 95th_IVT-based labels (figures S58 and S59), which were about the same as the climatology (i.e. the decadal means in figures B1(a) and (c)) while those derived with the 75th-and 85th-labels (figures S54-S57) were more than half of the climatology (cf figures 6(a), (c) and S4(a), (c)). During the negative AO, AR frequency anomalies were mostly reversed (figures 15(b) and S54-S59(b), (d)).
There are distinct features of the AR frequency anomalies in each decade. Most notably, in the most recent decade (2010-2019) the anomalous patterns developed a significant branch between the North Pole and Beaufort Sea (figures 15(a) and (b)). Moreover, the occurrence of AO phase skewed toward negative in the first decade but migrated to be more positive in early 90s before shifting back toward symmetry after 2000, as shown in figure B7 and discussed in (Nakamura et al 2015). The positive AR anomalies over the northeast Canada and Greenland during the negative AO (figures 15(b) and S54-S59(b), (d)) were stronger during 1980-1999 than during 2000-2019, and attained stronger magnitudes than negative anomalies during the positive phase (figures 15(a) and S54-S59(a), (c)) in the first two decades.

North Atlantic Oscillation (NAO)
The distributions of anomalous AR occurrence modulated by NAO in figures 15(c), (d) and S60-S65 were similar to those modulated by AO (figures 15(a), (b) and S54-S59), but with enhanced signals over the Arctic-North Atlantic Oceans. In the positive phase of NAO ( figure 15(c)), the westerly jet and storm track associated with the strong Icelandic Low and the Azores High moved northward (Hurrell et al 2003), setting up the ideal condition for ARs to proceed across the Atlantic sector into the Arctic. Consequently, the area of positive anomalous AR occurrence (shaded in blue in figure 15(c)) moved to the Greenland, Barents, and Kara Seas as opposed to the Eurasian continent during the positive AO ( figure 15(a)). Meanwhile, anomalous negative AR occurrence was observed over northeast Canada and Greenland, which was similar to effect of the positive AO modulation. During the negative NAO (figure 15(d)), the Greenland blocking regime (Cheng andWallace 1993, Kimoto andGhil 1993) diverted the jet stream and enhanced AR occurrence was observed over the northeast Canada and Greenland while reduced occurrence was detected over the Greenland, Barents, and Kara Seas. This seesaw pattern was consistent across all the 12 AR labels (figures S60-S65).
Similar to AR anomalies related to AO in 2010-2019, those related to NAO had a significant polar extension; however, the pattern bridged from the Atlantic to the Chukchi or Bering Seas in the Pacific (figures 15(c) and (d)), implying a decoupling between NAO and AO in the warming Arctic but different from that in (Hamouda et al 2021). The magnitudes of AR occurrence anomalies related to NAO were stronger than those related to AO. Especially, with the 95th-labels (figures S64 and S65), the former were more than twice of the latter over northeast Canada and Greenland. Over time, the NAO phases underwent similar shifts as AO phases. Regardless, the magnitudes of positive AR anomalies over northeast Canada-Greenland during the negative NAO (in dark blue in figure 15(d)) appeared stronger than those of the negative anomalies in the opposite phase (figure 15(c)) without clear trends.

Pacific North America (PNA)
During the positive phase of PNA, two positive geopotential height anomalies were centered over Hawaii and northwestern Canada, and two negative anomalies were over the Aleutian Islands and the southeastern US Gutzler 1981, Feldstein 2000). The associated jet stream increased, shifted eastward, and exited towards North America, which favored AR moisture transport along the eastern Pacific to the Arctic. Consequently, ARs occurred more frequently over Canada but less so over Eastern Siberia ( figure 15(e)). Central Siberia also experienced increased ARs during each analysis period. During the negative PNA, the anomalous height fields over the four centers were reversed, so were the related AR occurrence anomalies. Instead, AR entered the Arctic through western North Pacific and Bering Sea (figure 15(f)). Moreover, northwestern Canada, Beaufort and Chukchi Seas, and Eastern Siberia received increasingly more AR occurrence over the four decades. That is, the strength of positive AR occurrence anomalies generally exceeded those of negative anomalies in these regions (figures 15(e) and (f)).
Note that the AR frequency in the eastern Beaufort Sea had clearly been positively affected by the negative PNA (figure 15(f)), and in the most recent decade, by the positive AO and NAO (figures 15(a) and (c)). Similarly, the AR frequency over the Arctic Archipelago had been positively affected by the negative NAO, even though the surface warming in this area was weak to negative in the reanalyses ( figure S47). Moreover, poleward moisture transport increased in these areas, especially in the latest decade (figures 9(a) and S7-S10), along with increased AR contributions to the transport most evident in the colder seasons with both IVT-and IWV-based labels in both reanalyses (figures S19-S40). Since teleconnections are the most active in the boreal winter, future investigations on the time-dependent mechanisms associated with the changing teleconnections in the Arctic are imperative to clarify the relationship between the physical facets of ARs and Arctic sea ice decline.

Conclusions
Using 3-hourly ERA5 and the MERRA-2 data for ARTMIP from January 1980 to December 2019, a total of 12 Arctic AR labels were created based on combinatory conditions of either the IVT or the IWV field exceeding one of the monthly climate thresholds (75th, 85th, and 95th percentiles) at each location in the Arctic. An AR event is strung by extracted objects satisfying the geometric constraint (length ⩾ 1500 km and length-to-width ratio ⩾ 2), the spatial proximity and morphological similarity between adjacent time steps, and persisted at least 18 hours altogether. Consistent across all AR labels, there was an increasing trend of AR events with a striking quasi-biennial oscillatory signal over the 40 years. The increasing trend of AR event counts was the most prominent in the recent decade. The quasi-biennial oscillatory signal might be due to the influence of QBO on AR activities through AO, and deserves further study.
Further spatial exploratory analysis of all 12 labels consistently showed that the locations of maximum Arctic AR occurrence frequency shifted poleward, from overland in 1980-1999 to over the Arctic Ocean and its outlying Seas in 2000-2019, as regions across the North Atlantic, the Arctic, to the North Pacific Oceans trending toward higher AR occurrence, warmer surface temperature, and more column-integrated moisture. In the meantime, the AR contribution to the increasing poleward IVT also increased, estimated to reach up to 75%, 65%, and 35% with the 75th-, 85th-, and 95th-percentile labels in both reanalyses. The maximum ratios were located over the regions connected with the North Pacific and North Atlantic storm tracks.
The long-term trends of heightened AR frequencies over the ocean coincided with areas of warming hot spots. Using a scaling approach similar to that in Ma et al (2020), we decomposed the AR trend obtained by the 85th_IVT-based labels into the long-term trend and the scaled remainder trend. The former was more likely directly associated with Arctic surface warming, and the latter was more likely to show the effects of atmospheric circulation changes such as Arctic cyclones or teleconnection patterns. Even though the increase of Arctic AR occurrence was primarily associated with long-term Arctic surface warming and moistening, the decomposition also showed positive scaled trends confined to the Chukchi Sea ( figure 14(c)). The AR occurrence was indeed strongly modulated by the changing teleconnections including AO, NAO, and PNA, with notable anomalies in the Arctic-Pacific sector during the latest decade.
The paper's main data visualization features the 85th-percentile-based labels, while the 75th-and 95th-percentile labels are shown in the supportive information. This is because we expect the 85th-percentile labels will be the most commonly adopted in future research. The 75th-percentile labels likely could capture the broad presence and accumulation of surface hydrometeorological effects as we found in Zhang et al (2021). Moreover, the 75th IVT-and IWV-based labels, as well as the 85th IWV-based labels, are likely more sensitive to Arctic surface warming and moistening. Therefore, these labels detected some broadening of AR-affected areas over the study period, while the rest of the labels did not. It was also noted that the ERA5 IVT field was more precise than the coarser-resolution MERRA-2 IVT for AR detection and tracking.
In the interests of studying extreme events, however, the 95th-percentile labels exhibited the largest increase rates of AR occurrence and AR contribution to poleward IVT, especially from 2000-2009 to 2010-2019 in ERA5. When scaled by the corresponding climatology, the occurrence of extreme events identified by the 95th-percentile labels was the most influenced by the teleconnection patterns and exhibited the greatest relative trends in both reanalyses.
In conclusion, we applied a data-driven approach using two moisture fields, three levels of climate thresholds, and two reanalysis data sets to create an ensemble of 12 AR labels to present an atlas of Arctic AR climatology. However, the choice of Arctic AR detection labels should be adaptive to the problem to be addressed. For example, questions like which set of AR labels is most in synergy with the extratropical storms, which set of labels captures the cloud-radiative effects the best, or which set of labels is most predictive of surface energy balance changes remain open and should be explored in the future.

Data availability statement
The data that support the findings of this study will be openly available following an embargo at the following URL: https://purr.purdue.edu/publications/4322/1 (Tung et al 2023).

Plain language summary
Atmospheric Rivers (ARs), the narrow filaments of enhanced moisture transport in the lower troposphere, are conspicuous pathways for poleward moisture transport. We hypothesize that the changes in Arctic AR's occurrence frequency, strength, and location are integral to Arctic warming. We approached the data analysis in the Arctic, a data-sparse region compared with the mid-latitudes, with a statistical design to identify and track Arctic ARs based on the climate significant-event or extreme-event criteria. We created twelve sets of AR labels for 1980-2019 using the ERA5 and MERRA-2 reanalysis data and presented the Arctic AR climatology. All labels consistently showed that there had been an increasing number of Arctic AR events in recent decades. The labels more sensitive to Arctic surface warming and moistening could also detect some broadening of the areas affected by AR events. The locations of maximum Arctic AR frequency of occurrence shifted poleward from the Eurasian and American continents in 1980-1999 to over the oceanic areas in 2000-2019, with regions across the Atlantic, the Arctic, to the Pacific Oceans trending higher AR occurrence, warmer surface temperature, and more air moisture. These increased AR activities coincided with surface-warming hot spots but not entirely with the September sea ice retreat. ARs were increasingly responsible for the rising moisture transport into the Arctic. Even though the increase in Arctic AR activities was primarily associated with long-term Arctic surface warming and moistening, changes in atmospheric circulation could also affect AR occurrence locally, such as over the Chukchi Sea. The changing teleconnection patterns, such as Arctic Oscillation, North Atlantic Oscillation, and Pacific/North America patterns, strongly modulated AR activities, especially their extreme events. The Arctic AR labels and the detailed graphics can help navigate the uncertainty of detecting and quantifying Arctic ARs and their associated effects in current and future studies.

Appendix A. AR detection and tracking algorithm
The creation of an ensemble of twelve Arctic AR labels-using factors including two moisture fields (IVT and IWV) and three relative climate thresholds (75th, 85th, and 95th percentiles)-with ERA5 and the MERRA-2 source data for ARTMIP from January 1980 to December 2019 was based on the detection and tracking algorithm of (Zhang et al 2021) modified with a Lagrangian AR tracking routine such as the one used in (Guan and Waliser 2019). Step 1: AR object selection based on climate thresholds The first step was to create monthly climate threshold values for the moisture fields. For each moisture field (IVT or IWV) to the north of 35 • N, we initially aggregated at each grid point the datum at 1200 UTC of each day during neutral or weak El Niño-Southern Oscillation events. Then, we calculated their monthly 75th, 85th and 95th percentiles as the climate thresholds (figures A1 and A2). In addition to the common 85th percentile climate thresholds used for AR detection in the mid-latitude, the 95th percentile climate thresholds represent the extreme AR events in the Arctic, while the 75th percentile thresholds can broadly capture AR effects. The choice of the neutral or weak ENSO events was based on the bi-monthly NOAA Multivariate ENSO index (MEI.v2, (e.g. Winton 2006)) whose values were within ±1, consistent with the criteria used previously in Zhang et al (2021).
Figures A1 and A2 plot these three levels of monthly climate thresholds of (a) IVT and (b) IWV fields in the Arctic (above 60 • N) in January and July based on ERA5 and MERRA-2, respectively. Here, January is representative of the boreal winter, and July is representative of the summer as a way of summarizing the threshold instead of listing all twelve months. Overall, the spatial characteristics of the IVT and IWV climate thresholds in January and July were consistent between ERA5 and MERRA2. Figures A1 and A2 show that from the 75th to the 95th percentiles, the threshold value at each grid point elevated and became successively more restrictive to AR detection. In agreement with the Clausius-Clapeyron relation, IVT and IWV thresholds were generally stronger in the summer than in the winter. The minima of IVT and IWV were located over Greenland and more prominent in the summer. For both reanalyses, two maximum IVT regions were evident over the Arctic Pacific (Bering Sea) and Arctic Atlantic Ocean (extending from the Labrador Sea to the Norwegian Sea) in both January and July, where the end of the two storm tracks resided (figures A1(a) and A2(a)). While from the winter to the summer, the maximum IWV switched from the Pacific and Atlantic Oceans to the Eurasian and North American continents (figures A1(b) and A2(b)). This might be due to the different heat capacities between ocean and land, which leads to more evaporation over the continent in the summer (Vázquez et al 2016) but less in the winter. Figure A3 shows the differences between ERA5 and MERRA-2 in terms of the three levels of climate thresholds for IVT and IWV. Basically, in the IVT fields, the three climate thresholds calculated in ERA5 were smaller over the American and Eurasian continents but larger over a fractional part of the Arctic Ocean than Figure A3. Difference of the three levels of climate thresholds between ERA5 and MEERA-2 (ERA5−MERRA2) for IVT (a, in kg m −1 s −1 ) and IWV (b, in kg m −2 ) in the Arctic for the months of January and July. those calculated in MERRA-2 ( figure A3(a)). For the climate threshold values in the IWV fields, ERA5 exhibited smaller overall values, especially in the two continents during summer, than those in MERRA-2. Although in winter, slightly larger values were observed over the Arctic Ocean for the 75th and 85th percentile climate thresholds.
It is noted that figures A1(a) and A2(a) outline the fixed lower limit of the IVT field (100 kg m −1 s −1 ) used by Waliser (2015, 2019)and (Fearon et al 2021) in the polar areas to supplement the seasonal percentile thresholds in their AR detection method. In January, the majority of our monthly climate threshold values over the Arctic (excluding Arctic Atlantic and Bering Sea) were below 100 kg m −1 s −1 . Therefore, figure A1 shows that northward of 60 • N, the fixed lower limit could result in fewer detected ARs from Eurasia, across the entire Arctic Ocean, to North America. Furthermore, as the thresholds decreased from the 95th to the 75th percentiles, this under-detected region enlarged to almost the entire Arctic. In July, the under-detected regions shrunk significantly, but our monthly climate thresholds still permitted more frequent ARs to be detected over or around Greenland. Overall, for both reanalyses, our monthly percentile threshold yielded more detected ARs throughout the entire Arctic, especially during the winter months and over Greenland in the summer, compared with the previous results of the global AR detection in the Arctic (Guan and Waliser 2015).
Next, at each time step, we identified the grid points whose moisture values equaled or exceeded a threshold. The principal curves method (Hastie and Stuetzle 1989) was then applied to determine the length of the curvy patterns formed by aggregating the maximum IVT or IWV values at each latitude and longitude within the identified spatial pattern. The algorithm had the periodic boundary condition about each latitude. The width of the pattern was calculated as the surface area divided by the length. A pattern with a length ⩾ 1500 km and a length-to-width ratio ⩾ 2 was kept as a potential AR object.
Step 2: AR event tracking Then, we examined the time sequence of identified AR objects to construct AR events. We defined an AR event as a run of AR objects that persisted for at least 18 hours in the Arctic region. We modified the Eulerian AR event tracking algorithm suited for regional studies in the United States in Zhang et al (2021) with the Lagrangian tracking framework based on the spatial proximity and morphological similarity in Guan and Waliser (2019).
In this Lagrangian framework, for a given time step, an AR object was deemed to be a continuation of an AR object from the previous time step as long as they overlapped spatially, i.e. possessing spatial proximity. In the cases of AR object separation or merger, when multiple AR objects were continued from one object or vice versa, only one object was selected. In (Guan and Waliser 2019), this selection was based on the maximum morphological similarity between AR shapes at two consecutive time steps, determined by the Jaccard similarity index (J; (Deza and Deza 2014)) between the area-weighted IVT values of two identified AR objects as follows: with i a grid point in the analysis domain, A i the mesh area encasing only the ith grid point, the superscripts (1) and (2) notating the two AR objects from consecutive time steps. The point i may be inside an area affected by AR (1) , AR (2) , both ARs, or neither ARs. As seen in the equation (A1), J IVT is zero if the two AR objects do not intersect in space; it is 1 if the two AR objects completely overlap in space and have identical IVT values. We calculated the Jaccard similarity index for both IVT-and IWV-based labels. For the latter, the formula was: As in Guan and Waliser (2019), after a selection based on the maximum J took place, the rejected AR objects were considered as current genesis of new AR events in the case of separation or terminations of previous events in the case of merger. We performed the event tracking using the IWV and IVT fields that passed the thresholds at all horizontal grid points north of 35 • N in the reanalyses. At the same time, we defined Arctic AR events as those penetrating beyond 60 • N and persisted for at least 18 hours in the Arctic. The outcomes were 12 Arctic AR labels expressing spatial and temporal AR information based on the factors of moisture fields and climate thresholds from January 1980 to December 2019 from the ERA5 and MERRA-2 reanalyses. The detection and tracking algorithm was implemented via distributed-parallel computing, specifically, the divide-and-recombine approach using the R-based DeltaRho software back-ended by a Hadoop system (Cleveland andHafen 2014, Tung et al 2018). Figure B2. Spatial distribution of relative percentage increase (in %) of mean AR occurrence frequency during each decade compared against the previous decade (labeled on each row, for example, 'D2−D1' denotes the percentage increase in the second decade of 1990-1999 compared against the first decade of [1980][1981][1982][1983][1984][1985][1986][1987][1988][1989] according to the ARs identified by the 75th-, 85th-, and 95th-percentile climate thresholds applied to the IVT and IWV moisture fields in ERA5. Figure B3. Spatial distribution of the decadal mean poleward moisture transport averaged over each decade (1980-1989, 1990-1999, 2000-2009, 2010-2019) and the entire 40 years (IVTP in kg m −1 s −1 ), using (a) ERA5 and (b) MERRA-2. Figure B4. Spatial distribution of relative percentage increase (in %) of the fraction of AR contribution to the total poleward moisture transport (total IVT) during each decade compared against the previous decade (labeled on each row, for example, 'D2−D1' denotes the percentage increase in the second decade of 1990-1999 compared against the first decade of [1980][1981][1982][1983][1984][1985][1986][1987][1988][1989] according to the ARs identified by the 75th-, 85th-, and 95th-percentile climate thresholds applied to the IVT and IWV moisture fields in ERA5.