Change in the variability in the Western Pacific pattern during boreal winter: roles of tropical Pacific sea surface temperature anomalies and North Pacific storm track activity

Using multiple reanalysis datasets, this study reveals that the variability in the western Pacific pattern (WP) in boreal winter has shown notable changes during recent decades. The variability in the winter WP exhibited a marked weakening trend before the early 2000s and increased slightly thereafter. Two epochs with the highest and lowest WP variabilities are selected for a comparative analysis. Winter WP-related meridional dipole atmospheric anomalies over the North Pacific were stronger and had a broader range during the high-variability epoch than during the low-variability epoch. Correspondingly, the winter WP had larger impacts on surface temperature variations over Eurasia and North America during the high-variability epoch than that during the low-variability epoch. We find that the shift in the winter WP variability is closely related to changes in the connection between the winter WP and the El Niño-Southern Oscillation (ENSO) and to changes in the amplitude of the North Pacific storm track. Specifically, ENSO had a closer connection with the WP during the high-variability epoch, when the amplitude of the North Pacific storm track was also stronger. During the high-variability epoch, the extratropical atmospheric anomalies generated by the tropical ENSO shifted westward and projected more on the WP-related atmospheric anomalies, thus contributing to an increase in WP variability. In addition, the larger amplitude of the North Pacific storm track that occurred during the high-variability epoch led to the stronger feedback of synoptic-scale eddies to the mean flow and contributed to stronger WP variability. Further analysis indicates that the change in the connection of ENSO with the WP may be partly related to the zonal shift of the sea surface temperature anomaly in the tropical Pacific associated with ENSO.


