North American Winter Dipole: Observed and Simulated Changes in Circulations

: In recent years, a pair of large-scale circulation patterns consisting of an anomalous ridge over northwestern North America and trough over northeastern North America was found to accompany extreme winter weather events such as the 2013–2015 California drought and eastern U.S. cold outbreaks. Referred to as the North American winter dipole (NAWD), previous studies have found both a marked natural variability and a warming-induced ampliﬁcation trend in the NAWD. In this study, we utilized multiple global reanalysis datasets and existing climate model simulations to examine the variability of the winter planetary wave patterns over North America and to better understand how it is likely to change in the future. We compared between pre- and post-1980 periods to identify changes to the circulation variations based on empirical analysis. It was found that the leading pattern of the winter planetary waves has changed, from the Paciﬁc–North America (PNA) mode to a spatially shifted mode such as NAWD. Further, the potential inﬂuence of global warming on NAWD was examined using multiple climate model simulations.


Introduction
Cold-air intrusion events in eastern North America are common during boreal winter, but severe events, with costly consequences, used to be rare. According to the National Centers for Environmental Information report (NOAA 2019) [1], during 1985-2005, Florida was the only state suffering billion-dollar losses related to freeze events. After 2006, the number of eastern states experiencing billion-dollar freeze-related disasters has increased threefold. This trend is concurrent with a marked increase in the number of billion-dollar drought disasters in western states. The exceptionally dry conditions in the 2013-2014 winter and the associated high-amplitude atmospheric ridge in the Northeast Pacific was accompanied by an abnormally deepened trough in the eastern USA with strong cold-air outbreaks. This striking division of severe drought in the west and extreme cold in the east has been referred to as the "North American Winter Temperature Dipole" (Singh et al., 2016) [2], while its accompanying circulation pattern is called the "North American Winter Dipole" (NAWD)   [3].
The polarity and location of this circulation NAWD coincide with the wintertime stationary waves over North America, which feature a high-pressure ridge in the west and a low-pressure trough in the east (Wang et al., 2017, Stuivenvolt Allen and Wang, 2019) [4,5]. Oscillating in sync with the stationary waves, the positive-phase NAWD is associated with an anomalous ridge over the Gulf of Alaska and a deepened trough near the Great Lakes, thereby enhancing the east-west temperature contrast in North America including prolonged freezing in the Great Lakes (Voelker et al., 2019) [6]. Past studies suggested that the ridge in the western U.S. (Swain et al., 2014) [7] and the temperature dipole itself have amplified (Singh et al., 2016) [2], owing to internal and external variabilities such as tropical Pacific heating (Hartmann 2015; Schulte and Lee 2017) [8,9] and Arctic warming (Francis and Varus 2012;Overland et al., 2016) [10,11]. The intense reversal of NAWD seen in the winter of 2016-2017 and associated extreme precipitation in California suggests that the NAWD variability has also increased (Wang et al., 2017;Swain et al., 2018; [4,12,13]. This study explored the wintertime leading pattern evolution across North America and sought to identify potential forcing sources through diagnostic analyses. While the 2013/2014 extreme NAWD event has been extensively studied, the role of large-scale climate variability in terms of the wintertime stationary waves has not been satisfactorily addressed. The overarching goal here was to examine the dominant patterns during the recent decades and to compare that with the known dominant climate mode, including the Pacific-North American (PNA) pattern. Given the reportedly weak correlation between the PNA and NAWD   [13,14], the two phenomena may be independent and their collective changes in amplitude over time deserves research attention.
In this study, multiple global reanalysis and observational datasets were used alongside global climate model simulations, which are introduced in Section 2. In Section 3, we provide discussions of the analysis results. Concluding remarks are presented in Section 4.

Observations
The primary dataset used to study the atmospheric circulation features of interest is the NCEP-NCAR reanalysis (Kalnay et al., 1996) [15] from 1948 to the present, given its consistent data assimilation utilized throughout the entire period. Although satellite data were lacking before 1979, the temporal coverage of this reanalysis is longer than newer reanalysis that ingests satellite data. We also used European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-Interim; Dee et al., 2011) [16] from 1979 to 2015 and extend back to 1958 by concatenating with ERA-40 (Uppala et al., 2005) [17]. Additionally, we used a retrospective reanalysis dataset to depict the global atmospheric circulation: the ECMWF global reanalysis products for the twentieth century (ERA-20C) (Poli et al., 2016) [18] from 1871 to 2012 and the National Oceanic and Atmospheric Administration (NOAA) 20-Century Reanalysis data. Based on large ensembles, 20CR (ERA-20C) utilized a coupled Atmosphere/Land-surface/Ocean model to assimilate surface observations of sea level pressure (surface pressure) to produce the atmospheric conditions. It should be noted that the 20CR dataset is developed from the model simulation forced by surface conditions only; therefore, it is not a full reanalysis. In this study, we utilized the monthly mean North Hemisphere geopotential height from November to February on 2.5 • × 2.5 • long grid spacing. Regarding sea surface temperature (SST), NOAA Extended Reconstructed Sea Surface Temperature (ERSST) V4 dataset was used in this study, starting in 1854 with a 2 • × 2 • spatial resolution (Huang et al., 2014) [19].

