Atmospheric River Variability Over the Last Millennium Driven by Annular Modes

Atmospheric rivers (ARs) are filamentary conduits of intense moisture transport crucial for water delivery to mid‐latitude coastal regions. How ARs have responded to extratropical climate variability remains poorly understood despite ARs being features of the extratropical atmosphere. Here, using “Last Millennium” simulations, we characterize the role of annular modes of extratropical variability on ARs and moisture transport. We find that positive (negative) phases of the annular modes intensify (weaken) and weaken (intensify) ARs over the subpolar and subtropical latitudes, respectively, with up to ∼20–25 mm/month associated changes in precipitation. Importantly, the annular modes comprise the primary mode of AR variability over the last ∼1,000 years. We also separately examine the annular modes' influence on storm track activity and find it distinct from that on ARs and moisture transport, despite the storm tracks being associated with ARs and overlapping with strong moisture transport. Lastly, our results provide a robust paleoclimate baseline from which to contextualize projected 21st century AR intensification.


of 14
The relationship between ARs and extratropical variability, however, requires further understanding, despite ARs being features of the extratropical atmosphere. Although several studies have examined the regional role of (a) northern extratropical variability on ARs over Europe and North America (e.g., Brands et al., 2017;Guan et al., 2013;Lee et al., 2022;Zavadoff & Kirtman, 2020) and (b) southern extratropical variability on ARs over Australia and Antarctica (e.g., Meneghini et al., 2007;Pohl et al., 2021;Reid et al., 2022;Wille et al., 2021), respectively, the global relationship between ARs and extratropical variability is not well known. In particular, the global role of the dominant modes of extratropical variability-known as annular modes due to their zonal symmetry-on ARs has not been investigated. The annular modes account for 35% or more of variance in anomalous extratropical zonal winds (i.e., deviations in the jet stream), with positive and negative phases of the modes reflecting poleward and equatorward movement of the jet streams, respectively (Gerber & Thompson, 2017;Thompson & Wallace, 2000). These north-south displacements are of primary importance for the location of extratropical waves and associated extratropical cyclones that contribute to ARs (the so-called storm tracks; Shaw et al., 2016). Yet the relationship between the annular modes and the intensity of the storm tracks nevertheless remains ambiguous, with observational studies indicating that the two are largely decoupled (Simmons & Hoskins, 1978;Thompson & Li, 2015;Thompson & Woodworth, 2014).
Characterizing AR variability remains a challenge in part because satellite measurement and data assimilation systems required to directly observe ARs have only become widely available since roughly 1980 (Blamey et al., 2018;Gelaro et al., 2017;Lavers et al., 2012;Neiman, Ralph, Wick, Lundquist, & Dettinger, 2008). Climate model simulations, which do not depend on observations, can circumvent these sampling challenges by extending assessments of ARs beyond the past few decades (e.g., Lora et al., 2017), though model performance and biases must be considered. Here, we explore how the annular modes affect ARs, and moisture transport more broadly, using the "Last Millennium Ensemble" Project (LME; Kay et al., 2015;Otto-Bliesner et al., 2016). The LME simulates five million days of AR activity over the last millennium across 12 experiment members, with each member subject to identical radiative forcing changes but with a different realization of natural variability. The LME not only provides robust AR sampling but chronicles multiple realizations of internal climate variability across the full range of radiative forcings (both natural and anthropogenic) that humans have experienced over the past millennium. This is not possible, for instance, with a pre-industrial control experiment, making the LME an ideal tool from which to investigate how the annular modes affect global moisture transport and to place potential modern changes in a historical context.