Introduction
The western Pacific pattern (WP), one of the most prominent teleconnections in the Northern Hemisphere, features a marked north-south seesaw pattern in 500-hPa geopotential height anomalies over the Far East and western North Pacific (Wallace and Gutzler 1981;Mo and Livezey 1986;Barnston and Livezey 1987). A growing body of literature highlights the distinct impacts of the WP on weather and climate anomalies in Eurasia and North America. For example, studies have indicated that the WP exerts pronounced impacts on the East Asian winter monsoon, the occurrence frequency of the East Asian cold surge, and surface air temperature (SAT) anomalies over East Asia and North America (e.g. Nakamura 2005a, b, 2013;Linkin and Nigam 2008;Wang and Chen 2014;Chen and Song 2018;Aru et al. 2021). In addition, several studies have reported that the WP has a notable modulating effect on the connection between the Arctic Oscillation (AO) and SAT anomalies over East Asia (Cheung et al. 2016;Oh et al. 2017;Lim and Kim 2013). In particular, the spatial distributions and amplitudes of AO-generated winter SAT anomalies over East Asia strongly depend on the phase of the WP (Park and Ahn 2016;Oh et al. 2017). Studies have argued that the WP-induced SAT anomalies over East Asia and North America are comparable to, or are even stronger than, those associated with the El Niño-Southern Oscillation (ENSO) and the Pacific-North American teleconnection pattern (Cheung et al. 2016;Oh et al. 2017;Lim and Kim 2013). Furthermore, variability in the WP has a marked effect on the jet stream over the North Pacific Tanaka et al. 2016;Ma and Zhang 2018), as well as on sea ice conditions in the Arctic Ocean, especially in the Bering and Okhotsk Seas (Linkin and Nigam 2008;Yuan et al. 2015). Due to the crucial impact of the WP, investigating the factors responsible for WP variability is of great importance.
There are two main lines of thought concerning the mechanisms of the formation and maintenance of the WP. It is generally accepted that the WP is an intrinsic atmospheric mode that can form and sustain itself in the absence of external forcing. Studies have demonstrated that the interaction between storm tracks (i.e., synoptic-scale eddies) and mean flows, as well as the associated feedback of synoptic-scale eddies, are the essential mechanisms of the formation and maintenance of the WP (Hsu and Wallace 1985;Nakamura et al. 1987;Lau 1988;Lau and Nath 1991;Wettstein and Wallace 2010;Athanasiadis et al. 2010;Li and Wettstein 2012). Rivière (2010) highlighted a dynamic link between Rossby wave breaking and the WP that hinges on the latitudinal fluctuations of the Pacific jet stream. In addition, studies have indicated that the extraction of kinetic and available potential energy from the basic flow via barotropic and baroclinic energy conversion processes plays a crucial role in the maintenance of the WP (Tanaka et al. 2016;Sung et al. 2019Sung et al. , 2020. In light of the recent work by Dai and Tan (2019), the formation of WP can also be attributed to the precursory disturbances over the North Pacific and East Asia. There is also a possibility that a wintertime WP could originate from a PNA-like disturbance over the North Pacific. Consistent with Dai and Tan (2019), diagnostic analysis of the relative roles of various energy conversion processes in the formation and maintenance of the WP further indicated that both baroclinic energy conversion from the background flow and transient eddies feedback forcing facilitate the formation of the WP (Zhuge and Tan 2021).
On the other hand, studies have argued that external forcings (such as tropical sea surface temperatures (SSTs), Eurasian snow cover, and Arctic sea ice anomalies) also play a role in the excitation and maintenance of the WP (Horel and Wallace 1981;Kodera 1998;Koide and Kodera 1999;Linkin and Nigam 2008;Hirose et al. 2009;Frankignoul et al. 2011;Oshika et al. 2014;Nakamura et al. 2015). Horel and Wallace (1981) revealed that SST anomalies in the tropical Pacific related to ENSO contributed to the formation of the WP via a heating-induced atmospheric Rossby wave train. Kodera (1998) also claimed that a link exists between the WP and the upstream circulation conditions over the Eurasian continent during early winter. Evidence suggests that the Arctic sea ice anomalies generated by the North Atlantic Oscillation in the preceding winter could induce a WP-like atmospheric anomaly pattern in the following winter (Nakamura et al. 2015). Kodera (1998) argued that snow cover over Siberia can affect the WP by modulating the WP-ENSO connection. Some studies have revealed that the relationships between the WP and these external forcings are unstable (Furtado et al. 2012;Park et al. 2018;Xu and Fan 2020).
As highlighted above, many studies have investigated the interannual variability in the WP as well as the impacts of the WP on climate anomalies over Eurasia and North America. At present, few studies have discussed whether the historical variabilities (i.e., amplitude) in the WP are stable. The impacts of the WP on climatic anomalies may depend strongly on the amplitude of the WP. For example, recent studies have indicated that an extremely strong winter WP event was the largest contributor to the occurrence of recordbreaking cold events over North America during the winter of 2013/2014 (Yu and Zhang 2015;Baxter and Nigam 2015;Lee et al. 2015). Cheung et al. (2015) reported that a WPlike pattern notably accounted for the persistent cold spells that occurred in Hong Kong over 1948/1949. In addition, Sung et al. (2019 demonstrated that changes in climate extremes over North America are attributable to changes in WP variations. Hence, it is critical to examine interdecadal changes in the historical variability in the WP, as these changes have important implications for the prediction and projection of mean and extreme climate conditions over Eurasia and North America.
The remainder of this study is structured as follows. Section 2 describes the data and methods. Section 3 examines the interdecadal changes in variability in the WP. In addition, in this section, we compare atmospheric and surface air temperature anomalies associated with the winter WP before and after the interdecadal changes occurred. Section 4 examines the factors responsible for changes in WP variability. The summary and discussion are shown in Sect. 5.

Data
Monthly and daily geopotential height and horizontal wind data were extracted from the National Centers for Environmental Prediction (NCEP)-National Center for Atmospheric Research (NCAR) reanalysis dataset (Kalnay et al. 1996). The NCEP-NCAR reanalysis dataset is available from 1948 to the present with atmospheric variables at various pressure levels gridded at a 2.5° × 2.5° horizontal resolution. For a cross-check, we also compared the results obtained from other datasets, including atmospheric data with a horizontal resolution of 1.0° × 1.0° spanning the period from January 1979 to August 2019 provided by the European Centre for Medium-Range Weather Forecasts Interim (ERA-Interim) reanalysis product (Dee et al. 2011) and atmospheric data with a horizontal resolution of 1.25° × 1.25° produced by the second Japanese global atmospheric reanalysis project (JRA55) (Kobayashi et al. 2015) from January 1958 to the present. For validation based on a longer time period, we adopted the NOAA-CIRES-DOE Twentieth Century Reanalysis (V3) data with a resolution of 1.0° × 1.0° for monthly geopotential height and skin surface temperature, provided by the NOAA/OAR/ESRL PSL, Boulder, Colorado, USA. The NOAA-CIRES-DOE Twentieth Century Reanalysis (V3) data, covering the period from 1836 to 2015, uses a newer, higher-resolution model and assimilates a larger set of observations (Slivinski et al. 2019).
Monthly mean surface air temperature data were collected from the University of Delaware (Willmott and Matsuura 2001) with a horizontal resolution of 0.5° × 0.5° and availability from January 1900 to December 2017. The monthly SST data used in this study were derived from the National Oceanic and Atmospheric Administration (NOAA) Extended Reconstructed SST version 5 (ERSST.v5) dataset (Huang et al. 2017). This dataset covers the period from January 1854 to the present and has a 2° × 2° horizontal resolution. The datasets employed in this study are summarized in Table 1.

Identification of the employed indices
Following previous studies (Nishii et al. 2010;Aru et al. 2021), the WP is described by the first empirical orthogonal function mode (EOF1) of the 500-hPa geopotential height anomalies (Z500) over the western North Pacific region (20°N-70°N, 120°E-180°). The principal component (PC) time series corresponding to EOF1 is considered the WP index. Notably, the geopotential height anomalies were weighted by the square root of the cosine latitude prior to the EOF analysis to account for area differences occurring at different latitudes (Thompson and Wallace 2000). As indicated in previous studies, storm tracks can be highlighted by various eddy variance and covariance quantities (Blackmon et al. 1977;Chang and Fu 2002;Wettstein and Wallace 2010). In this study, the storm track activity (STA) is described based on the 2-8-day bandpassfiltered 500-hPa geopotential height variance (referred to as Z ′ Z ′ 500 , where the prime indicates the filtering anomalies and the overbar denotes the winter mean). Following previous studies Zhao et al. 2019Zhao et al. , 2020Piao et al. 2020), the Niño3.4 SST index was used to characterize the ENSO variability; this index is defined as the area-average SST over the area comprising 5° S-5° N and 170°-120° W.
The analysis period spanned from 1949 to 2019, including seventy-one winters. Winter is considered the average of December, January, and February (DJF). Here, the winter of 1949 refers to December in 1948 and January and February in 1949. The interannual variations in all variables were obtained by subjecting the original variables to a 9-year high-pass Lanczos filter (Duchon 1979). Significant composite difference, correlation and regression coefficient levels were estimated based on two-tailed Student's t tests. In addition, a two-tailed F test was applied to test the significance levels of the differences between two standard deviations. The statistical F test uses a F statistic to compare two variances or two standard deviations developed by Ronald A. Fisher (Von Storch and Zwiers 2001) and is also a measure of the ratio variance. The F statistic is defined as: where, in this context, x 1 and x 2 denote the same variable at different time period, respectively. If the F is larger than the value in the F table, it can be assumed that the variances/ standard deviations of the variable have varied considerably over time.

Geopotential height tendency
The geopotential height tendency is used to quantitatively describe the feedback of synoptic-scale eddies to low-frequency flows (Lau and Holopainen 1984;Lau 1988;Cai et al. 2007;Chen et al. 2015Chen et al. , 2020aHu et al. 2021) and can be expressed as follows:  Kalnay et al. (1996) 1948 to present JRA55 Monthly geopotential height Kobayashi et al. (2015) 1958 to present ERA-Interim Monthly geopotential height Dee et al. (2011) 1979-2019 Monthly SST Huang et al. (2017) 1854 to present NOAA-CIRES-DOE Monthly geopotential height and SST Slivinski et al. (2019Slivinski et al. ( ) 1836Slivinski et al. ( -2015 where g is the acceleration of gravity, f denotes the Coriolis parameter, and ��� ⃗ V ′ and ′ correspond to the synoptic-scale winds and vorticity, respectively. Moreover, ′ is calculated according to synoptic-scale winds. Synoptic-scale variations are obtained by subjecting the daily mean data to a 2-8-day bandpass Lanczos filter (Duchon 1979).

Temporal variation features
Before investigating the changes in WP variability, we first examined the spatiotemporal features of the WP during boreal winter based on the NCEP-NCAR reanalysis data over the 1949-2019 period. Here, the spatial pattern of the WP was represented by the EOF1 of the Z500 anomalies over 20° N-70° N and 120° E-180°. The EOF1 of the Z500 anomalies in boreal winter over 1949-2019 accounted for 44.3% of the total variance. Figure 1a depicts a map of Z500 anomalies regressed onto the winter WP index. The spatial pattern of the WP features a meridional seesaw pattern in the Z500 anomalies between the subtropics and mid-latitudes of the western North Pacific (WNP) (Fig. 1a). Concurrent with the dipole atmospheric anomaly pattern over the WNP, marked positive Z500 anomalies were observed north of the North American Great Lakes (Fig. 1a). In addition, largescale positive Z500 anomalies appeared in the tropics, suggesting a connection of the WP with tropical SST anomalies, which was later investigated. It is known that the WP may vary in spatial structure due to the differences in the methods and variables utilized to define the WP (Aru et al. 2021). The spatial pattern of the 250-hPa geopotential height anomalies associated with the WP index in this study has been found to share a strong similarity with the spatial pattern of WP defined by Dai and Tan (2019), but with all centres of action shifted slightly eastward (figure not shown). In addition, significant WP-related signals are also existed over the tropical Pacific, which is not evident as shown in Dai and Tan (2019). The discrepancies may be attributed to different timescales focused in Dai and Tan (2019) and our study. In particular, our study emphasized the interannual variability of seasonal-averaged WP with the 9-year high-pass filtered data, whilst Dai and Tan (2019) focused primarily on the sub-seasonal variability. We also examined the WP based on the JRA55 and ERAinterim reanalysis datasets. The winter WP indices derived from the three reanalysis datasets were highly correlated with each other during the overlapping period (Fig. 1b). The correlation coefficient between the WP indices obtained from the JRA55 (ERA-interim) and NCEP-NCAR datasets over 1959-2019 (1979-2019) reached 0.94 (0.89), which was significant at the 99.9% confidence level. In addition, Fig. 1b shows that the amplitude of the winter WP index seems to be larger before the early 1980s than after. In the following text, we only present the results derived from the NCEP-NCAR dataset. The above analysis indicates that the winter WP variability exhibited notable changes in the past. To help understand the characteristics and factors associated with these interdecadal changes in WP variability, we selected the two 25-year epochs with the highest and lowest WP variabilities for a subsequent comparative analysis. The standard deviation of the winter WP was highest during 1950-1974 (corresponding  In particular, the standard deviation of the winter WP index during 1950-1974 was 1.12, which is approximately two times greater than that during 1989-2013 (0.66). Notably, the difference in the standard deviations of the winter WP index between 1950 and 1974 and 1989-2013 was significant at the 95% confidence level according to a two-tailed F-test. Hence, 1950Hence, -1974Hence, and 1989Hence, -2013 were selected to represent high and low epochs of WP variability in the comparative analysis.

Changes in atmospheric and SAT anomalies related to the WP
In this subsection, we first compare the differences in the spatial structure of the WP between 1950-1974 and 1989-2013. Then, we examine the WP-related SAT anomalies during the two periods. Figure 3 shows the composite atmospheric circulation anomalies between the positive and negative phases of the winter WP index during the two studied epochs. The positive and negative phases of the winter WP index were identified according to the application of ± 0.5 standard deviations (Table 2). During both epochs, the positive phase of the WP was associated with a significant increase in zonal winds over the North Pacific at approximately 40° N-55° N, indicative of a northward shift of the eddy-driven jet stream over the northwestern Pacific ( Fig. 3a and b); this result is consistent with the findings of Li and Wettstein (2012) and Wettstein and Wallace (2010).
In addition, significant decreases in 200-hPa zonal winds were observed around the Russian Far East with an eastward extension to Alaska ( Fig. 3a and d). In comparison, the zonal wind anomalies over the North Pacific related to the winter WP index covered larger areas and extended farther westward during 1950-1974 than during 1989-2013. Simultaneously, significant decreases in the 200-hPa zonal wind also occurred to the south of Japan, in the tropical eastern Pacific and in the Great Lakes region during 1950-1974 that could not be observed during 1989-2013. The composited 500-hPa geopotential height anomalies during 1950-1974 exhibited a northwest-southeast-oriented dipole pattern that stretched zonally over the subtropics and mid-and high-latitude of the North Pacific (Fig. 3b). During 1989-2013, a marked dipole pattern also appeared over the western North Pacific, but its intensity was weaker and its latitudinal extension was smaller than that of its counterpart during 1950-1974). In addition, prominent dipole sea level pressure (SLP) and 850-hPa atmospheric anomaly patterns occurred over the western Pacific in the two epochs (Fig. 3c, f), with much stronger amplitudes in the earlier epoch. The above results indicate that larger WP variability tends to correspond to stronger atmospheric anomaly amplitudes over the North Pacific. Studies have demonstrated that WP-related atmospheric anomalies significantly impact SAT variations over Eurasia and North America by modulating atmospheric anomalies (Linkin and Nigam 2008;Lim and Kim 2013;Takaya and Nakamura 2013;Yu and Zhang 2015;Baxter and Nigam 2015;Shi et al. 2019). In the following text, we further examine differences in the WP-related SAT anomalies recorded between the 1950-1974 and 1989-2013 periods. Figure 4a and b display composited difference of wintertime SAT anomalies between the positive and negative phases of the winter WP index during the two studied epochs, respectively. During 1950-1974, significant decreases in the SAT occurred around the Russian Far East; these decreases were attributed to the northerly wind anomalies in this region (Fig. 4a). Under these conditions, anomalous northerly winds brought colder air from higher latitudes, leading to SAT cooling (Fig. 4a). The SATs in northern Mexico and southeastern America also presented significant decreases in this period. In addition, significant SAT warming was seen over central Canada and many parts of south and southeast Asia (Fig. 4a). Marked increases in SAT also appeared over the northern part of South America and a small patch of Alaska (Fig. 4a). Conversely, during 1989-2013, the composited SAT anomalies associated with WP are not that much evident (Fig. 4b). Smaller patches of SAT warming . c, f As in (a, d) but for SLP (shading, unit: hPa) and 850-hPa winds (vectors, unit: m s −1 ). The shaded anomalies are significant at the 95% confidence level. Wind anomalies less than 0.1 m s −1 are not shown Table 2 List of the positive and negative years of the 9-year high-pass-filtered PC1 time series (defined as the WP index) pertinent to the EOF1 of 500-hPa geopotential height anomalies over the western North Pacific during 1949-2019 The positive and negative cases are selected as the years in which the absolute normalized values of the PC time series are larger than 0.5 1950-19741989Positive WPs 1952, 1954, 1958, 1959, 1960, 1964, 1966, 1967, 1969, 19731989, 1992, 1998, 1999, 2012Negative WPs 1950, 1956, 1957, 1962, 1963, 1965, 1968, 1971, 19741991, 1997, 2000, 2005, 2011 anomalies over central Canada and southeast Asia during 1950-1974, we have examined the low-level wind anomalies associated with the winter WP index. The result indicates that, to some extent, the warming both in central Canada and southeast Asia can be partly ascribed to the southerly wind anomalies there, which carrying warmer air northward from lower latitudes (figures not shown). In addition, central Canada and southeast Asia are covered by positive geopotential height anomalies, and the total cloud cover tends to be decreased. Thus, enhanced downward shortwave radiation may also contribute to increase in the surface air temperature there.

Factors affecting the interdecadal variations in the winter WP
In this section, we examine the factors affecting the changes in WP variability. We first examine the role of tropical Pacific SST anomalies. Then, we discuss the role of North Pacific storm track activity.

Role of tropical Pacific SST anomalies
Studies have suggested that the tropical SST anomalies that are related to ENSO significantly contribute to WP variability by triggering an extratropical atmospheric teleconnection (Horel and Wallace 1981; Barnston and Livezey 1987;Koide and Kodera 1999;Wang et al. 2007;Furtado et al. 2012;He et al. 2013;Park et al. 2018;Aru et al. 2021;Yang and Huang 2021). Figure 5 presents the lead-lag correlation coefficients of the wintertime WP index (D(-1)JF(0)) with the Niño3.4 index from the preceding summer (JJA(-1)) to the following summer (JJA (0)). Here, the time notations "− 1" and "0" denote the order of the years. A remarkable difference was found in the WP-ENSO correlation between the two epochs (Fig. 5). During 1950During -1974, the wintertime WP index had a significant positive correlation with the Niño3.4 index from the previous summer (JJA(− 1)) to the following spring (MAM (0)). In contrast, the correlation of the winter WP index with the Niño3.4 index was weak during the 1989-2013 period. Hence, ENSO has a strong (weak) relationship with the WP during periods of strong (weak) WP variability. The above results imply that the observed changes in WP variability may be partly related to changes in the ENSO-WP relation. To confirm the above speculation, we examined the 25-year moving standard deviations of the winter WP index and the 25-year running correlation coefficients of the winter WP index with the Niño3.4 index (Fig. 6a). As is clear in Fig. 6a, the changes in the winter WP-ENSO relation are in good alignment with the changes in WP variability. Figure 6b presents a scatter plot of the 25-year moving standard deviations of the winter WP index versus the 25-year running correlations of the winter WP index with the Niño3.4 index. The two quantities shown in Fig. 6b are highly correlated with each other, with correlation coefficients as high as 0.94. Therefore, the above evidence indicates that changes in the WP-ENSO relationship are closely related to variations in WP variability. Specifically, during period with strong (weak) WP-ENSO connections, the variability in the WP is strong (weak). We only show the anomalies that are significant at the 95% confidence level  1950-1974and 1989. During 1950-1974, the geopotential height anomalies in the upper troposphere (200 hPa) featured distinct negative anomalies over the Russian Far East extending eastward to the west of Alaska, accompanied by large-scale positive anomalies over the tropics (Fig. 7a). The Z500, 850-hPa wind and SLP anomalies (Fig. 7b, c) were similar in this period to those related to the WP (Fig. 1a), with positive Z500, SLP and anticyclonic anomalies observed over the subtropical WNP and anomalies with opposite signs observed over the midhigh-latitude WNP. Quantitatively, the pattern correlation coefficient between the WP-related Z500 anomalies (Fig. 3b) and the ENSO-related Z500 anomalies (Fig. 7b) over region of 25° N-70° N, 120° E-180° E is as high as 0.99. Hence, during this period, the ENSO-related SST anomalies in the tropics strongly impacted the WP and thus contributed to stronger WP variability. The atmospheric anomalies related to the Niño3.4 index over the North Pacific that occurred during 1989-2013 ( Fig. 7e-g) show notable differences from those recorded during 1950-1974 (Fig. 7a-c). In particular, the atmospheric anomalies related to the Niño3.4 index exhibited a PNA-like pattern, with strong negative geopotential height, SLP and cyclonic anomalies over the midlatitude North Pacific and the southeastern United States, along with anomalies with opposite signs over the subtropical northeastern Pacific and northern Canada (Fig. 7e-g). More specifically, we have also calculated the pattern correlation coefficient between the WP-related (Fig. 3e) and the ENSO-related Z500 anomalies (Fig. 7f) over 25° N-70° N, 120° E-180° E. The pattern correlation is only − 0.39 during 1989-2013, which is much less than that during 1950-1974. The above evidence indicates that the atmospheric anomalies generated by tropical ENSO projected more onto the PNA pattern during 1989-2013, in sharp contrast to the WP pattern during 1950-1974. Therefore, the weak contribution of the tropical SST anomalies related to ENSO to the WP may partly explain the weak WP variability recorded during 1989-2013.
The above evidence indicates that ENSO-related atmospheric anomalies showed notable differences between the two periods. A question naturally arose: what leads to changes in the spatial structure of atmospheric anomalies that are related to ENSO? Fig. 7d and h depict the regression patterns of winter SST anomalies upon the winter Niño3.4 index. The spatial patterns of the SST anomalies in the tropics appeared similar during the two periods, with significant SST warming occurring in the tropical central and eastern Pacific and in the tropical Indian Ocean, accompanied by SST cooling in the tropical western Pacific. In addition, the amplitudes of the SST anomalies in the tropical central and eastern Pacific were similar between the two periods. In particular, we examined a scatter plot between the 25-year running correlations of the winter WP index and Niño3.4 index against the 25-year running standard deviations of the Niño3.4 index (not shown). The correlation coefficient between the above two variables was weak, indicating that a change in the ENSO intensity could not explain changes in the spatial structure of ENSO-induced extratropical atmospheric anomalies and thus could not explain changes in the ENSO-WP relation.
Studies have demonstrated that the spatial distributions of atmospheric anomalies induced by tropical SST anomalies are highly sensitive to the zonal locations of SST anomalies in the tropics (Trenberth et al. 1998;Hoerling and Kumar 2002;Barsugli and Sardeshmukh 2002;Zhou et al. 2014;Jo et al. 2015;Soulard et al. 2019;Chen et al. 2021). In particular, the ENSO-induced atmospheric anomalies over the North Pacific shifted eastward (westward) when the ENSO-related SST anomalies in the tropical Pacific were located farther eastward (westward) (Jo et al. 2015;Soulard et al. 2019;Chen et al. 2021). Following Chen et al. (2021), a quantification of the zonal shift of the ENSOrelated SST anomalies in the tropics was defined as the longitude of the 0 °C isotherms of the SST anomalies in the equatorial western Pacific. For accuracy, we interpolated the SST data onto a higher-resolution grid of 0.5° × 0.5°, taking the mean longitude of the 0 °C isotherms between 5°S and 5°N as the index describing the east-west movement of the ENSO-related SST anomaly pattern, termed the "zero line". Figure 8a presents a scatter plot of the zero line of SST anomalies in the tropical Pacific versus the 25-yr moving standard deviations of the WP index. Figure 8b shows a scatter plot of the zero line versus the 25-yr running correlation of the winter WP index with the Niño3.4 index. The correlation coefficients between the two quantities shown in Fig. 8a and b were − 0.50 and − 0.38, respectively, both of which are significant at the 99% confidence These results suggest that changes in the zonal location of tropical SST anomalies related to ENSO may partly contribute to changes in the spatial structure of the atmospheric anomalies induced by ENSO; this result is consistent with previous findings (Jo et al. 2015;Soulard et al. 2019;Chen et al. 2021).

Impact of ENSO on the WP variability in a coupled model simulation
The above analyses indicate that the observed changes in the variability in the WP were partly attributed to changes in the ENSO-WP relationship. In this subsection, we use a long historical simulation of a coupled climate model over 1850-2014 to confirm this conclusion. This long simulation was conducted using the National Aeronautics and Space Administration (NASA) Goddard Institute for Space Studies (GISS) climate model GISS-E2.1-G, which is participating in CMIP6 (Kelley et al. 2020). A more comprehensive description of GISS-E2.1-G can be found in Kelley et al. (2020). The atmospheric variables obtained from the GISS-E2.1-G historical simulation were interpolated to a common horizontal resolution of 2.5° × 2.5° to facilitate comparisons with the observations. The WP index in GISS-E2.1-G was obtained by projecting the model-simulated Z500 anomalies onto the observed EOF1 loading pattern. We first examined the ability of the GISS-E2.1-G model to simulate the spatial patterns of the WP and ENSO. Figure 9a displays the winter 500-hPa geopotential height anomalies obtained by regression on the normalized WP index for the 1850-2014 period. This figure depicts a canonical meridional dipole pattern over the subtropical western North Pacific and mid-to high-latitude North Pacific and the Far East, bearing a close resemblance to the WP shown in the observations (Fig. 1a). In addition, ENSO-related SST anomalies are simulated well by the model, featuring prominent SST warming in the tropical central and eastern Pacific and Indian Ocean, together with SST cooling in the tropical western Pacific (Fig. 9b). Moreover, the GISS-E2.1-G model simulated the spatial distribution of the WP-related SAT anomalies over Eurasia and the North America well, as has been reported in previous studies (Sung et al. 2019(Sung et al. , 2020Aru et al. 2021). Thus, it is reasonable to use the GISS-E2.1-G model to examine changes in the variability in the WP and changes in the WP-ENSO relation. Figure 10a shows the 25-yr moving standard deviations of the WP index and the 25-year running correlation coefficients between the winter WP index and Niño3.4 index obtained with the GISS-E2.1-G model. The change in the variability in the wintertime WP corresponded well with the variation in the WP-ENSO relation, except during the early twentieth century (Fig. 10a). The time series of the 25-year moving standard deviation of the WP index has a correlation coefficient of 0.62 with the time series of the 25-year running correlation between the WP index and the Niño3.4 index, which was significant at the 99% confidence level (Fig. 10b). Hence, the GISS-E2.1-G-derived long-term simulations further confirmed that changes in the winter WP variability are closely related to changes in the WP-ENSO relation. In addition, a statistically significant correlation coefficient of 0.35 was found between the variation in the WP-ENSO relation and the variation in the zero line of the ENSO-related SST anomalies in the tropical Pacific (Fig. 10c). This finding confirmed the observed results in which changes in the wintertime WP-ENSO relation were determined to be partly related to the zonal shifts of ENSOrelated SST anomalies in the tropics.
To further validate the robustness of aforementioned result, longer time period data from the NOAA-CIRES-DOE Twentieth Century Reanalysis (V3) and monthly SST covering 1836-2015 are adopted. We first examine the ability of the NOAA-CIRES-DOE V3 data to reproduce the key features of WP and ENSO. A classic dipole of WP is observed which shares a close resemblance with the WP in Fig. 1a. Meanwhile, the NOAA reanalysis data can reproduce the WP-related surface air temperature anomalies over the Eurasia and North America. In addition, the ENSO-related SST anomalies can also be reproduced with prominent SST warming in the tropical central and eastern Pacific (Please see Fig. S1 in the supplementary material). Hence, it is reasonable to utilize the NOAA-CIRES-DOE Version 3 reanalysis data to detect whether the WP decadal variability is robust over a longer time period. The time series of the 25-year moving standard deviation of the WP index has a correlation coefficient of 0.75 with that of the 25-year running correlation between the WP index and the Niño3.4 index, which is significant at the 99% confidence level (Please see Fig. S2 in the supplementary material). The results confirm that the connection between the interdecadal variations in the WP variability and the WP-ENSO relationship is robust. However, it should be noted that further validation of the relationship between the variation in storm track intensity and the interdecadal WP variability cannot be made due to lack of daily data.

Role of North Pacific storm track activity
Apart from tropical SSTs, studies have suggested that the interaction between the low-frequency mean flow and synoptic-scale eddies is a crucial factor affecting the formation and maintenance of atmospheric circulation patterns associated with the WP (Lau 1988;Chang and Fu 2002;Nakamura et al. 2002;Linkin and Nigam 2008;Tanaka et al. 2016;Wettstein and Wallace 2010). This, therefore, suggests that changes in the wave-mean flow interactions and associated synoptic-scale eddy feedbacks may also contribute to changes in the WP variability. Figure 11a and c show regression maps of storm track anomalies against the winter WP index for the 1950-1974and 1989periods, respectively. During 1950-1974 increases in storm track activity occurred in the area of approximately 40°-60° N over the North Pacific (Fig. 11a), corresponding well with westerly anomalies in the same region (Fig. 3c). Studies have demonstrated that increased (decreased) storm track activity is associated with westerly (easterly) anomalies and is immediately accompanied by positive anticyclonic vorticity forcing to the south and cyclonic vorticity forcing to the north (Lau et al. 1988;Chen et al. 2014), facilitating the maintenance of a WP-like meridional atmospheric dipole pattern (Fig. 3c). Positive storm track anomalies were also observed over the North Pacific during 1989-2013, but these anomalies were less statistically significant (Fig. 11c). We further calculated the geopotential height tendency to quantitatively examine changes in the feedback of the synoptic-scale eddy to the mean flow (e.g., Lau and Holopainen 1984;Lau 1988;Cai et al. 2007). Figure 11b and d show the regression maps of geopotential height tendency at 500 hPa onto the winter WP index during the two studied periods. During 1950During -1974, the spatial structure of the geopotential height tendency anomalies (Fig. 11b) was similar to those of the atmospheric anomalies related to the WP index (Fig. 3b).
In particular, notable negative geopotential height tendency anomalies occurred over the Russian Far East, extending eastward to the Aleutian region, and positive geopotential height tendency anomalies were found over East Asia, extending eastward to the central North Pacific (Fig. 11b).   , during 1950-1974, the feedback of synoptic eddies to the mean flow was conducive to the formation and maintenance of the WP. During 1989-2013, the geopotential height tendency anomalies also exhibited a meridional dipole pattern over the North Pacific but with much weaker amplitudes and less significance compared to those that occurred during 1950-1974. Hence, a stronger wave-mean flow interaction and stronger feedback of synoptic-scale eddies to the mean flow may also partly contribute to the stronger WP variability that occurred during . What contributes to changes in the feedback strength of synoptic-scale eddies? Studies have indicated that the feedback strength of synoptic-scale eddies is related to changes in the mean-state circulation and amplitude of storm track activity (Lau 1988;Chang and Fu 2002;Jin et al. 2006a, b;Kug et al. 2009;Jin 2010). We examined the difference in the mean flow between the 1950-1974 and 1989-2013 periods. The difference was weak over the North Pacific (not shown). This implies that the change in the feedback strength of the synoptic-scale eddies to the mean flow is not related to changes in the mean-state circulation. Then, we further examined the change in the amplitude of the storm track activity. Following previous studies (Lee et al. 2012;Wettstein and Wallace 2010;Ma and Zhang 2018), we used EOF analysis to obtain the dominant mode of storm track activity over the North Pacific (20°-90° N, 120° E-100° W). The long-term trend was removed prior to the EOF analysis. Figure 12a depicts the EOF1 of the winter storm track activity over the North Pacific, and Fig. 12b shows the corresponding PC time series. The EOF1 of the North Pacific storm track activity accounted for 19.62% of the total variance and was characterized by a monopole pattern over the North Pacific. The center of EOF1 corresponded well with the maximum center of the mean winter storm track, as has been reported in previous studies (Lee et al. 2012;Wettstein and Wallace 2010). This suggests that EOF1 represents variations in the wintertime North Pacific storm track intensity. Thus, the corresponding PC1 time series was used to represent the intensity index of the North Pacific storm track. The variability in the PC1 time series was stronger during 1950-1974 than during 1989-2013, indicating that the amplitude of the North Pacific storm track activity was stronger during 1950-1974 than during 1989-2013, corresponding to the stronger feedback of synoptic-scale eddies to the mean flow measured during the earlier period. Figure 13 shows a scatter plot of the 25-year moving standard deviations of the North Pacific storm track intensity index against the 25-year moving standard deviations of the winter WP index. The correlation coefficient between the two variables shown in Fig. 12 reached 0.69, which was significant at the 99.9% confidence level. This result indicates that changes in the intensity of the North Pacific storm track also play a role in modulating changes in the WP intensity.

Conclusion and discussion
This study investigates changes in the variability in the wintertime WP. We found that the variability in the WP showed a significantly weakening trend from the 1950s to the late 1990s. After the early 2000s, the variability in the WP exhibited a slightly increasing trend. The two epochs with the strongest and weakest WP variabilities were selected for a comparative analysis. Figure 14 summarizes the differences in the atmospheric and SAT anomalies related to the WP between the two epochs, as well as the factors responsible for the changes in WP variability.
During the high-variability epoch, a striking meridional atmospheric anomaly pattern occurred over the North Pacific, with significant anticyclonic anomalies over the western and central subtropical regions of the North Pacific and cyclonic anomalies over the Russian Far East and the Bering Strait. The atmospheric circulation anomalies during the low-variability epoch were much weaker and covered a smaller area than those recorded during the high-variability epoch (Fig. 14). Due to the stronger atmospheric anomalies, the WP had a more pronounced impact on the SAT anomalies over most parts of the Eurasian continent and North America during the high-variability epoch than during the low-variability epoch (Fig. 14).
During the high-variability epoch, an explicit link was observed between ENSO and the WP (Fig. 14). In contrast, during the low-variability epoch, the ENSO-induced extratropical atmospheric circulation anomalies were closely aligned with those of the PNA pattern but had a weaker relation with the WP (Fig. 14). The stronger contribution of ENSO during the high-variability epoch partly contributed to the stronger variability in the winter WP. Further analysis suggested that the observed changes in the ENSO-WP relationship were partly related to the zonal shift in ENSOrelated SST anomalies over the tropical Pacific (Fig. 14). The SST anomalies in the tropical Pacific and ENSO-related extratropical North Pacific atmospheric anomalies were located farther westward during the high-variability epoch than during the low-variability epoch, causing a stronger projection on the WP and indicating a closer ENSO-WP connection (Fig. 14). The changes in the variability in the WP were also connected to changes in the intensity of the North Pacific storm track. The stronger intensity of the North Pacific storm track during the high-variability epoch led to stronger wave-mean flow interactions and stronger feedback of synoptic-scale eddies to the mean flow, also contributing to larger variability in the WP during this period (Fig. 14).
This study has revealed that the WP-ENSO relation and the North Pacific storm track intensity are two important factors which may exert large impacts on the WP variability. Previous studies have reported that ENSO could modulate North Pacific storm track variation (Trenberth and Hurrell 1994;Straus and Shukla 1997;Harnik and Chang 2004;Yang et al. 2020). One may ask whether there exists a link between ENSO and the North Pacific storm track intensity. To address this issue, we further examine the connection between the Niño3.4 index and the North Pacific storm track intensity index (which is represented by the PC1 of the North Pacific storm track as described above). The correlation coefficient between the two indices is only 0.15, implicating barely any significant linkage between them.
In recent decades, emerging scientific interests also focused on how ENSO affects the WP (Kodera 1998;Dai and Tan 2016;Tanaka et al. 2016). SST anomalies in the central equatorial Pacific have been suggested to modulate the frequency of WP occurrence (Horel and Wallace 1981;Kodera 1998;Tanaka et al. 2016). A recent study by Dai and Tan (2019) sheds light on how ENSO impacts different types of WP events. Dai and Tan (2019) pointed out that the PNA-preceded WP has a larger amplitude and stronger climate impacts than the PNA-free WP. Simultaneously, ENSO exerts great impacts on the amplitude and frequency of the occurred PNA-preceded WP (Dai and Tan 2019). Our study found that changes in the variability of WP may depend on the WP-ENSO relationship. Accordingly, it is worth investigating influence of the PNA-preceded WP-ENSO relation on the WP variability in future.
Previous studies have indicated that Eurasian snow cover and Arctic Sea ice anomalies (Kodera 1998;Linkin and Nigam 2008;Yuan et al. 2015;Oshika et al. 2014;Nakamura et al. 2015) also play roles in modulating the WP. Due to the limited length of available Eurasian snow cover data and the limited reliability of Arctic sea ice data before the satellite era, the roles of Eurasian snow cover and Arctic sea ice anomalies in modulating interdecadal changes in WP variability remain to be explored. In addition, the present analysis suggests that the zonal shift of SST anomalies in the tropical Pacific related to ENSO may partly contribute to changes in the ENSO-generated atmospheric circulation anomalies over the North Pacific and to changes in the ENSO-WP relation. However, this mechanism cannot fully explain the change in the spatial structure of the ENSOinduced extratropical atmospheric anomalies. The other factors that also contribute to changes in the ENSO-WP relation require further investigation.