Simulations
Model simulation outputs were analyzed from two sources: (1) the Community Earth System Model version 1 (CESM1) Large-Ensemble Project with 40 members (Kay et al., 2015) [20] from 1950 to 2005 for the historical period and from 2040 to 2080 for the projected perspectives driven by Representative Concentration Pathways (RCP) 8.5 scenario; and (2) the 21,000-year-long Paleoclimate simulation from 21,000 years before present (1950 A.D.) produced by the model ECBilt-CLIO, which was forced with time-dependent ice-sheet topography, orbital forcing, and greenhouse gas forcing.
Based on a quasi-geostrophic adiabatic core with T21 resolution (64 points for longitude and 32 points for latitude) and three vertical layers (200 mb, 500 mb, and 800 mb), the ECBilt-CLIO includes the prognostic vorticity equation and the thermodynamic equation with a set of physical parameterizations of diabatic processes and coupled with an ocean model. The Paleoclimate simulation started from a 2000-year-long Last Glacial Maximum equilibrium simulation with pre-industrial greenhouse gases. Here, the geopotential height is only available at 500 mb for wintertime (DJF) means.

NAWD-Related Climate Indices
We further considered a number of documented large-scale climate oscillation indices that were linked to the significant 2013/2014 and 2014/2015 winter anomalies in North America (Table 1). Western North Pacific (WNP) is defined as a specific SSTA pattern that forms one year before a full-fledged El Nino/La Nina, with a growing connection between the WNP pattern and the development of El Niño-Southern Oscillation (ENSO) in the following year since 1960 (Pegion and Selman 2017) [21]. By using the Niño 4 index in the year prior an ENSO event, denoted as Nino4 (Y + 1) fluctuation, Wang et al. (2014) [13] found a close association of the NAWD index with Nino4 (Y + 1). Next, the Blob index describes the Pacific warm water patch in the northeast Pacific during the 2014/2015 California drought event (Bond et al., 2015) [22] owing to the anomalous high-pressure system that works against the prevailing surface westerlies and subsequent reduction in local surface evaporation. Further, Lorenzo and Mantua (2016) [23] identified the Gulf of Alaska (GOA) SST anomaly as an index to depict warming in the Pacific coastal boundary of North America. Schulte and Lee (2017) [9] proposed the East Pacific/North Pacific (EP/NP) pattern citing its relationship with USA temperature variability since the 1950s. The Arctic influence was also considered in the 2013-2014 winter event, when the Arctic Oscillation (AO) reached a strong positive phase   [13] and apparently affected the eastern and southern part of the US (Lorenzo and Mantua 2016). Finally, we follow Singh et al.'s (2016) definition of the North American Winter Temperature Dipole index by computing the difference of surface temperature between the western and eastern USA. All these indices were documented to have certain connection with the extreme 2013-2014 winter anomalies, representing many processes that involve air-sea interaction and internal climate variability (Table 1).
For the index of NAWD, we followed Wang et al. (2017) [4] by subtracting the monthly anomalies of geopotential height at 250 hPa values between the ridge and trough centers from the 5 • × 5 • domain averages centered at (53 • W, 54 • N) and (103 • W, 63 • N), respectively, to depict the standing pattern of the upper-level circulations that is in-phase with the winter stationary waves. Note that the NAWD index in the paleoclimate simulation data differed from the reanalysis because of discrepancy of the vertical level. To make these data comparable, we used least squares regression, where the NAWD index from the 20CR reanalysis dataset was the independent variable and the NAWD index from the paleoclimate dataset was the dependent variable, computed over the overlapping time period of 1950-1981 for the two datasets.  (2017) The positive phase is associated with an amplified jet stream configuration in which a ridge is situated over Alaska and a trough is situated over the eastern US.