Data and Methods
We employ 12 members of the LME conducted with Community Earth System Model 1 (CESM1), a fully coupled global climate model that successfully simulates AR features (e.g., Baek & Lora, 2021;Shields et al., 2018;Skinner et al., 2020). The LME forces CESM1 with reconstructed estimates of natural solar, volcanic, and orbital radiative forcing, as well as pre-industrial and historical radiative forcing from greenhouse gases, aerosols, and land use/land change. In addition to these "Full forcing" LME experiments, we examine five members of the parallel "Volcanic only" experiment, which includes volcanic forcing but excludes all other natural and anthropogenic forcings. Note that we use "LME" throughout our study to refer to the "Full forcing" experiment (and not the "Volcanic only" experiment). Outputs from the LME and "Volcanic only" experiments span 850-2005 and are available on daily resolved, ∼1.9° latitude by 2.5° longitude grid resolution.
To calculate the annular modes, we remove (a) 10-year running averages from longitudinally averaged zonal daily winds at every latitude at 850, 500, and 200 hPa (the three vertical levels available), and then (b) the seasonal cycle from the detrended data. We perform Empirical Orthogonal Functions (EOFs) analysis on the resultant wind field anomalies (i.e., wind anomalies in a [daily time by three vertical levels] array), multiplied by the square root of the cosine of latitude, defining Northern Annular Mode (NAM) and Southern Annular Mode (SAM) as the normalized principal component of the leading pattern of variability over 20°-70°N and 20°-70°S, respectively (Hartmann & Lo, 1998;Lorenz & Hartmann, 2001;Thompson & Woodworth, 2014). These first EOFs are well-separated from subsequent patterns of variability (North et al., 1982). We define the signs of the annular modes such that positive phases of NAM and SAM correspond to the jets shifting poleward. We examine the influences of the two annular modes on ARs (see below for description of the detection algorithm), integrated vapor transport (IVT), zonal winds, and eddy kinetic energy (EKE), removing the seasonal cycle from each respective field beforehand. Note that we remove seasonal influences from these fields in addition to also removing the seasonal cycle from the annular modes, as the distribution of extrema in the fields can be seasonally dependent. To identify annular mode-related changes, we (a) regress these fields onto NAM/SAM principal components and (b) composite these fields over positive and negative phases of the annular modes, which we define as when the principal components are one standard deviation (SD) above and below (e.g., Ambaum et al., 2001;Doddridge & Marshall, 2017) their 850-2005 means, respectively. For regressions, significance (p < 0.01) is determined using a two-sided Wald Test (null hypothesis is that slope equals zero). For composites, significance (p < 0.01) is determined using a two-sided t-test on the respective fields (a) during positive or negative phases of the annular modes and (b) during all phases of the annular modes.
We identify ARs using the same method employed in Baek and Lora (2021); the following description of ARs is derived from therein with minor modifications. Our AR identification algorithm is originally derived from Lora et al. (2017) and described in Skinner et al. (2020) in the context of paleoclimate changes; it is specifically optimized for climate-change studies, taking into consideration characteristics of the background climate state in each case. The algorithm has been evaluated as part of the Atmospheric River Tracking Method Intercomparison Project (which compared AR-detection algorithms on a common data set; Shields et al., 2019;Rutz et al., 2019) and identifies moderate-to-strong ARs globally similarly to many other methods .
The AR detection relies on IVT and IWV, which are defined as vertically integrated horizontal moisture transport and vertically IWV content, respectively. We estimate daily IVT and IWV using the three available vertical levels at 850, 500, and 200 hPa available in the LME ensemble: Fields from the 850 hPa vertical level do not contribute to IVT or IWV for regions where monthly surface pressure does not exceed 675 hPa. Prior studies show that most of the moisture and moisture transport associated with ARs are found in the lower troposphere (Baek & Lora, 2021;Neiman, Ralph, Wick, Kuo, et al., 2008;Ralph & Dettinger, 2012;Ralph et al., 2005Ralph et al., , 2017 terms therefore likely account for the bulk of the IWV and IVT integrals, respectively, making our approach a reasonably good approximation.
ARs are defined as continuous regions, at least 2,000 km in length, where IVT exceeds a minimum threshold dependent on the background IWV: where IVT min = 225 kg m −1 s −1 and IWV is the 30-day running mean of zonal-mean IWV (Skinner et al., 2020). Our IVT threshold approximates those of various prior studies of ARs; moreover, our AR algorithm effectively distinguishes ARs (even those of relatively low-magnitude IVT) from high-moisture tropical regions. We apply our AR identification method on every ensemble member of the LME experiment; the ensemble mean of the LME is calculated after identifying ARs for each ensemble member. We define AR-driven IVT and precipitation as IVT and precipitation co-located and co-occurring with our identified ARs.
EKE is defined as: where u′ and v′ are bandpass filtered zonal and meridional winds, respectively. For the purposes of our study, we isolate u′ and v′ by applying a 10-day Butterworth high-pass filter to daily winds. We estimate EKE with the following: We compare the LME moisture transport against observations using the European Centre for Medium-Range Weather Forecasts Reanalysis fifth Generation data set (ERA5; Hersbach et al., 2020). The ERA5 Reanalysis is generated at hourly timesteps at 0.25° (∼30 km) horizontal resolution and over 137 vertical levels for the period 1979-present. For this study, the data set is sampled at every 3 hr at 0.5° (∼60 km) horizontal resolution on 30 vertical levels for 1980-2020. We identify ARs and calculate the annular modes using similar steps as those taken with the LME. When calculating the annular modes, however, we do not remove the 10-year running averages from longitudinally averaged zonal winds (we still remove seasonality); this is because the 40-year period analyzed is too short to experience meaningful long-term decadal variability. We also weight the zonal-mean of the zonal wind by the mass represented at each level (in addition to by the square root of the cosine of latitude).

Results
Positive phases of the NAM and SAM intensify and weaken ARs over the subpolar and subtropical latitudes, respectively, in the CESM1 LME ( Figure 1). Positive phases are defined such that the jet streams are displaced poleward; the poleward shifts in the storm tracks and primary AR activity are thus consistent with expectations. Negative phases of the annular modes show similar but opposite influences, with ARs weakening over the subpolar latitudes but intensifying over the subtropics (consistent with an equatorward shift of the jet stream). The annular modes also modulate AR frequency, with increases (decreases) in the number of ARs where ARs intensify (weaken). Notably, the magnitude of NAM's influence on AR activity exceeds that of SAM in intensity.
Though the primary influence of NAM is in the Northern Hemisphere, NAM also modulates ARs in the Southern Hemisphere. For instance, positive phases of NAM induce broad weakening of AR-driven IVT and precipitation in the Southern Hemisphere, with approximately opposite influences seen for negative phases of NAM.
Positive phases of NAM induce one dipole of AR activity over the North Pacific and another over the North Atlantic ( Figure 1). In western North America, the dipole over the North Pacific (a) reduces AR-driven IVT and precipitation over California and the Pacific Northwest but (b) increases AR-driven IVT and precipitation over British Columbia; in eastern Asia, the dipole reduces (increases) AR activity in the southern (northern) regions of Japan. Notably, the dipole over the North Atlantic drastically increases the frequency and intensity of landfalling ARs over Northern Europe during positive phases of NAM (with up to ∼15 mm/month increases in AR-driven precipitation) but decreases AR activity over southern Europe and southeast US. On the other hand, positive phases of SAM induce one zonal dipole of positive and negative AR frequency and intensity. The wet/dry dipoles flip in sign for negative phases of SAM, with approximately opposite but equal magnitude changes in AR activity.
Importantly, CESM1 successfully replicates observed NAM and SAM influences on ARs (Figure 2). Regressions of AR-driven IVT onto the annular modes from ERA5 show annular structures that are nearly identical across the LME and observational record, with positive (negative) phases of the annular modes intensifying (weakening) ARs over the subpolar latitudes but weakening (intensifying) ARs over the subtropical latitudes. Notwithstanding the stronger influence of NAM in observations, the CESM1 simulations successfully reproduce a NAM spatial structure that is less zonally symmetric than that of SAM. The CESM1 results furthermore successfully simulate SAM influences that are nearly identical both in structure and magnitude to those of observations. Collectively, our observational analyses indicate that the LME results are robust and representative of real-world influences.
To further investigate the role of the annular modes in modulating ARs, we perform EOF analyses on AR frequency (Figure 3). The EOF is performed (a) separately over the northern and southern extratropics (20°-70°N/S) and (b) member-by-member across the LME over the 850-2005 period. The leading pattern of AR frequency (accounting for 11% and 11%-12% of variance in the northern and southern extratropics, respectively) strongly resembles the AR influences shown in Figure 1. The first EOF reproduces the north-south dipole patterns present over the northern extratropical Pacific and Atlantic Oceans (the pattern correlation between EOF1 and regression of AR frequency onto NAM is 0.76 over the northern extratropics). It also reproduces the annular dipole structure spanning most of the southern extratropics (albeit with one annular band that is less robust than the other; pattern correlation between EOF1 and regression of AR frequency onto SAM is 0.62 in the southern extratropics). This leading pattern is furthermore stable across all 12 members of the LME, as demonstrated by pattern amplitudes corresponding to the first EOF that are uniform across the entire ensemble space. Nevertheless, the first EOF may not capture all of the annular modes' influences. For instance, the third EOF in the northern hemisphere also displays a dipole pattern over the North Atlantic, presumably capturing Atlantic jet shifts independent from the larger and stronger Pacific jet undulations (Hallam et al., 2022;Wettstein & Wallace, 2010), while the third EOF in the southern hemisphere displays a robust annular band over the southern subpolar latitudes that is only modestly represented in the first EOF. That the leading pattern of AR frequency variability in every member of the LME matches the spatial structure of NAM's (SAM's) influence on ARs in the Northern (Southern) Hemisphere is nevertheless a strong indication that the annular modes represent the primary modes by which ARs have varied over the last 1,000 years. Importantly, these EOFs (a) corroborate the spatial structures of the first two EOFs of winter AR variability in the Northern Hemisphere found by Ma and Chen (2022) and (b) do not change appreciably when anthropogenic influences over 1920-2005 are removed, indicating that our LME results provide a robust pre-industrial baseline.
We next consider the annular modes' influence on global moisture transport more broadly (of which ARs are a subset; Figure 4). SAM shows two distinct annular bands of moisture transport in the Southern Hemisphere, with one band in the subpolar latitudes (∼60°) and another in the subtropics (∼40°) that align with the latitudinal centers of the annular modes (Hartmann & Lo, 1998;Thompson & Woodworth, 2014). NAM also shows two bands of moisture transport that are less annular than SAM counterparts, as well as a third, zonally asymmetric band over the tropics (∼10°). These moisture transport patterns appear wind-driven: SAM and NAM display a dipolar and tripolar wind structure, respectively, corresponding to the two and three bands of moisture transport (and consistent with the spatial structure of the present-day modes; Thompson & Li, 2015;Thompson & Woodworth, 2014). Mean north-south displacements in zonal winds induced by annular modes (corresponding to latitudinal shifts in jet stream position) thus appear to be the primary mode of variability in moisture transport and AR activity. Notably, SAM's influence on moisture transport is negligible in the Northern Hemisphere, while NAM's influence, though concentrated in the Northern Hemisphere, is global in scope and affects the Southern Hemispheric subtropics.
The more annular structure of SAM is due to a comparative lack of continents in the southern extratropical latitudes, which creates longitudinally symmetric conditions and modest north-south meanders of the southern jet stream relative to its northern counterpart (Fogt & Marshall, 2020;Gerber & Thompson, 2017;Kidson et al., 2009). Topography also may contribute to NAM exhibiting greater sensitivity to the seasonal cycle than SAM. The SD of NAM during winter (SD = 1.31), for instance, is nearly twice that during summer (SD = 0.67), indicating that the extrema of NAM are both far more common and far more extreme during winter ( Figure 5; note that the seasonal cycle has already been removed from the NAM and SAM indices). This is likely because disruptions of the jet stream by topography-induced waves intensify in winter (due to a strengthening of the jet itself), contributing to larger undulations in winter than summer (Held et al., 2002;H. Wang & Ting, 1999;Woolings et al., 2010). SAM, in contrast, demonstrates much lower sensitivity to the seasonal cycle, with distributions that are more similar throughout the seasons.
Because of their statistical and dynamical links to extratropical cyclones, ARs are sometimes characterized in the literature as a storm track process (as opposed to a distinct phenomenon; Shaw et al., 2016). We separately  examine the annular modes' influence on storm track activity (summarized by EKE) to see if there are distinctions in how the annular modes modulate ARs and extratropical cyclones ( Figure 6). We find that poleward (equatorward) shifts in the jets are accompanied by increases (decreases) in storm track activity. For SAM, storm track activity manifests as one clear annular band over the southern extratropics and a less robust, more regional band over the subtropical Pacific. NAM also shows one annular band (though less annular than SAM) in the northern extratropics, and much weaker regional eddy activity over the subtropical North Atlantic. Notably, NAM's influence on storm track processes is limited to the Northern Hemisphere (unlike with AR or moisture transport). These monopolar storm track changes (as opposed to the dipolar or tripolar AR and moisture transport changes) are furthermore not latitudinally located between the dipolar or tripolar AR structures outlined above. The structural differences between a monopolar EKE signature and dipolar AR signature suggest that the annular modes' influence on storm track intensity is distinct from that on ARs and moisture transport, despite the storm tracks being associated with ARs and moisture transport.
Finally, we examine how poleward-equatorward shifts in mid-latitudinal moisture transport induced by the annular modes have changed in response to forced variability over the past millennium (Figure 7). Annual averages of meridional IVT at 50°N and 50°S show forced variability (represented by the LME ensemble mean) that is strongly correlated in the poleward direction (r = 0.51, p < 0.01). Such tandem latitudinal shifts in moisture transport are consistent with poleward (equatorward) shifts of the jet streams occurring with global warming (cooling) (Barnes & Polvani, 2013;Barnes & Screen, 2015;Yin, 2005). Global-mean temperature changes from volcanic activity, for instance, correlate strongly with poleward meridional IVT at 50°N (r = 0.59, p < 0.01) and 50°S (r = −0.32, p < 0.01), with sharp global cooling events induced by volcanic episodes (such as the three largest volcanic cooling events in 1259, 1453, and 1816), inducing relatively equatorward moisture transport in both hemispheres. Consistently, anthropogenic global warming post ∼1920 also induces stronger poleward moisture transport in both hemispheres (but particularly so in the Southern Hemisphere). Indeed, forced meridional IVT at 50°N/S appears to reach a maximum around 2005, indicating a strong role for anthropogenic forcings. Collectively, our analyses of forced meridional moisture transport provide a background context from which to evaluate past and future moisture transport change.

Discussion and Conclusion
Using the LME, we have characterized how ARs have responded to extratropical variability (both natural and anthropogenically forced) over the last millennium. Prior investigations into the relationship between ARs and extratropical variability have focused almost exclusively on the regional role of (a) northern extratropical variability on European and North American ARs (e.g., Brands et al., 2017;Guan et al., 2013;Lee et al., 2022;Zavadoff & Kirtman, 2020) and (b) southern extratropical variability on Australian and Antarctic ARs (e.g., Meneghini et al., 2007;Pohl et al., 2021;Reid et al., 2022;Wille et al., 2021). We add to these regional studies by providing new perspective on the global role of the NAM and SAM on ARs: mean northsouth displacements in zonal winds induced by the annular modes potentially comprise the primary driver of AR variability, with positive (negative) phases of the annular modes intensifying (weakening) and weakening (intensifying) ARs over the subpolar and subtropical latitudes, respectively. Notably, our results demonstrate that NAM has near global influences on ARs that extend to the Southern Hemisphere, with positive (negative) phases of NAM weakening (strengthening) ARs in the Southern Hemisphere. SAM's influence on ARs, on the other hand, is largely confined to the Southern Hemisphere, with only muted fingerprints in the Northern Hemisphere.
Although many studies have examined AR activity over the observational period (e.g., Baggett et al., 2017;Blamey et al., 2018;Lavers et al., 2011Lavers et al., , 2012Neiman, Ralph, Wick, Lundquist, & Dettinger, 2008;Rutz et al., 2019) and AR intensification in the future interval (e.g., Espinoza et al., 2018;Gao et al., 2015Gao et al., , 2016Hagos et al., 2016;O'Brien et al., 2022;Payne & Magnusdottir, 2015;Payne et al., 2020), few have investigated ARs into the paleoclimate past (e.g., Lora et al., 2017;Menemenlis et al., 2021;Shields et al., 2021;Skinner et al., 2020), particularly for the last millennium. Past assessments are nevertheless necessary in order to place modern and future AR changes in paleoclimatic context. In examining the daily available, 12-member LME over the entire 850-2005 interval, we provide over 5 million days of AR sampling across the full range of natural and anthropogenic radiative forcings over the last millennium (note, for instance, that this is not possible with a pre-industrial control run, which does not have varying forcings). Our work thus serves as a useful addition to the relative dearth of literature characterizing past AR variability.
Importantly, our efforts also provide a robust baseline from which to contextualize projected 21st century AR changes. Key among our findings is the revelation that phase changes in the annular modes induce up to ∼20-25 mm/month changes in AR-driven precipitation. While climatological AR activity has not changed much over the historical interval (due to influences from aerosols; Baek & Lora, 2021), ARs are widely expected to intensify in the future due to global warming (Espinoza et al., 2018;Gao et al., 2015Gao et al., , 2016Hagos et al., 2016;Payne & Magnusdottir, 2015;Payne et al., 2020). AR-driven precipitation increases under high-emission global warming (i.e., Representative Concentration Pathway 8.5 W/m 2 ), for instance, are projected to be approximately 20 mm/month over the 2006-2080 interval (Baek & Lora, 2021). The ∼20-25 mm/month changes in AR-driven precipitation induced by the annular modes thus (a) potentially comprise the largest natural changes in AR variability and (b) are comparable in magnitude to projected AR-driven precipitation increases in the coming decades. This is alarming, as it indicates that AR intensification from severe global warming over the next few decades may rival in magnitude the most prominent natural extremes of AR activity.
The primary role of the annular modes elucidated in our results furthermore holds important implications for (a) subseasonal-to-seasonal (S2S) predictability of ARs and, by extension, (b) water resource and flood risk management. Given the outsized importance of ARs as mechanisms of (sometimes extreme) water delivery to coastal regions (Corringham et al., 2019;Newman et al., 2012;Ralph et al., 2017;Zhu & Newell, 1998), the potential to forecast AR events is crucial to managing on-the-ground impacts of ARs. Skillful S2S prediction of ARs, for instance, can be accomplished with advanced knowledge of prominent modes of tropical variability (such as the Madden-Julian Oscillation and quasi-biennial oscillation; Baggett et al., 2017;Mundhenk et al., 2018). The annular modes too can be skillfully predicted on S2S timescales, though more advanced forecasts are difficult because of chaotic weather-scale atmospheric motions (Albers & Newman, 2021;Klavans et al., 2021;Marshall et al., 2012;Seviour et al., 2014). We demonstrate here that AR activity is intimately tied to the annular modes: this unequivocal role indicates annular mode forecasts may also skillfully predict AR activity with S2S lead times.
Finally, our finding that the annular modes imprint a dipolar signature on ARs but a monopolar signature on EKE suggests that ARs and extratropical cyclones are related but distinct processes. This is corroborated by the large correlations between the primary structures of variability of the ARs themselves and the NAM/SAM (Figure 3). Although the storm tracks are guided by the jet location (Shaw et al., 2016), the primary mode of EKE comprises amplitude variations (intensification/dampening) of the entire storm tracks as opposed to dipolar latitudinal shifts induced by shifts in the jets Thompson & Li, 2015;Thompson & Woodworth, 2014). Though not directly investigated in this study, the disconnect in why the primary mode described in the zonal-mean winds (a dipolar structure) is not the same as that described by the zonal-mean EKE (a monopolar structure) remains an open source of inquiry Lorenz & Hartmann, 2001;L. Wang & Nakamura, 2015Robinson, 2000;Thompson & Woodworth, 2014).