Long-Term Change in Leading Circulation Pattern across North America
To depict the stationary eddies, we computed monthly anomalies of geopotential height (Z) at 250 hPa from November to February, with zonal means removed, herein ∆Z E 250. We then subjected the Empirical orthogonal function (EOF) to running 30-year period of ∆Z E 250, in which each EOF analysis contained 120 realizations (30 years × 4 months), covering the Northern Hemisphere. This month-by-month arrangement of ∆Z E 250 reflects the sub-seasonal variability of winter climate anomalies, which is robust over North America and North Pacific (Higgins et al., 2000) [24]. NCEP-NCAR reanalysis was used as the initial analysis of ∆Z E 250.
By equally dividing the 1948-2018 period, we set 1980 as the mid-point to distinguish whether circulation patterns have changed. The first EOF mode of ∆Z E 250 during the earlier period  is shown in Figure 1a, explaining 18.6% of the variance. This EOF mode features a wave train emanating from the central Pacific to the U.S. and is coincident with the PNA pattern, which is consistent with the literature (Horel and Wallace 1980) [25]. By superimposing PNA contours (produced by correlating ∆Z E 250 with the PNA index from the NOAA Climate Prediction Center), both EOF1 and PNA patterns are spatially in-phase. By comparison, the post-1980 EOF1 (Figure 1b shading) shows a wave train of similar extent, but the wave centers are shifted from the PNA pattern by a quarter phase. It appears that the post-1980 EOF1 becomes phase-coincident with the NAWD centers marked with × and + in North America. This result suggests that the leading mode of ∆Z E 250 variability may have changed, from a PNA-dominant pattern to one that resembles NAWD. We further tested the rule of thumb provided by North et al. [26] for determining the eigenvalue with an approximate 95% confidence interval. Figure 2 shows the variance of the four modes for each period. The gap between the first and second modes in the early period is greater than in the late period, suggesting that the PNA was initially the dominant driver with a shift toward the NAWD becoming more of a factor (or at least as important as the PNA) in recent decades.
We note that, when applying EOF analysis on seasonal means (Nov-Feb average; Figure 3) instead of monthly intervals, the latter period's EOF1 remains similar to the PNA (Figure 3b), whereas NAWD lies closer to the second mode (Figure 3d). These results imply that the transition of the leading mode from PNA to NAWD is related to sub-seasonal variability, an emerging subject of intense research (National Academies of Sciences, Engineering, and Medicine. 2016). However, we also note that EOF2 of the earlier period did not reveal the NAWD pattern (Figure 3c), suggesting that its seasonal mode has grown stronger nonetheless in recent decades.
To delineate possible oceanic influences on these leading modes of atmospheric circulation variability, we correlated the first principal component (PC1) time series with monthly SST anomalies (Nov-Feb). As shown in Figure 1c, the pre-1980 EOF1 corresponds to a broad La Niña-like pattern resembling the cold-phase of the Pacific Decadal Oscillation (PDO). After 1980, the PC1 correlation pattern of SST changed dramatically, absent the PDO features alongside most of the tropical signatures (Figure 1d). This post-1980 SST pattern of PC1 reveals the oceanic "blob" along the US West Coast (Kintisch 2015) [27] and robust negative anomalies in the far end of the Western North Pacific and Bay of Bengal. This latter SST feature coincides with an emerging ENSO precursor, called the Western North Pacific pattern (WNP; Wang et al. 2012) [28], that saw an increased association with ENSO development in the following winter during recent decades (Wang et al., 2013; Capotondi and Sardeshmukh 2015) [29,30]. Marked difference between these SST patterns accompanying EOF1 (∆Z E 250) of the two periods delineates two distinct modes of subtropical variability that influences North American winter climate. These factors are discussed further in Section 3.3.

Evolution of Pattern Change
It is important to examine the timing in which the leading mode of ∆Z E 250 started to change. Adopting the running-EOF method of Zhang et al. (2008) [31], a series of EOF analysis was conducted with a 30-year window and repeated every five years. The leading EOFs were then subjected to spatial correlation analyses with the PNA and NAWD patterns, thus forming a series of correlation coefficients from each 30-year period. To validate the NCEP-NCAR reanalysis, these analyses were performed on the ERA20C and the ERA-Interim as well. As shown in Figure 4a, the running-EOFs of all three reanalysis datasets reveal a decreasing trend in the correlations between EOF1 and PNA from 1970 to 2000; meanwhile, the correlations increased between EOF1 and the NAWD pattern. Correspondingly, spatial correlations of the second running-EOF (EOF2) with the PNA and NAWD (Figure 4b) show opposite trends, suggesting that the NAWD, which used to be the second mode, intensified during the 1990s to the extent that it overtook the PNA as the leading mode.
The mechanism behind this reported decadal shift in the prevailing modes of variation is manifold. Previous studies proposed that an ENSO-forced PNA would shift eastward in response to warming of SST (Zhou et al., 2014) [32]. Interdecadal variability of the North Pacific sea level pressures can also induce a shift in the PNA (Johnson and Feldstein 2010) [33]. Moreover, one could relate NAWD to the maintenance of wintertime stationary waves in the Northern Hemisphere (e.g., Chang 2009) [34]. The western ridge of the North American stationary wave is primarily linked to orographic forcing of the Tibetan Plateau, whereas the eastern/downwind trough is amplified by orographic forcing imposed by the Rocky Mountains. Diabatic heating from the subtropical Western Pacific (i.e., the warm pool) and midlatitude North Pacific (i.e., the storm track) further enhances this ridge-trough pattern over North America (Chang 2009) [34]. Therefore, the observed NAWD amplification could be related to such changes in the combined forcing, i.e., jet stream-terrain interactions and diabatic heating. This aspect of the changing winter circulation maintenance over North America has not been explored in the wake of studies devoted to the 2013-2014 extreme anomalies. To reveal the associated forcing of teleconnection involved in the circulation pattern change, we conducted the analysis of Rossby Wave Activity Flux (WAF) derived by Takaya and Nakamura (2001) [35] at 250 hPa, regressed upon the PC1 and PC2 time series of Figure 1. A wave-like structure in the upper troposphere ( Figure 5) is associated with the wave activity fluxes that suggest two forcing sources: (a) one emanated from the tropical central Pacific during the pre-1980 period; and (b) one from the subtropical western Pacific during the post-1980 period. The downstream WAF over the eastern U.S. regenerate another wave train towards Europe, which is known to the PNA. The upstream WAF after 1980 shifts and expends northward to near Japan. The post-1980 wave train has a more robust cyclonic anomaly in the western Pacific, conveying the energy that enhances the ridge over the Gulf of Alaska while the wave train tracks through northeastern North America. Notably, this wave train is phase-coincident with NAWD (marked by × and + in Figure 5b) and does not re-generate downwind waves towards Europe.

Past NAWD
To gain the historical and projected perspectives for the transition between the leading EOF modes of atmospheric circulation, we next computed a 30-year running variance of the NAWD index by following the definition of Wang et al. (2014). Figure 6 shows the running variance of NAWD derived from the NOAA 20CR, the NCEP-NCAR and ERA-Interim reanalysis data, overlaid with the CESM ensemble simulations. NAWD has undergone pronounced inter-decadal variation on the order of 60 years and demonstrated the largest variance in the early 21st century. Similar analyses applied to each CESM ensemble member and averaged as the ensemble mean revealed an increase in the NAWD variance during the historical period (up to 2005 with increasing greenhouse gas), a trend that is projected to continue under the high-emission scenario (RCP 8.5) towards 2080. The projected ensemble-mean variance starts to decline after 2050, despite an increase in the ensemble spread. At this point, we can only attribute this late-21st-century decline of the NAWD variance to the natural variability revealed in the observed period.
As an alternative means to understand how the NAWD and stationary waves vary with respect to the global temperature change, we applied the same analytical approach to the paleoclimate simulations. Considering the much longer paleoclimate time scale, the running variance was calculated across a 100-year window, with a 10-year interval. Figure 7 displays the reconstructed paleoclimate NAWD index (blue line) and the global and Northern Hemispheric mean surface temperatures (green lines), superimposed with the 20CR's NAWD 100-year running variance from 1851 to 1950 (red dot in the right end) for comparison. The time series from 20CR and paleoclimate simulations did reveal similar fluctuations during their overlapping years. During the Pleistocene era, the global temperature is lower and NAWD is less variable; during the Holocene with increasing global temperature, the NAWD variance increases significantly and shows a larger fluctuation. The time period after the Younger Dryas event appears to be a turning period in which global temperature and NAWD variability change correspondingly for almost ten centuries. We should note that the paleoclimate dataset has a lower resolution so the NAWD index calculated may not be compatible in magnitude with that of the modern reanalysis data; further, a least-squares regression was applied to approximate the 500-hPa geopotential height to the 250-hPa NAWD index (see Section 2). During the last glacial period, North America was characterized by a thick layer of ice covering the area that now constitutes the eastern dipole of the NAWD; this would certainly have affected atmospheric circulation patterns beyond that expected from temperature alone. Large ice sheets remained just to the north of the eastern dipole region of central Canada until between 9000 and 7000 years BP and much of the Holocene vegetation types we are familiar with did not occur until 5000 years BP (Williams et al., 2004) [36]. As such, comparisons of NAWD variability to Northern Hemisphere temperature trends over the period from 5000 years BP until the present likely presents the most robust paleo analog for how temperature may affect NAWD into the future. From 5000 years BP to the present, ice core stable isotopes and palynological studies indicate that temperatures across the Northern Hemisphere have slowly declined, and that NAWD variance has also declined ( Figure 7). Overall, reconstructions of paleo temperature and simulations support an interpretation of greater energy and temperatures across the Northern hemisphere promoting greater NAWD variability.
We should caution that this paleoclimate model lacks feedback processes in the carbon cycle, vegetation, and terrestrial ice sheets; moreover, ice sheet height of Greenland is different in the past so the circulation-terrain interaction might have impacted the jet stream and subsequent NAWD. It is also possible that the stationary waves' position experienced some shifting. These aspects require more complete simulations to reveal in future research.

Possible Forcing Change
Here, we examine the documented large-scale mode of climate variability linked to the significant 2013-2015 winter anomalies in North America. Figure 8a shows the 30-year running correlation analysis of the NAWD index with the list of relevant climate indices as shown in Table 1 (see Section 2.3). It was found that the Blob, GOA, Niño 4 (Y + 1), and WNP indices are significantly correlated with NAWD (p < 0.05) during the latter 20th century, consistent with the literature. These significant yet episodic correlations imply that NAWD is mostly linked to local SST anomalies in the northeast and northwest Pacific, which is in good agreement with Figure 1d. The increasingly high correlation between the Blob index/GOA index and the Dipole index since 1960 indicates a strengthened relationship of NAWD with SST variations. As mentioned in the last section, the northeast Pacific warm SST patch (or Blob) is induced by an anomalous ridge of high-pressure, which brings cloud-free conditions and excess solar radiation that warms the sea surface. Lorenzo and Mantua (2016) [23] indicated that the correlation between the EP/NP index and 300 hPa stream function resembles NAWD. However, the correlation between the EP/NP index and NAWD here does not exceed the significance level. As their study suggested, the relationship between EP/NP and the circulation is clearer in the eastern and northern part of the US, but the NAWD index represents the higher-latitude ridge/trough system, which may lead to the insignificant correlation. Lorenzo and Mantua (2016) [23] and Wang et al. (2014) [4] also mentioned that the AO index has a stronger relationship with the southern U.S. This may explain why the correlations of AO and EP/NP with NAWD do not reach statistical significance, but nonetheless show an increasing trend after the mid-1970s. The other indices do not show significant correlations but some do exhibit a growing relationship with the NAWD; this observation infers that constructive interactions between these climate modes could be a compounding mechanism that strengthened NAWD in winter 2013/2014 (Lee et al., 2014) [37].
Here we provide a provisional discussion on the possible role of inter-decadal climate variability in modulating Figure 8a's low-frequency variations in the sliding correlations: the Atlantic multidecadal oscillation (AMO) and the Pacific decadal oscillation (PDO). It has been suggested that the AMO can impact climate variations in the Pacific Ocean (e.g., Zhang and Delworth 2007) [38]. Figure 8b shows the 5-year running mean of the AMO and PDO indices and, upon visual examination, it appears that the AMO has an inverse, yet corresponding, phasing with the correlations of GOA and the Blob (green lines in Figure 8a). This analysis, though preliminary, suggests a possible role of the natural variability (i.e., the AMO) in modulating the SST-NAWD co-variability in the long run; this feature requires further research.

Conclusions
The diagnostics undertaken here suggest that the leading mode of the wintertime atmospheric stationary waves in the Northern Hemisphere has undergone a notable change. Since the 1980s, the first leading mode of anomalous winter stationary eddies in the subseasonal timescale has alternated from the PNA to a comparably strong NAWD. Given that EOF describes the variance of individual patterns, the observed increase in the amplitude of NAWD suggests that it may have overtaken the PNA as the more common type of winter variability, with increased influence from western North Pacific SSTs. The CESM large-ensemble simulations forced with increasing greenhouse gas indicate that the NAWD variance will gradually amplify alongside its low-frequency natural variability. This result implies that the subseasonal variation of the atmospheric circulations over North America could continue to be dominated by the NAWD mode, with the potential to sharpen the east-west temperature and precipitation division across North America.
Further comparison of EOF1/EOF2 with PNA/NAWD using spatial correlation reveals a pattern change that is consistent among different reanalysis datasets. Analysis of the paleoclimate simulations suggests that NAWD does co-vary with the global temperature variation at centennial timescales: the NAWD variance was weak in the cooler Pleistocene and gradually increased and stabilized in the warmer Holocene. We further examined large-scale forcing of both atmospheric and oceanic origins for the amplified NAWD and found that most of the documented climate indices, which showed a connection to the 2013-2014 extreme winter of North America, are insignificantly correlated with NAWD with the exception of certain oceanic features like the WNP and Niño 4 (Y + 1); these imply that regional SST anomalies in the subtropical West Pacific may provoke teleconnections that affect NAWD. Future examination of the dynamic processes leading to the NAWD amplification should consider the stationary waves maintenance, namely the jet-terrain interactions, tropical and extratropical diabatic heating, and the effect of Arctic amplification.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: